1 |
47cdfebe
|
Jim Regetz
|
#!/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))
|