Python Forum

Python Forum > Python Coding > Data Science > monthly composite of LST comes out incorrect, with error "mean of empty slices"
Full Version: monthly composite of LST comes out incorrect, with error "mean of empty slices"
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

ns423

Mar-18-2018, 01:58 PM
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

ns423

Mar-18-2018, 06:40 PM
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\\xxx\\Documents\\xxx\\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

    #Growing season index 
    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'

        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 you 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 like the original files
        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)
[output]lastday = None
        avgtif = None
        #counttif = None
        print "Done!"
        del allarrays
[/output]
Apologies for the mistake, hopefully this is more useful. I am not sure why the "output" does not appear correctly at the bottom of the script. To confirm, the script runs without an error, but with incorrect outputs (which is not what I wrote in the title). thank you.
Python Forum > Python Coding > Data Science > monthly composite of LST comes out incorrect, with error "mean of empty slices"