Revision 5fb982c8
Added by Jim Regetz over 12 years ago
- ID 5fb982c8f6832a514c96d74f37841dc600f77f9e
climate/extra/aggregate-daily-lst.py | ||
---|---|---|
16 | 16 |
# - deal with nightly LST too? |
17 | 17 |
# - deal with more than just first two bits of QC flags? |
18 | 18 |
# e.g.: r.mapcalc 'qc_level2 = (qc>>2) & 0x03' |
19 |
# - record all downloads to a log file? |
|
19 | 20 |
# |
20 | 21 |
# Jim Regetz |
21 | 22 |
# NCEAS |
... | ... | |
36 | 37 |
return date.strftime('%Y.%m.%d') |
37 | 38 |
|
38 | 39 |
# quick function to return list of dirs in wd |
39 |
def list_contents(): |
|
40 |
def list_contents(ftp):
|
|
40 | 41 |
listing = [] |
41 | 42 |
ftp.dir(listing.append) |
42 | 43 |
contents = [item.split()[-1] for item in listing[1:]] |
43 | 44 |
return contents |
44 | 45 |
|
45 |
# create raster of pixelwise N for a list of input rasters |
|
46 |
def calc_nobs(maplist, name, overwrite=False): |
|
47 |
denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % m |
|
48 |
for m in maplist]) |
|
49 |
gs.mapcalc('nobs_%s = %s' % (name, denominator), overwrite=overwrite) |
|
46 |
def download_mod11a1(destdir, tile, start_doy, end_doy, year, ver=5): |
|
47 |
"""Download into destination directory the MODIS 11A1 HDF files |
|
48 |
matching a given tile, year, and day range. If for some unexpected |
|
49 |
reason there are multiple matches for a given day, only the first is |
|
50 |
used. If no matches, the day is skipped with a warning message. |
|
51 |
|
|
52 |
Arguments: |
|
53 |
destdir -- path to local destination directory |
|
54 |
tile -- tile identifier (e.g., 'h08v04') |
|
55 |
start_doy -- integer start of day range (0-366) |
|
56 |
end_doy -- integer end of day range (0-366) |
|
57 |
year -- integer year (>=2000) |
|
58 |
ver -- MODIS version (4 or 5) |
|
59 |
|
|
60 |
Returns list of absolute paths to the downloaded HDF files. |
|
61 |
""" |
|
62 |
# connect to data pool and change to the right directory |
|
63 |
ftp = ftplib.FTP('e4ftl01.cr.usgs.gov') |
|
64 |
ftp.login('anonymous', '') |
|
65 |
ftp.cwd('MOLT/MOD11A1.%03d' % ver) |
|
66 |
# make list of daily directories to traverse |
|
67 |
available_days = list_contents(ftp) |
|
68 |
desired_days = [yj_to_ymd(year, x) for x in range(start_doy, end_doy+1)] |
|
69 |
days_to_get = filter(lambda day: day in desired_days, available_days) |
|
70 |
if len(days_to_get) < len(desired_days): |
|
71 |
missing_days = [day for day in desired_days if day not in days_to_get] |
|
72 |
print 'skipping %d day(s) not found on server:' % len(missing_days) |
|
73 |
print '\n'.join(missing_days) |
|
74 |
# get each tile in turn |
|
75 |
hdfs = list() |
|
76 |
for day in days_to_get: |
|
77 |
ftp.cwd(day) |
|
78 |
files_to_get = [file for file in list_contents(ftp) |
|
79 |
if tile in file and file[-3:]=='hdf'] |
|
80 |
if len(files_to_get)==0: |
|
81 |
# no file found -- print message and move on |
|
82 |
print 'no hdf found on server for tile', tile, 'on', day |
|
83 |
continue |
|
84 |
elif 1<len(files_to_get): |
|
85 |
# multiple files found! -- just use the first... |
|
86 |
print 'multiple hdfs found on server for tile', tile, 'on', day |
|
87 |
file = files_to_get[0] |
|
88 |
local_file = os.path.join(destdir, file) |
|
89 |
ftp.retrbinary('RETR %s' % file, open(local_file, 'wb').write) |
|
90 |
hdfs.append(os.path.abspath(local_file)) |
|
91 |
ftp.cwd('..') |
|
92 |
# politely disconnect |
|
93 |
ftp.quit() |
|
94 |
# return list of downloaded paths |
|
95 |
return hdfs |
|
50 | 96 |
|
51 | 97 |
def get_hdf_paths(hdfdir, tile, start_doy, end_doy, year): |
52 | 98 |
"""Look in supplied directory to find the MODIS 11A1 HDF files |
... | ... | |
75 | 121 |
elif 1<len(files): |
76 | 122 |
# multiple files found! -- just use the first... |
77 | 123 |
print 'multiple hdfs found for tile', tile, 'on', day |
78 |
hdfs.append(files[0])
|
|
124 |
hdfs.append(os.path.abspath(files[0]))
|
|
79 | 125 |
return hdfs |
80 | 126 |
|
81 | 127 |
# create raster of pixelwise means for a list of input rasters |
... | ... | |
87 | 133 |
gs.mapcalc('mean_%s = %s/%s' % (name, numerator, denominator), |
88 | 134 |
overwrite=overwrite) |
89 | 135 |
|
136 |
# create raster of pixelwise N for a list of input rasters |
|
137 |
def calc_nobs(maplist, name, overwrite=False): |
|
138 |
denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % m |
|
139 |
for m in maplist]) |
|
140 |
gs.mapcalc('nobs_%s = %s' % (name, denominator), overwrite=overwrite) |
|
141 |
|
|
90 | 142 |
# ...combined version of the above two functions... |
91 | 143 |
def calc_clim(maplist, name, overwrite=False): |
92 | 144 |
denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % m |
... | ... | |
126 | 178 |
# return name of qc-adjusted LST raster in GRASS |
127 | 179 |
return lstname |
128 | 180 |
|
129 |
#---------------------------------------- |
|
130 |
# set up a (temporary?) GRASS data store
|
|
131 |
#---------------------------------------- |
|
181 |
#---------------------------------------------
|
|
182 |
# test: download then aggregate for one month
|
|
183 |
#---------------------------------------------
|
|
132 | 184 |
|
133 |
## TODO ## |
|
134 |
|
|
135 |
# hack to fix datum??? |
|
136 |
gs.run_command('g.proj', flags='c', |
|
137 |
proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere') |
|
138 |
# proj4='+proj=sinu +R=6371007.181 +nadgrids=@null +wktext') |
|
139 |
|
|
140 |
#------------------------- |
|
141 |
# test download & process |
|
142 |
#------------------------- |
|
143 |
|
|
144 |
tile = 'h09v04' |
|
145 |
year = 2000 |
|
146 |
start_doy = 122 |
|
147 |
end_doy = 152 |
|
148 |
|
|
149 |
# connect to data pool and change to the right directory |
|
150 |
ftp = ftplib.FTP('e4ftl01.cr.usgs.gov') |
|
151 |
ftp.login('anonymous', '') |
|
152 |
ftp.cwd('MOLT/MOD11A1.005') |
|
153 |
|
|
154 |
# make list of daily directories to traverse |
|
155 |
day_dirs = listdirs() |
|
156 |
days_to_get = [yj_to_ymd(year, x) for x in range(start_doy, end_doy+1)] |
|
157 |
if not all [day in day_dirs for day in days_to_get]: |
|
158 |
error |
|
159 |
|
|
160 |
# using v004 data, this was taking 60 seconds for the 8 hdf/xml pairs on |
|
161 |
# atlas [16-May-2012] ... each individual hdf is ~23MB |
|
162 |
hdfs = list() |
|
163 |
for day in days_to_get: |
|
164 |
ftp.cwd(day) |
|
165 |
files_to_get = [file for file in list_contents() |
|
166 |
if tile in file and file[-3:]=='hdf'] |
|
167 |
if len(files_to_get)==0: |
|
168 |
# no file found -- print message and move on |
|
169 |
print 'no hdf found on server for tile', tile, 'on', day |
|
170 |
continue |
|
171 |
elif 1<len(files_to_get): |
|
172 |
# multiple files found! -- just use the first... |
|
173 |
print 'multiple hdfs found on server for tile', tile, 'on', day |
|
174 |
file = files_to_get[0] |
|
175 |
ftp.retrbinary('RETR %s' % file, open(file, 'wb').write) |
|
176 |
hdfs.append(file) |
|
177 |
ftp.cwd('..') |
|
178 |
|
|
179 |
# politely disconnect |
|
180 |
ftp.quit() |
|
181 |
|
|
182 |
|
|
183 |
# |
|
184 |
# test 31-day mean with local files |
|
185 |
# TODO: set up a (temporary?) GRASS database to use for processing? code |
|
186 |
# currently assumes it's being run within an existing GRASS session |
|
187 |
# using an appropriately configured LOCATION... |
|
185 | 188 |
# |
189 |
# note the following trick to fix datum for modis sinu; |
|
190 |
# TODO: check that this works properly...compare to MRT output? |
|
191 |
# gs.run_command('g.proj', flags='c', |
|
192 |
# proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere') |
|
193 |
## proj4='+proj=sinu +R=6371007.181 +nadgrids=@null +wktext') |
|
186 | 194 |
|
187 |
start_doy = 122 |
|
188 |
end_doy = 152 |
|
195 |
tile = 'h09v04' |
|
189 | 196 |
year = 2005 |
190 |
hdfdir = '/home/layers/data/climate/MOD11A1.004-OR-orig' |
|
197 |
start_doy = 1 |
|
198 |
end_doy = 31 |
|
199 |
hdfdir = '.' |
|
200 |
|
|
201 |
# download data |
|
202 |
### [atlas 17-May-2012] Wall time: 111.62 s |
|
203 |
hdfs = download_mod11a1(hdfdir, tile, start_doy, end_doy, year) |
|
191 | 204 |
|
192 |
hdfs = get_hdf_paths(hdfdir, tile, start_doy, end_doy, year) |
|
205 |
# note: now that we've downloaded the files into local directory |
|
206 |
# 'hdfdir', this should yield the same list... |
|
207 |
#hdfs = get_hdf_paths(hdfdir, tile, start_doy, end_doy, year) |
|
193 | 208 |
|
209 |
# generate monthly pixelwise mean & count of high-quality daytime LST |
|
210 |
# values |
|
211 |
### [atlas 17-May-2012] Wall time: 51.18 s |
|
194 | 212 |
gs.os.environ['GRASS_OVERWRITE'] = '1' |
195 |
# load and qc-adjust all daily LSTs, return list of the GRASS map names |
|
196 | 213 |
LST = [load_qc_adjusted_lst(hdf) for hdf in hdfs] |
197 |
# calculate mean and nobs rasters |
|
198 | 214 |
calc_clim(LST, 'LST_%s_%d_%03d_%03d' % (tile, year, start_doy, end_doy)) |
199 | 215 |
# clean up |
200 | 216 |
gs.run_command('g.remove', rast=','.join(LST)) |
Also available in: Unified diff
wrote function to download HDF files for given tile/year/day-range