1
|
#!/usr/bin/python
|
2
|
#
|
3
|
# Early version of assorted Python code to download MODIS 11A1 (daily)
|
4
|
# data associated with specified dates and tiles, then interface with
|
5
|
# GRASS to calculate temporally averaged LST in a way that accounts for
|
6
|
# QC flags.
|
7
|
#
|
8
|
# TODO:
|
9
|
# - functionalize to encapsulate high level procedural steps
|
10
|
# - as specific example of above, write a function that:
|
11
|
# - takes hdf path
|
12
|
# - reads in both lst and qc
|
13
|
# - nulls any lst values with bad qc scores
|
14
|
# - removes qc
|
15
|
# - returns name of stored qc-adjusted lst raster
|
16
|
# - create a 'main' function to allow script to be run as a command
|
17
|
# that can accept arguments?
|
18
|
# - put downloaded data somewhere other than working directory?
|
19
|
# - write code to set up and take down a temporary GRASS location/mapset
|
20
|
# - write code to export climatology rasters to file (geotiff? compressed?)
|
21
|
# - calculate other aggregate stats? (stddev, ...)
|
22
|
# - deal with nightly LST too?
|
23
|
# - deal with more than just first two bits of QC flags?
|
24
|
# e.g.: r.mapcalc 'qc_level2 = (qc>>2) & 0x03'
|
25
|
#
|
26
|
# Jim Regetz
|
27
|
# NCEAS
|
28
|
# Created on 16-May-2012
|
29
|
|
30
|
import os
|
31
|
import datetime
|
32
|
import ftplib
|
33
|
import grass.script as gs
|
34
|
|
35
|
#------------------
|
36
|
# helper functions
|
37
|
#------------------
|
38
|
|
39
|
# quick function to reformat e.g. 200065 to 2000.03.05
|
40
|
def yj_to_ymd(year, doy):
|
41
|
date = datetime.datetime.strptime('%d%d' % (year, doy), '%Y%j')
|
42
|
return date.strftime('%Y.%m.%d')
|
43
|
|
44
|
# quick function to return list of dirs in wd
|
45
|
def list_contents():
|
46
|
listing = []
|
47
|
ftp.dir(listing.append)
|
48
|
contents = [item.split()[-1] for item in listing[1:]]
|
49
|
return contents
|
50
|
|
51
|
# create raster of pixelwise N for a list of input rasters
|
52
|
def calc_nobs(maplist, name, overwrite=False):
|
53
|
denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % m
|
54
|
for m in maplist])
|
55
|
gs.mapcalc('nobs_%s = %s' % (name, denominator), overwrite=overwrite)
|
56
|
|
57
|
# create raster of pixelwise means for a list of input rasters
|
58
|
def calc_mean(maplist, name, overwrite=False):
|
59
|
numerator = '(%s)' % '+'.join(['if(isnull(%s), 0, %s)' % (m, m)
|
60
|
for m in maplist])
|
61
|
denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % m
|
62
|
for m in maplist])
|
63
|
gs.mapcalc('mean_%s = %s/%s' % (name, numerator, denominator),
|
64
|
overwrite=overwrite)
|
65
|
|
66
|
# ...combined version of the above two functions...
|
67
|
def calc_clim(maplist, name, overwrite=False):
|
68
|
denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % m
|
69
|
for m in maplist])
|
70
|
gs.mapcalc('nobs_%s = %s' % (name, denominator), overwrite=overwrite)
|
71
|
numerator = '(%s)' % '+'.join(['if(isnull(%s), 0, %s)' % (m, m)
|
72
|
for m in maplist])
|
73
|
#gs.mapcalc('mean_%s = round(float(%s)/nobs_%s)' % (name, numerator, name),
|
74
|
gs.mapcalc('mean_%s = %s/nobs_%s' % (name, numerator, name),
|
75
|
overwrite=overwrite)
|
76
|
|
77
|
#----------------------------------------
|
78
|
# set up a (temporary?) GRASS data store
|
79
|
#----------------------------------------
|
80
|
|
81
|
## TODO ##
|
82
|
|
83
|
# hack to fix datum???
|
84
|
gs.run_command('g.proj', flags='c',
|
85
|
proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere')
|
86
|
# proj4='+proj=sinu +R=6371007.181 +nadgrids=@null +wktext')
|
87
|
|
88
|
#-------------------------
|
89
|
# test download & process
|
90
|
#-------------------------
|
91
|
|
92
|
tile = 'h09v04'
|
93
|
year = 2000
|
94
|
start_doy = 122
|
95
|
end_doy = 152
|
96
|
|
97
|
# connect to data pool and change to the right directory
|
98
|
ftp = ftplib.FTP('e4ftl01.cr.usgs.gov')
|
99
|
ftp.login('anonymous', '')
|
100
|
ftp.cwd('MOLT/MOD11A1.005')
|
101
|
|
102
|
# make list of daily directories to traverse
|
103
|
day_dirs = listdirs()
|
104
|
days_to_get = [yj_to_ymd(year, x) for x in range(start_doy, end_doy+1)]
|
105
|
if not all [day in day_dirs for day in days_to_get]:
|
106
|
error
|
107
|
|
108
|
# using v004 data, this was taking 60 seconds for the 8 hdf/xml pairs on
|
109
|
# atlas [16-May-2012] ... each individual hdf is ~23MB
|
110
|
hdfs = list()
|
111
|
for day in days_to_get:
|
112
|
ftp.cwd(day)
|
113
|
files_to_get = [file for file in list_contents()
|
114
|
if tile in file and file[-3:]=='hdf']
|
115
|
if len(files_to_get)==0:
|
116
|
# no file found -- print message and move on
|
117
|
print 'no hdf found on server for tile', tile, 'on', day
|
118
|
continue
|
119
|
elif 1<len(files_to_get):
|
120
|
# multiple files found! -- just use the first...
|
121
|
print 'multiple hdfs found on server for tile', tile, 'on', day
|
122
|
file = files_to_get[0]
|
123
|
ftp.retrbinary('RETR %s' % file, open(file, 'wb').write)
|
124
|
hdfs.append(file)
|
125
|
ftp.cwd('..')
|
126
|
|
127
|
# politely disconnect
|
128
|
ftp.quit()
|
129
|
|
130
|
|
131
|
#-----------------------------------------------------------------------
|
132
|
gs.os.environ['GRASS_OVERWRITE'] = '1'
|
133
|
|
134
|
LST = list()
|
135
|
for hdf in hdfs:
|
136
|
hdflabel = '_'.join(hdf.split('.')[1:3])
|
137
|
LST.append('LST_' + hdflabel)
|
138
|
# NOTE: this gives 'datum not recognized' warning, ok to ignore?
|
139
|
# read in lst grid
|
140
|
gs.run_command('r.in.gdal',
|
141
|
input='HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:LST_Day_1km' % hdf,
|
142
|
output=LST[-1])
|
143
|
# read in qc grid
|
144
|
gs.run_command('r.in.gdal',
|
145
|
input = 'HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:QC_Day' % hdf,
|
146
|
output = 'qc_this_day')
|
147
|
gs.run_command('g.region', rast=LST[-1])
|
148
|
# null out any LST cells for which QC is not high [00]
|
149
|
gs.mapcalc('${lst} = if((qc_this_day & 0x03)==0, ${lst}, null())',
|
150
|
lst=LST[-1])
|
151
|
# calculate mean and nobs rasters
|
152
|
calc_clim(LST, 'LST_%s_%d_%03d_%03d' % (tile, year, start_doy, end_doy))
|
153
|
# clean up
|
154
|
gs.run_command('g.remove', rast=','.join(LST))
|
155
|
gs.run_command('g.remove', rast='qc_this_day')
|
156
|
gs.os.environ['GRASS_OVERWRITE'] = '0'
|
157
|
#-----------------------------------------------------------------------
|
158
|
|
159
|
|
160
|
#
|
161
|
# test 31-day mean with local files
|
162
|
#
|
163
|
|
164
|
# Note that some statements below are redundant with (and possibly
|
165
|
# slightly different) from code above -- still need to consolidate
|
166
|
|
167
|
start_doy = 122
|
168
|
end_doy = 152
|
169
|
|
170
|
gs.os.environ['GRASS_OVERWRITE'] = '1'
|
171
|
ordir = '/home/layers/data/climate/MOD11A1.004-OR-orig-2000'
|
172
|
LST = list()
|
173
|
for doy in range(start_doy, end_doy+1):
|
174
|
fileglob = 'MOD11A1.A%d%03d.%s*hdf' % (year, doy, tile)
|
175
|
try:
|
176
|
hdfpath = glob.glob(os.path.join(ordir, fileglob))[0]
|
177
|
except:
|
178
|
print '*** problem finding', fileglob, '***'
|
179
|
continue
|
180
|
hdfname = os.path.basename(hdfpath)
|
181
|
LST.append('LST_' + '_'.join(hdfname.split('.')[1:3]))
|
182
|
# TODO: handle failures here...
|
183
|
gs.run_command('r.in.gdal',
|
184
|
input='HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:LST_Day_1km' % hdfpath,
|
185
|
output=LST[-1])
|
186
|
gs.os.environ['GRASS_OVERWRITE'] = '0'
|
187
|
gs.run_command('g.region', rast=LST[0])
|
188
|
#numerator = '(%s)' % '+'.join(['if(isnull(%s), 0, %s)' % (lst, lst) for lst in LST])
|
189
|
#denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % lst for lst in LST])
|
190
|
#gs.mapcalc('LST_mean_%s = %s / %s' % (tile, numerator, denominator))
|
191
|
calc_mean(LST, 'LST_mean_%s_%d_%03d_%03d' % (tile, year, start_doy, end_doy))
|
192
|
# clean up
|
193
|
gs.run_command('g.remove', rast=','.join(LST))
|