1 |
82b0f65e
|
Benoit Parmentier
|
#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 |
b33bef6a
|
Benoit Parmentier
|
from datetime import datetime
|
12 |
|
|
import calendar
|
13 |
82b0f65e
|
Benoit Parmentier
|
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 |
b33bef6a
|
Benoit Parmentier
|
year_end=2010
|
27 |
82b0f65e
|
Benoit Parmentier
|
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 |
b33bef6a
|
Benoit Parmentier
|
|
35 |
82b0f65e
|
Benoit Parmentier
|
#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 |
b33bef6a
|
Benoit Parmentier
|
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 |
82b0f65e
|
Benoit Parmentier
|
#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 |
b33bef6a
|
Benoit Parmentier
|
|
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 |
82b0f65e
|
Benoit Parmentier
|
day=listfiles_wd2[:]
|
61 |
b33bef6a
|
Benoit Parmentier
|
table_data = numpy.ones((12,2))#declaring an array of size stop+1 and initialized with number 1
|
62 |
82b0f65e
|
Benoit Parmentier
|
|
63 |
b33bef6a
|
Benoit Parmentier
|
i=0;
|
64 |
82b0f65e
|
Benoit Parmentier
|
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 |
b33bef6a
|
Benoit Parmentier
|
day_no.sort()
|
66 |
|
|
|
67 |
|
|
nb_month=12
|
68 |
82b0f65e
|
Benoit Parmentier
|
|
69 |
b33bef6a
|
Benoit Parmentier
|
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 |
82b0f65e
|
Benoit Parmentier
|
|
79 |
b33bef6a
|
Benoit Parmentier
|
table_data = numpy.ones((12,2)) #declaring an array of size stop+1 and initialized with number 1
|
80 |
82b0f65e
|
Benoit Parmentier
|
|
81 |
b33bef6a
|
Benoit Parmentier
|
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 |
82b0f65e
|
Benoit Parmentier
|
f4.write(line0)
|
101 |
b33bef6a
|
Benoit Parmentier
|
for line in list_files[i-1]:
|
102 |
82b0f65e
|
Benoit Parmentier
|
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 |
b33bef6a
|
Benoit Parmentier
|
|
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 |
82b0f65e
|
Benoit Parmentier
|
f5.write(line0)
|
111 |
b33bef6a
|
Benoit Parmentier
|
for line in list_files2[i-1]:
|
112 |
82b0f65e
|
Benoit Parmentier
|
line = line[line.rfind('\\')+1:]
|
113 |
|
|
#line = line[line.rfind('\\')+1:-4]
|
114 |
|
|
line=line.replace('.rst','\n')
|
115 |
|
|
f5.write(line)
|
116 |
b33bef6a
|
Benoit Parmentier
|
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 |
82b0f65e
|
Benoit Parmentier
|
|
127 |
b33bef6a
|
Benoit Parmentier
|
# Calulating the average with TSTAT
|
128 |
|
|
rgf_name = wd1+'list_month_'+str(i)+'.rgf'
|
129 |
|
|
outfilename =wd1+'mean_month_'+str(i)
|
130 |
|
|
|
131 |
82b0f65e
|
Benoit Parmentier
|
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 |
b33bef6a
|
Benoit Parmentier
|
rgf_name2 = wd1+'list_month_valid_obs_'+str(i)+'.rgf'
|
139 |
|
|
outfilename2 =wd1+'mean_month_valid_obs_'+str(i)
|
140 |
82b0f65e
|
Benoit Parmentier
|
|
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 |
b33bef6a
|
Benoit Parmentier
|
parameters = wd1+'mean_month_'+str(i)+'_Mean.rst'+'*'+wd1+'mean_month'+str(i)+'_rescaled.rst'+'*3*'+'0.02'
|
152 |
82b0f65e
|
Benoit Parmentier
|
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 |
b33bef6a
|
Benoit Parmentier
|
i=i+1
|
159 |
|
|
|
160 |
|
|
#outfile1=wd1+'table_data.txt'
|
161 |
|
|
#numpy.savetxt(outfile1, table_data, fmt='%-7.2f')
|
162 |
82b0f65e
|
Benoit Parmentier
|
|
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 |
b33bef6a
|
Benoit Parmentier
|
# '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])
|