Mar-18-2018, 01:58 PM
(This post was last modified: Mar-18-2018, 02:02 PM by sparkz_alot.)
Hello!
I am a new to the Forum, and a complete beginner with Python coding, which I am using extensively for my dissertation. I am trying to make monthly composites from MODIS' daily land surface temperature by using a modified code from Guerreiro (2015). The code runs with no errors, but the outputs created are incorrect, ie. the scale range for the average temperature is not what it should be, compared to manual calculations. For example, the output for June 2000 should be 35 - 8 degrees celsius, but is coming out as 15 - 0.5 degrees celsius. Can anyone see where the problem may be?
from osgeo import gdal
from osgeo.gdalconst import *
import os,fnmatch
import numpy as np
from subprocess import call
#from itertools import izip
from shutil import copyfile
#calculate average between tiff files
#-------Choose the input folder and pattern tiff------#
#folders
inputFolder="C:\\Users\\xxxxx\\Documents\\xxxx\\Dis\\MASK_MODLST_DRY\\MOD11A1_lst"
#----------- Function of development ------------#
def calc_average_LST(allarrays):
#function to calculate the mean between the axis without considering the nodatavalue
meanarray = np.nanmean(allarrays,axis=0)
return meanarray
#--------------- Apply Function -----------------#
#Creates a list with all the tiff files that match the pattern
for root, dirs, files in os.walk(inputFolder):
#print root
#print dirs
#print files
# Index of the growing season
year = root[70:75]
print year
pattern_tif='QCMASK_MOD11A1.005_LST_Day_1km_doy'+year+'*_aid0001.tif'
pattern_avgtif='AVERAGE_MOD11A1.005_LST_Day_1km_doy'+year+'*_aid0001.tif'
count = 0
monlength=len(fnmatch.filter(files,pattern_tif))
for tifname in fnmatch.filter(files,pattern_tif):
count = count + 1
tif = gdal.Open(root+"/"+tifname,gdal.GA_ReadOnly)
arraytmp = [[np.array(tif.GetRasterBand(1).ReadAsArray(),dtype=float)]]
print tifname
if 'allarrays' in locals():
x_exists = 1
else:
x_exists = 0
if x_exists:
teste=np.reshape(arraytmp[0][0],-1)
allarrays=np.append(allarrays,[teste],axis=0)
print np.shape(allarrays)
else:
print np.shape(arraytmp[0][0])
rows = np.shape(arraytmp[0][0])[0]
cols = np.shape(arraytmp[0][0])[1]
print "rows = "+str(rows)
print "cols = "+str(cols)
teste=np.reshape(arraytmp[0][0],-1)
print np.shape(teste)
allarrays=[teste]
print 'SHAPE ALL ARRAYS'
print np.shape(allarrays)
tif = None
if (count == monlength):
lastdaytifname=tifname
if files:
print np.shape(allarrays)
#boolarray=allarrays!=0
#validpointsarray=np.sum(boolarray,axis=0)
meanarray = calc_average_LST(allarrays)
print 'SHAPE MEANARRAY'
print np.shape(meanarray)
monthname = root[65:68]
monname = monthname.lower()
print monname
tifname1 = 'AVERAGE_MOD11A1.005_LST_Day_1km_doy'+year+monname+'_aid0001.tif'
#tifname2 = 'COUNT_MOD11A1.005_LST_Day_1km_doy'+year+monname+'_aid0001.tif'
# Open the file for the last day of the month and copy it in order to get
lastday = gdal.Open(root+"/"+lastdaytifname,gdal.GA_ReadOnly)
call(["gdal_translate","-ot","Float64","-of","GTiff","-co","TILED=YES",root+"/"+lastdaytifname,root+"/"+tifname1])
# Open the file just copied for read/write
avgtif = gdal.Open(root+"/"+tifname1,gdal.GA_Update)
#counttif = gdal.Open(root+"\\"+tifname2[0],gdal.GA_Update)
#Change the 1D array containing the mean to be 2D
meanarray=np.reshape(meanarray,(rows,cols))
print np.shape(meanarray)
#validpointsarray=np.reshape(validpointsarray,(rows,cols))
#Updates values from existing files to match meanarray
avgtif.GetRasterBand(1).WriteArray(meanarray)
#counttif.GetRasterBand(1).WriteArray(validpointsarray)
lastday = None
avgtif = None
#counttif = None
print "Done!"
del allarrays
I am a new to the Forum, and a complete beginner with Python coding, which I am using extensively for my dissertation. I am trying to make monthly composites from MODIS' daily land surface temperature by using a modified code from Guerreiro (2015). The code runs with no errors, but the outputs created are incorrect, ie. the scale range for the average temperature is not what it should be, compared to manual calculations. For example, the output for June 2000 should be 35 - 8 degrees celsius, but is coming out as 15 - 0.5 degrees celsius. Can anyone see where the problem may be?
from osgeo import gdal
from osgeo.gdalconst import *
import os,fnmatch
import numpy as np
from subprocess import call
#from itertools import izip
from shutil import copyfile
#calculate average between tiff files
#-------Choose the input folder and pattern tiff------#
#folders
inputFolder="C:\\Users\\xxxxx\\Documents\\xxxx\\Dis\\MASK_MODLST_DRY\\MOD11A1_lst"
#----------- Function of development ------------#
def calc_average_LST(allarrays):
#function to calculate the mean between the axis without considering the nodatavalue
meanarray = np.nanmean(allarrays,axis=0)
return meanarray
#--------------- Apply Function -----------------#
#Creates a list with all the tiff files that match the pattern
for root, dirs, files in os.walk(inputFolder):
#print root
#print dirs
#print files
# Index of the growing season
year = root[70:75]
print year
pattern_tif='QCMASK_MOD11A1.005_LST_Day_1km_doy'+year+'*_aid0001.tif'
pattern_avgtif='AVERAGE_MOD11A1.005_LST_Day_1km_doy'+year+'*_aid0001.tif'
count = 0
monlength=len(fnmatch.filter(files,pattern_tif))
for tifname in fnmatch.filter(files,pattern_tif):
count = count + 1
tif = gdal.Open(root+"/"+tifname,gdal.GA_ReadOnly)
arraytmp = [[np.array(tif.GetRasterBand(1).ReadAsArray(),dtype=float)]]
print tifname
if 'allarrays' in locals():
x_exists = 1
else:
x_exists = 0
if x_exists:
teste=np.reshape(arraytmp[0][0],-1)
allarrays=np.append(allarrays,[teste],axis=0)
print np.shape(allarrays)
else:
print np.shape(arraytmp[0][0])
rows = np.shape(arraytmp[0][0])[0]
cols = np.shape(arraytmp[0][0])[1]
print "rows = "+str(rows)
print "cols = "+str(cols)
teste=np.reshape(arraytmp[0][0],-1)
print np.shape(teste)
allarrays=[teste]
print 'SHAPE ALL ARRAYS'
print np.shape(allarrays)
tif = None
if (count == monlength):
lastdaytifname=tifname
if files:
print np.shape(allarrays)
#boolarray=allarrays!=0
#validpointsarray=np.sum(boolarray,axis=0)
meanarray = calc_average_LST(allarrays)
print 'SHAPE MEANARRAY'
print np.shape(meanarray)
monthname = root[65:68]
monname = monthname.lower()
print monname
tifname1 = 'AVERAGE_MOD11A1.005_LST_Day_1km_doy'+year+monname+'_aid0001.tif'
#tifname2 = 'COUNT_MOD11A1.005_LST_Day_1km_doy'+year+monname+'_aid0001.tif'
# Open the file for the last day of the month and copy it in order to get
lastday = gdal.Open(root+"/"+lastdaytifname,gdal.GA_ReadOnly)
call(["gdal_translate","-ot","Float64","-of","GTiff","-co","TILED=YES",root+"/"+lastdaytifname,root+"/"+tifname1])
# Open the file just copied for read/write
avgtif = gdal.Open(root+"/"+tifname1,gdal.GA_Update)
#counttif = gdal.Open(root+"\\"+tifname2[0],gdal.GA_Update)
#Change the 1D array containing the mean to be 2D
meanarray=np.reshape(meanarray,(rows,cols))
print np.shape(meanarray)
#validpointsarray=np.reshape(validpointsarray,(rows,cols))
#Updates values from existing files to match meanarray
avgtif.GetRasterBand(1).WriteArray(meanarray)
#counttif.GetRasterBand(1).WriteArray(validpointsarray)
lastday = None
avgtif = None
#counttif = None
print "Done!"
del allarrays