Revision b33bef6a
Added by Benoit Parmentier almost 12 years ago
| climate/research/oregon/modis-lst/Average_QC_MOD11A1_Oregon_case.py | ||
|---|---|---|
| 8 | 8 |
import win32com.client |
| 9 | 9 |
import shutil |
| 10 | 10 |
from datetime import date |
| 11 |
from datetime import datetime |
|
| 12 |
import calendar |
|
| 11 | 13 |
import numpy |
| 12 | 14 |
|
| 13 | 15 |
#variables, input parameters |
| ... | ... | |
| 21 | 23 |
min_nb_obs =1 |
| 22 | 24 |
nb_years =10 |
| 23 | 25 |
year_start=2001 |
| 26 |
year_end=2010 |
|
| 24 | 27 |
bi_sextile_years='2004,2008' |
| 25 | 28 |
|
| 26 | 29 |
t1 = 2 #minimum threshold to consider a pixel as good quality in the LST QC flag level1 |
| ... | ... | |
| 28 | 31 |
|
| 29 | 32 |
#START OF THE SCRIPT |
| 30 | 33 |
api = win32com.client.Dispatch('idrisi32.IdrisiAPIServer') # api is the handle for IDRISI api
|
| 31 |
##wd = api.GetWorkingDir(); #Is variable containing the path of the working directory |
|
| 32 |
##wd = 'H:\Benoit_Backup\Paper3_STA_07202011\\' |
|
| 33 |
|
|
| 34 |
##PART applies the QC mask to the |
|
| 35 |
#Idrisi_filename2='Oregon_2004_366_MOD11A1_Reprojected_QC_Day' |
|
| 36 |
#MODIS QC extraction |
|
| 37 |
|
|
| 38 |
listfiles_wd1 = glob.glob(wd2+'*'+'Reprojected'+'*LST*.rst') #list the all the LST files of interest |
|
| 39 |
listfiles_wd3 = glob.glob(wd2+'*'+'Reprojected'+'*QC*.rst') #list the all the LST files of interest |
|
| 40 |
nb_files = len(listfiles_wd1) |
|
| 41 |
loop_number3=0 |
|
| 42 |
i=0 |
|
| 43 |
## |
|
| 44 |
##for i in range(1,nb_files+1): |
|
| 45 |
## |
|
| 46 |
## loop_number3 = loop_number3 + 1 |
|
| 47 |
## filename1= listfiles_wd1[i-1] |
|
| 48 |
## idrisi_filename1= filename1[filename1.rfind('\\')+1:-4] #Remove the path and the extension
|
|
| 49 |
## |
|
| 50 |
## filename3= listfiles_wd3[i-1] |
|
| 51 |
## idrisi_filename3= filename3[filename3.rfind('\\')+1:-4] #Remove the path and the extension
|
|
| 52 |
## |
|
| 53 |
## #Conversion of MODISQC into quality information: level 1 |
|
| 54 |
## parameters ='2'+'*'+wd2+idrisi_filename3+'.rst'+'*'+wd1+idrisi_filename3+'_L1'+'.rst'+'*1' |
|
| 55 |
## module = 'MODISQC' |
|
| 56 |
## print('Running ' + module + ' module with parameters ' + parameters)
|
|
| 57 |
## success = api.RunModule(module, parameters, True, '', '', '', '', True) #the output of api.runmodule is a value '1' if it is sucessful and '0' if not. |
|
| 58 |
## if not success: |
|
| 59 |
## print('Error running ' + module + '!')
|
|
| 60 |
## |
|
| 61 |
## #Conversion of MODISQC into quality information: level 2 |
|
| 62 |
## parameters ='2'+'*'+wd2+idrisi_filename3+'.rst'+'*'+wd1+idrisi_filename3+'_L2'+'.rst'+'*2' |
|
| 63 |
## module = 'MODISQC' |
|
| 64 |
## print('Running ' + module + ' module with parameters ' + parameters)
|
|
| 65 |
## success = api.RunModule(module, parameters, True, '', '', '', '', True) #the output of api.runmodule is a value '1' if it is sucessful and '0' if not. |
|
| 66 |
## if not success: |
|
| 67 |
## print('Error running ' + module + '!')
|
|
| 68 |
## |
|
| 69 |
## #OVERLAY to eliminate bad pixels: this is using QC level 1 |
|
| 70 |
## |
|
| 71 |
## #RECLASS I*C:\Data\Benoit\Clark_University\Python\dissertationunburnt_pixels_selection\OVERLAY_ID_83_399_144_TEST_BURNT_83_144_399_allocating_bool_Class_83overlay_random.rst*C:\Data\Benoit\Clark_University\Python\dissertationunburnt_pixels_selection\selection_83.rst*3*C:\Data\Benoit\Clark_University\Python\dissertationunburnt_pixels_selection\idrtemp.rcl*1 |
|
| 72 |
## parameters = 'I'+'*'+wd1+idrisi_filename3+'_L1'+'.rst'+'*'+wd1+idrisi_filename3+'_L1'+'_rec_'+str(t1)+'.rst'+'*2'+'*0*0*'+str(t1)+'*1*'+str(t1)+'*5'+'*-9999*1' #option 2 is an avl file with option 1 , count in cells in textfile '*.avl' |
|
| 73 |
## module = 'reclass' |
|
| 74 |
## print('Running ' + module + ' module with parameters ' + parameters)
|
|
| 75 |
## success = api.RunModule(module, parameters, True, '', '', '', '', True) #the output of api.runmodule is a value '1' if it is sucessful and '0' if not. |
|
| 76 |
## if not success: |
|
| 77 |
## print('Error running ' + module + '!')
|
|
| 78 |
## |
|
| 79 |
## #OVERLAY to eliminate bad pixels: this is using QC level 2 |
|
| 80 |
## #RECLASS I*C:\Data\Benoit\Clark_University\Python\dissertationunburnt_pixels_selection\OVERLAY_ID_83_399_144_TEST_BURNT_83_144_399_allocating_bool_Class_83overlay_random.rst*C:\Data\Benoit\Clark_University\Python\dissertationunburnt_pixels_selection\selection_83.rst*3*C:\Data\Benoit\Clark_University\Python\dissertationunburnt_pixels_selection\idrtemp.rcl*1 |
|
| 81 |
## parameters = 'I'+'*'+wd1+idrisi_filename3+'_L2'+'.rst'+'*'+wd1+idrisi_filename3+'_L2'+'_rec_'+str(t2)+'.rst'+'*2'+'*0*0*'+str(t2)+'*1*'+str(t2)+'*5'+'*-9999*1' #option 2 is an avl file with option 1 , count in cells in textfile '*.avl' |
|
| 82 |
## module = 'reclass' |
|
| 83 |
## print('Running ' + module + ' module with parameters ' + parameters)
|
|
| 84 |
## success = api.RunModule(module, parameters, True, '', '', '', '', True) #the output of api.runmodule is a value '1' if it is sucessful and '0' if not. |
|
| 85 |
## if not success: |
|
| 86 |
## print('Error running ' + module + '!')
|
|
| 87 |
## |
|
| 88 |
## #"This file was created by the MODISQC module with the command line:",2*H:\Data\IPLANT_project\outrst_data\Oregon_2004_366_MOD11A1_Reprojected_QC_Day.rst*H:\Data\IPLANT_project\IPLANT_idrisi_project_03062012\test2.rst*2 |
|
| 89 |
## ###OVERLAY MODULE### |
|
| 90 |
## parameters = '3*'+wd1+idrisi_filename3+'_L1'+'_rec_'+str(t1)+'.rst'+'*'+wd1+idrisi_filename3+'_L2'+'_rec_'+str(t2)+'.rst'+'*'+wd1+idrisi_filename3+'_mask_'+str(t1)+'_'+str(t2)+'.rst' |
|
| 91 |
## module = 'overlay' |
|
| 92 |
## print('Running ' + module + ' module with parameters ' + parameters)
|
|
| 93 |
## success = api.RunModule(module, parameters, True, '', '', '', '', True) #the output of api.runmodule is a value '1' if it is sucessful and '0' if not. |
|
| 94 |
## if not success: |
|
| 95 |
## print('Error running ' + module + '!')
|
|
| 96 |
## |
|
| 97 |
## #This file was created by the MODISQC module with the command line:",2*H:\Data\IPLANT_project\outrst_data\Oregon_2004_366_MOD11A1_Reprojected_QC_Day.rst*H:\Data\IPLANT_project\IPLANT_idrisi_project_03062012\test2.rst*2 |
|
| 98 |
## ###OVERLAY MODULE### |
|
| 99 |
## parameters = '3*'+wd1+idrisi_filename3+'_mask_'+str(t1)+'_'+str(t2)+'.rst'+'*'+wd2+idrisi_filename1+'.rst'+'*'+wd1+idrisi_filename1+'_masked_'+str(t1)+'_'+str(t2)+'.rst' |
|
| 100 |
## module = 'overlay' |
|
| 101 |
## print('Running ' + module + ' module with parameters ' + parameters)
|
|
| 102 |
## success = api.RunModule(module, parameters, True, '', '', '', '', True) #the output of api.runmodule is a value '1' if it is sucessful and '0' if not. |
|
| 103 |
## if not success: |
|
| 104 |
## print('Error running ' + module + '!')
|
|
| 105 |
## |
|
| 106 |
## ###Reclass valid values ### |
|
| 107 |
## parameters = 'I'+'*'+wd1+idrisi_filename1+'_masked_'+str(t1)+'_'+str(t2)+'.rst'+'*'+wd1+idrisi_filename1+'_masked_bool_'+str(t1)+'_'+str(t2)+'.rst'+'*2'+'*0*0*'+str(min_range)+'*1*'+str(min_range)+'*'+str(max_range)+'*-9999*1' #option 2 is an avl file with option 1 , count in cells in textfile '*.avl' |
|
| 108 |
## module = 'reclass' |
|
| 109 |
## print('Running ' + module + ' module with parameters ' + parameters)
|
|
| 110 |
## success = api.RunModule(module, parameters, True, '', '', '', '', True) #the output of api.runmodule is a value '1' if it is sucessful and '0' if not. |
|
| 111 |
## if not success: |
|
| 112 |
## print('Error running ' + module + '!')
|
|
| 113 |
## |
|
| 114 |
## |
|
| 34 |
|
|
| 115 | 35 |
#PART 2 Calculating average and number of valid observations |
| 116 | 36 |
|
| 117 | 37 |
listfiles_wd2 = glob.glob(wd1+'*'+'Reprojected'+'*LST*'+'_masked_'+str(t1)+'_'+str(t2)+'.rst') #list the all the LST files of interest |
| 38 |
listfiles_wd2.sort() #sorting the list |
|
| 39 |
listfiles_wd4 = glob.glob(wd1+'*'+'*Reprojected'+'*LST*'+'_masked_bool_'+str(t1)+'_'+str(t2)+'.rst') |
|
| 40 |
listfiles_wd4.sort() #sorting the list |
|
| 41 |
|
|
| 118 | 42 |
#To find the number of valid points used in the average calculation, one must reclass the pixel values in IDRISI, this is inefficient and |
| 119 | 43 |
#can be improved a lot later on. |
| 120 | 44 |
|
| ... | ... | |
| 129 | 53 |
f4.close() |
| 130 | 54 |
|
| 131 | 55 |
#EXTRACTING THE DOY from the files names: Note that this code will only work with a specific format: Oregon_2003_054_MOD11A1_Reprojected_LST_Day_1km.rst |
| 56 |
|
|
| 57 |
date_start= date(year_start,1,1) #This creates a date object for the first day of MODIS LST |
|
| 58 |
date_end=date(year_end,12,31) |
|
| 59 |
|
|
| 132 | 60 |
day=listfiles_wd2[:] |
| 133 |
i=0; |
|
| 134 |
for filename in listfiles_wd2: |
|
| 135 |
#filename = listfiles_wd2[i] #list the files of interest |
|
| 136 |
day[i] = filename[filename.rfind('_MOD')-3:filename.rfind('_MOD')]
|
|
| 137 |
i=i+1 |
|
| 61 |
table_data = numpy.ones((12,2))#declaring an array of size stop+1 and initialized with number 1 |
|
| 138 | 62 |
|
| 63 |
i=0; |
|
| 139 | 64 |
day_no =list(set(day)) #this is the unique set of day in the list. It must be 366 if there are bisextiles years |
| 140 |
day_no.sort() |
|
| 65 |
day_no.sort() |
|
| 66 |
|
|
| 67 |
nb_month=12 |
|
| 141 | 68 |
|
| 69 |
list_files=list() |
|
| 70 |
list_files[:]=[] |
|
| 71 |
list_files2=list() |
|
| 72 |
list_files2[:]=[] |
|
| 73 |
monthlist=list() |
|
| 74 |
monthlist[:]=[] |
|
| 75 |
monthlist2=list() |
|
| 76 |
monthlist2[:]=[] |
|
| 77 |
i=1 |
|
| 142 | 78 |
|
| 143 |
table_data = numpy.ones((len(day_no),2)) #declaring an array of size stop+1 and initialized with number 1
|
|
| 79 |
table_data = numpy.ones((12,2)) #declaring an array of size stop+1 and initialized with number 1
|
|
| 144 | 80 |
|
| 145 |
i=0 |
|
| 146 |
#Grouping files per day from 1 to 366 (or |
|
| 147 |
for day in day_no: |
|
| 148 |
i=i+1 |
|
| 149 |
#listfiles_wd3=glob.glob(wd1+'*20**'+day+'*Reprojected'+'*LST*'_masked_'+str(t1)+'_'+str(t2)+'.rst') #list the files of interest: finding the day of each file |
|
| 150 |
listfiles_wd3=glob.glob(wd1+'*20**'+day+'*Reprojected'+'*LST*'+'_masked_'+str(t1)+'_'+str(t2)+'.rst') |
|
| 151 |
table_data[i-1,0]= day |
|
| 152 |
table_data[i-1,1]= len(listfiles_wd3) #number of available data for a specific date |
|
| 153 |
f4 = open( wd1+'list_day_'+day+'.rgf','w') #This is the output rgf masked |
|
| 154 |
line0 = str(len(listfiles_wd3))+'\n' #This line contains the number of files per day |
|
| 81 |
for i in range(1,nb_month+1): |
|
| 82 |
#i=i+1 |
|
| 83 |
#filename = listfiles_wd2[i] #list the files of interest |
|
| 84 |
monthlist[:]=[] #empty the list |
|
| 85 |
monthlist2[:]=[] #empty the list |
|
| 86 |
j=0 |
|
| 87 |
for filename in listfiles_wd2: |
|
| 88 |
date = filename[filename.rfind('_MOD')-8:filename.rfind('_MOD')]
|
|
| 89 |
d=datetime.strptime(date,'%Y_%j') #This produces a date object from a string extracted from the file name. |
|
| 90 |
month=d.strftime('%m')
|
|
| 91 |
if int(month)==i: |
|
| 92 |
monthlist.append(listfiles_wd2[j]) |
|
| 93 |
monthlist2.append(listfiles_wd4[j]) |
|
| 94 |
j=j+1 |
|
| 95 |
list_files.append(monthlist) #This is a list of a list containing a list for each month |
|
| 96 |
list_files2.append(monthlist2) |
|
| 97 |
#NOW Write RGF |
|
| 98 |
f4 = open( wd1+'list_month_'+str(i)+'.rgf','w') #This is the output rgf masked |
|
| 99 |
line0 = str(len(list_files[i-1]))+'\n' #This line contains the number of files per day |
|
| 155 | 100 |
f4.write(line0) |
| 156 |
for line in listfiles_wd3:
|
|
| 101 |
for line in list_files[i-1]:
|
|
| 157 | 102 |
line = line[line.rfind('\\')+1:]
|
| 158 | 103 |
#line = line[line.rfind('\\')+1:-4]
|
| 159 | 104 |
line=line.replace('.rst','\n')
|
| 160 | 105 |
f4.write(line) |
| 161 | 106 |
f4.close() |
| 162 |
|
|
| 163 |
#listfiles_wd3=glob.glob(wd1+'*20**'+day+'*Reprojected'+'*LST*'_masked_'+str(t1)+'_'+str(t2)+'.rst') #list the files of interest: finding the day of each file |
|
| 164 |
listfiles_wd4=glob.glob(wd1+'*20**'+day+'*Reprojected'+'*LST*'+'_masked_bool_'+str(t1)+'_'+str(t2)+'.rst') |
|
| 165 |
#f4 = open( wd2+'ps_'+inputf3+'_'+str(nb_files)+'.rgf','w') |
|
| 166 |
f5 = open( wd1+'list_day_valid_obs_'+day+'.rgf','w') #This is the output rgf masked |
|
| 167 |
line0 = str(len(listfiles_wd4))+'\n' |
|
| 107 |
|
|
| 108 |
f5 = open( wd1+'list_month_valid_obs_'+str(i)+'.rgf','w') #This is the output rgf masked |
|
| 109 |
line0 = str(len(list_files2[i-1]))+'\n' |
|
| 168 | 110 |
f5.write(line0) |
| 169 |
for line in listfiles_wd4:
|
|
| 111 |
for line in list_files2[i-1]:
|
|
| 170 | 112 |
line = line[line.rfind('\\')+1:]
|
| 171 | 113 |
#line = line[line.rfind('\\')+1:-4]
|
| 172 | 114 |
line=line.replace('.rst','\n')
|
| 173 | 115 |
f5.write(line) |
| 174 |
f5.close() |
|
| 175 |
|
|
| 176 |
rgf_name = wd1+'list_day_'+str(day)+'.rgf' |
|
| 177 |
outfilename =wd1+'mean_day_'+day |
|
| 116 |
f5.close() |
|
| 117 |
|
|
| 118 |
# i=i+1 |
|
| 119 |
|
|
| 120 |
|
|
| 121 |
#Grouping files per month from 1 to 12 |
|
| 122 |
#for day in day_no: |
|
| 123 |
# i=i+1 |
|
| 124 |
#table_data[i-1,0]= month |
|
| 125 |
#table_data[i-1,1]= len(list_files[i-1]) #number of available data for a specific date |
|
| 178 | 126 |
|
| 179 |
# Calulating the average with TSTAT |
|
| 127 |
# Calulating the average with TSTAT |
|
| 128 |
rgf_name = wd1+'list_month_'+str(i)+'.rgf' |
|
| 129 |
outfilename =wd1+'mean_month_'+str(i) |
|
| 130 |
|
|
| 180 | 131 |
parameters ='1'+'*'+rgf_name+'*'+outfilename+'*'+str(min_range)+'*'+str(max_range)+'*'+str(background)+'*'+str(min_nb_obs) |
| 181 | 132 |
module = 'TStats' |
| 182 | 133 |
print('Running ' + module + ' module with parameters ' + parameters)
|
| ... | ... | |
| 184 | 135 |
if not success: |
| 185 | 136 |
print('Error running ' + module + '!')
|
| 186 | 137 |
|
| 187 |
rgf_name2 = wd1+'list_day_valid_obs_'+str(day)+'.rgf'
|
|
| 188 |
outfilename2 =wd1+'mean_day_valid_obs_'+day
|
|
| 138 |
rgf_name2 = wd1+'list_month_valid_obs_'+str(i)+'.rgf'
|
|
| 139 |
outfilename2 =wd1+'mean_month_valid_obs_'+str(i)
|
|
| 189 | 140 |
|
| 190 | 141 |
# Calulating the number of valid observation with TSTAT |
| 191 | 142 |
parameters ='100000'+'*'+rgf_name2+'*'+outfilename2+'*'+str(0)+'*'+str(1)+'*'+str(background)+'*'+str(min_nb_obs) |
| ... | ... | |
| 197 | 148 |
|
| 198 | 149 |
# Rescaling the values with scalar |
| 199 | 150 |
#H:\Data\IPLANT_project\python_working_dir\Oregon_2008_366_MOD11A1_Reprojected_LST_Day_1km_masked_1_3.rst*H:\Data\IPLANT_project\IPLANT_idrisi_project_03062012\test.rst*3*0.02 |
| 200 |
parameters = wd1+'mean_day_'+day+'_Mean.rst'+'*'+wd1+'mean_day'+day+'_rescaled.rst'+'*3*'+'0.02'
|
|
| 151 |
parameters = wd1+'mean_month_'+str(i)+'_Mean.rst'+'*'+wd1+'mean_month'+str(i)+'_rescaled.rst'+'*3*'+'0.02'
|
|
| 201 | 152 |
module = 'scalar' |
| 202 | 153 |
print('Running ' + module + ' module with parameters ' + parameters)
|
| 203 | 154 |
success = api.RunModule(module, parameters, True, '', '', '', '', True) #the output of api.runmodule is a value '1' if it is sucessful and '0' if not. |
| 204 | 155 |
if not success: |
| 205 | 156 |
print('Error running ' + module + '!')
|
| 206 | 157 |
|
| 207 |
outfile1=wd1+'table_data.txt' |
|
| 208 |
numpy.savetxt(outfile1, table_data, fmt='%-7.2f') |
|
| 158 |
i=i+1 |
|
| 159 |
|
|
| 160 |
#outfile1=wd1+'table_data.txt' |
|
| 161 |
#numpy.savetxt(outfile1, table_data, fmt='%-7.2f') |
|
| 209 | 162 |
|
| 210 | 163 |
#REPORT ON THE NUMBER OF MISSING VALUES FOR EACH DATE!!! and each pixel?? |
| 211 | 164 |
#GDALcommand='gdal_translate F:\Data\Paper2_STA_03182012\dem\dem.bil F:\Data\Paper2_STA_03182012\dem_output2.rst -of RST' |
| ... | ... | |
| 232 | 185 |
|
| 233 | 186 |
# today = datetime.datetime.now() |
| 234 | 187 |
# today.strftime('%j') Returns the day of the year in a string.
|
| 235 |
# '083' |
|
| 188 |
# '083' |
|
| 189 |
# |
|
| 190 |
|
|
| 191 |
#THIS CAN BE USED TO CREATE A SEQUENCE |
|
| 192 |
##def date_range(start, end): |
|
| 193 |
## r = (end+datetime.timedelta(days=1)-start).days |
|
| 194 |
## return [start+datetime.timedelta(days=i) for i in range(r)] |
|
| 195 |
## |
|
| 196 |
##start = datetime.date(2007,01,01) |
|
| 197 |
##end = datetime.date(2008,02,01) |
|
| 198 |
##dateList = date_range(start, end) |
|
| 199 |
##print '\n'.join([str(date) for date in dateList]) |
|
Also available in: Unified diff
LST climatology, using IDRIS API modification using datetime module, OR task#375 and task#416