Python Forum
monthly composite of LST comes out incorrect, with error "mean of empty slices"
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
monthly composite of LST comes out incorrect, with error "mean of empty slices"
#1
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
Reply
#2
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.
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  How to split monthly data smoothly into days? AlekseyPython 1 1,512 Jan-27-2022, 10:56 AM
Last Post: Larz60+
  selecting only october to march from monthly data u Staph 1 1,628 Jul-14-2019, 10:15 AM
Last Post: Larz60+
  selecting customized seasons from monthly time series Staph 3 2,814 Jul-14-2019, 09:40 AM
Last Post: scidam
  error: INCLUDE environment variable is empty xwchen 4 7,002 Mar-22-2018, 03:40 AM
Last Post: xwchen
  TypeError: list indices must be integers or slices, not str TheDovah7 4 13,345 Mar-06-2018, 03:00 PM
Last Post: TheDovah7
  Modulus operator giving incorrect result Chaos 1 3,014 Jun-16-2017, 12:20 PM
Last Post: Ofnuts

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020