1
|
#This script was written by Benoit Parmentier on March 27, 2012.
|
2
|
#This script computes the daily average for MOD11A1 data. The code takes into account the QC and computes the number of valid observation.
|
3
|
#Note that this code will only work with a specific format: Oregon_2003_054_MOD11A1_Reprojected_LST_Day_1km.rst
|
4
|
#This code should be reviewed to use the datetime package in python. As of now, the code relies on the IDRISI software's API.
|
5
|
|
6
|
import os, glob, sys
|
7
|
from osgeo import gdal, gdalconst
|
8
|
import win32com.client
|
9
|
import shutil
|
10
|
from datetime import date
|
11
|
from datetime import datetime
|
12
|
import calendar
|
13
|
import numpy
|
14
|
|
15
|
#variables, input parameters
|
16
|
|
17
|
wd1='H:\Data\IPLANT_project\python_working_dir3\\' #output folder
|
18
|
wd2='H:\Data\IPLANT_project\outrst_data\\' #input folder containg the QC and LST files
|
19
|
|
20
|
min_range =7500 #This value can be changed, it correspond to the valid range in LST MOD11A1
|
21
|
max_range =65535
|
22
|
background =0
|
23
|
min_nb_obs =1
|
24
|
nb_years =10
|
25
|
year_start=2001
|
26
|
year_end=2010
|
27
|
bi_sextile_years='2004,2008'
|
28
|
|
29
|
t1 = 2 #minimum threshold to consider a pixel as good quality in the LST QC flag level1
|
30
|
t2 = 3 #minimum threshold to consider a pixel as good quality in the LST QC flag level2
|
31
|
|
32
|
#START OF THE SCRIPT
|
33
|
api = win32com.client.Dispatch('idrisi32.IdrisiAPIServer') # api is the handle for IDRISI api
|
34
|
|
35
|
#PART 2 Calculating average and number of valid observations
|
36
|
|
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
|
|
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
|
43
|
#can be improved a lot later on.
|
44
|
|
45
|
f4 = open( wd1+'output_list'+'.rgf','w')
|
46
|
line0 = str(len(listfiles_wd2))+'\n'
|
47
|
f4.write(line0)
|
48
|
for line in listfiles_wd2:
|
49
|
line = line[line.rfind('\\')+1:]
|
50
|
#line = line[line.rfind('\\')+1:-4]
|
51
|
line=line.replace('.rst','\n')
|
52
|
f4.write(line)
|
53
|
f4.close()
|
54
|
|
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
|
|
60
|
day=listfiles_wd2[:]
|
61
|
table_data = numpy.ones((12,2))#declaring an array of size stop+1 and initialized with number 1
|
62
|
|
63
|
i=0;
|
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
|
65
|
day_no.sort()
|
66
|
|
67
|
nb_month=12
|
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
|
78
|
|
79
|
table_data = numpy.ones((12,2)) #declaring an array of size stop+1 and initialized with number 1
|
80
|
|
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
|
100
|
f4.write(line0)
|
101
|
for line in list_files[i-1]:
|
102
|
line = line[line.rfind('\\')+1:]
|
103
|
#line = line[line.rfind('\\')+1:-4]
|
104
|
line=line.replace('.rst','\n')
|
105
|
f4.write(line)
|
106
|
f4.close()
|
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'
|
110
|
f5.write(line0)
|
111
|
for line in list_files2[i-1]:
|
112
|
line = line[line.rfind('\\')+1:]
|
113
|
#line = line[line.rfind('\\')+1:-4]
|
114
|
line=line.replace('.rst','\n')
|
115
|
f5.write(line)
|
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
|
126
|
|
127
|
# Calulating the average with TSTAT
|
128
|
rgf_name = wd1+'list_month_'+str(i)+'.rgf'
|
129
|
outfilename =wd1+'mean_month_'+str(i)
|
130
|
|
131
|
parameters ='1'+'*'+rgf_name+'*'+outfilename+'*'+str(min_range)+'*'+str(max_range)+'*'+str(background)+'*'+str(min_nb_obs)
|
132
|
module = 'TStats'
|
133
|
print('Running ' + module + ' module with parameters ' + parameters)
|
134
|
success = api.RunModule(module, parameters, True, '', '', '', '', True) #the output of api.runmodule is a value '1' if it is sucessful and '0' if not.
|
135
|
if not success:
|
136
|
print('Error running ' + module + '!')
|
137
|
|
138
|
rgf_name2 = wd1+'list_month_valid_obs_'+str(i)+'.rgf'
|
139
|
outfilename2 =wd1+'mean_month_valid_obs_'+str(i)
|
140
|
|
141
|
# Calulating the number of valid observation with TSTAT
|
142
|
parameters ='100000'+'*'+rgf_name2+'*'+outfilename2+'*'+str(0)+'*'+str(1)+'*'+str(background)+'*'+str(min_nb_obs)
|
143
|
module = 'TStats'
|
144
|
print('Running ' + module + ' module with parameters ' + parameters)
|
145
|
success = api.RunModule(module, parameters, True, '', '', '', '', True) #the output of api.runmodule is a value '1' if it is sucessful and '0' if not.
|
146
|
if not success:
|
147
|
print('Error running ' + module + '!')
|
148
|
|
149
|
# Rescaling the values with scalar
|
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
|
151
|
parameters = wd1+'mean_month_'+str(i)+'_Mean.rst'+'*'+wd1+'mean_month'+str(i)+'_rescaled.rst'+'*3*'+'0.02'
|
152
|
module = 'scalar'
|
153
|
print('Running ' + module + ' module with parameters ' + parameters)
|
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.
|
155
|
if not success:
|
156
|
print('Error running ' + module + '!')
|
157
|
|
158
|
i=i+1
|
159
|
|
160
|
#outfile1=wd1+'table_data.txt'
|
161
|
#numpy.savetxt(outfile1, table_data, fmt='%-7.2f')
|
162
|
|
163
|
#REPORT ON THE NUMBER OF MISSING VALUES FOR EACH DATE!!! and each pixel??
|
164
|
#GDALcommand='gdal_translate F:\Data\Paper2_STA_03182012\dem\dem.bil F:\Data\Paper2_STA_03182012\dem_output2.rst -of RST'
|
165
|
#output = os.popen(GDALcommand).read()
|
166
|
|
167
|
#print('END OF SCRIPT')
|
168
|
### ==================================== #
|
169
|
### Load GDAL Drivers:
|
170
|
##dst_frmts = 'rst' #data format driver for gdal
|
171
|
##dst_driver = gdal.GetDriverByName( dst_frmts )
|
172
|
##
|
173
|
##if dst_driver == None:
|
174
|
## Usage()
|
175
|
##
|
176
|
###src_ds = gdal.Open( inputf1 )
|
177
|
##src_ds = gdal.Open( 'G:\\Benoit_Backup\\STA_PAPER2_02262011\\test5_ID.rst' )
|
178
|
##src_ds = gdal.Open('G:\Benoit_Backup\STA_PAPER2_02262011\GAP_ALB_07182011_filled_t3\\GAP_ALB_07182011_filled_t3_GAP_ALBEDO_16D_MASKED2_FILL4_6_MEAN.rst')
|
179
|
##src_rb = src_ds.GetRasterBand( 1 );
|
180
|
##
|
181
|
##x_size = src_rb.XSize #nb_cols
|
182
|
##y_size = src_rb.YSize #nb_rows
|
183
|
##
|
184
|
##src_rb.ReadAsArray(0,6,None,None,24,None)
|
185
|
|
186
|
# today = datetime.datetime.now()
|
187
|
# today.strftime('%j') Returns the day of the year in a string.
|
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])
|