Revision 04d1bd09
Added by Benoit Parmentier over 11 years ago
climate/research/oregon/modis-lst/LST_downloading_and_climatology.py | ||
---|---|---|
1 |
{\rtf1\ansi\ansicpg1252\cocoartf1138\cocoasubrtf320 |
|
2 |
{\fonttbl\f0\fnil\fcharset0 Menlo-Regular;} |
|
3 |
{\colortbl;\red255\green255\blue255;\red0\green116\blue0;\red170\green13\blue145;\red196\green26\blue22; |
|
4 |
\red28\green0\blue207;\red63\green105\blue30;} |
|
5 |
\margl1440\margr1440\vieww28360\viewh15140\viewkind0 |
|
6 |
\deftab560 |
|
7 |
\pard\tx560\pardeftab560\pardirnatural |
|
8 |
|
|
9 |
\f0\fs28 \cf2 \CocoaLigature0 #!/usr/bin/python\cf0 \ |
|
10 |
\cf2 #\cf0 \ |
|
11 |
\cf2 # Early version of assorted Python code to download MODIS 11A1 (daily)\cf0 \ |
|
12 |
\cf2 # data associated with specified dates and tiles, then interface with\cf0 \ |
|
13 |
\cf2 # GRASS to calculate temporally averaged LST in a way that accounts for\cf0 \ |
|
14 |
\cf2 # QC flags.\cf0 \ |
|
15 |
\cf2 #\cf0 \ |
|
16 |
\cf2 # TODO:\cf0 \ |
|
17 |
\cf2 # - functionalize to encapsulate high level procedural steps\cf0 \ |
|
18 |
\cf2 # - create a 'main' function to allow script to be run as a command\cf0 \ |
|
19 |
\cf2 # that can accept arguments?\cf0 \ |
|
20 |
\cf2 # - put downloaded data somewhere other than working directory?\cf0 \ |
|
21 |
\cf2 # - write code to set up and take down a temporary GRASS location/mapset\cf0 \ |
|
22 |
\cf2 # - write code to export climatology rasters to file (geotiff? compressed?)\cf0 \ |
|
23 |
\cf2 # - calculate other aggregate stats? (stddev, ...)\cf0 \ |
|
24 |
\cf2 # - deal with nightly LST too?\cf0 \ |
|
25 |
\cf2 # - deal with more than just first two bits of QC flags?\cf0 \ |
|
26 |
\cf2 # e.g.: r.mapcalc 'qc_level2 = (qc>>2) & 0x03' \ |
|
27 |
# 0x03?? that means 0000_0011 and >>2 means shit twice on the right for \cf0 \ |
|
28 |
\cf2 # - record all downloads to a log file?\cf0 \ |
|
29 |
\cf2 #\cf0 \ |
|
30 |
\cf2 # Jim Regetz\cf0 \ |
|
31 |
\cf2 # NCEAS\cf0 \ |
|
32 |
\cf2 # Created on 16-May-2012\ |
|
33 |
# Modified by Benoit Parmentier\ |
|
34 |
# NCEAS on 18-October-2012\ |
|
35 |
# \cf0 \ |
|
36 |
\ |
|
37 |
\cf3 import\cf0 os, glob\ |
|
38 |
\cf3 import\cf0 datetime, calendar\ |
|
39 |
\cf3 import\cf0 ftplib\ |
|
40 |
\cf3 import\cf0 grass.script \cf3 as\cf0 gs\ |
|
41 |
\ |
|
42 |
\cf2 #------------------\cf0 \ |
|
43 |
\cf2 # helper functions\cf0 \ |
|
44 |
\cf2 #------------------\cf0 \ |
|
45 |
\ |
|
46 |
\cf3 def\cf0 yj_to_ymd(year, doy):\ |
|
47 |
\cf4 """Return date as e.g. '2000.03.05' based on specified year and\ |
|
48 |
numeric day-of-year (doy) """\cf0 \ |
|
49 |
date = datetime.datetime.strptime(\cf5 '%d%03d'\cf0 % (year, doy), \cf5 '%Y%j'\cf0 )\ |
|
50 |
\cf3 return\cf0 date.strftime(\cf5 '%Y.%m.%d'\cf0 )\ |
|
51 |
\ |
|
52 |
\cf3 def\cf0 get_doy_range(year, month):\ |
|
53 |
\cf4 """Determine starting and ending numeric day-of-year (doy)\ |
|
54 |
asscociated with the specified month and year.\ |
|
55 |
\ |
|
56 |
Arguments:\ |
|
57 |
year -- four-digit integer year\ |
|
58 |
month -- integer month (1-12)\ |
|
59 |
\ |
|
60 |
Returns tuple of start and end doy for that month/year.\ |
|
61 |
"""\cf0 \ |
|
62 |
last_day_of_month = calendar.monthrange(year, month)[\cf5 1\cf0 ]\ |
|
63 |
start_doy = int(datetime.datetime.strptime(\cf5 '%d.%02d.%02d'\cf0 % (year,\ |
|
64 |
month, \cf5 1\cf0 ), \cf5 '%Y.%m.%d'\cf0 ).strftime(\cf5 '%j'\cf0 ))\ |
|
65 |
end_doy = int(datetime.datetime.strptime(\cf5 '%d.%02d.%02d'\cf0 % (year,\ |
|
66 |
month, last_day_of_month), \cf5 '%Y.%m.%d'\cf0 ).strftime(\cf5 '%j'\cf0 ))\ |
|
67 |
\cf3 return\cf0 (int(start_doy), int(end_doy))\ |
|
68 |
\ |
|
69 |
\cf2 # quick function to return list of dirs in wd\cf0 \ |
|
70 |
\cf3 def\cf0 list_contents(ftp):\ |
|
71 |
\cf4 """Parse ftp directory listing into list of names of the files\ |
|
72 |
and/or directories contained therein. May or may not be robust in\ |
|
73 |
general, but seems to work fine for LP DAAP ftp server."""\cf0 \ |
|
74 |
listing = []\ |
|
75 |
ftp.dir(listing.append)\ |
|
76 |
contents = [item.split()[-\cf5 1\cf0 ] \cf3 for\cf0 item \cf3 in\cf0 listing[\cf5 1\cf0 :]]\ |
|
77 |
\cf3 return\cf0 contents\ |
|
78 |
\ |
|
79 |
\cf3 def\cf0 download_mod11a1(destdir, tile, start_doy, end_doy, year, ver=\cf5 5\cf0 ):\ |
|
80 |
\cf4 """Download into destination directory the MODIS 11A1 HDF files\ |
|
81 |
matching a given tile, year, and day range. If for some unexpected\ |
|
82 |
reason there are multiple matches for a given day, only the first is\ |
|
83 |
used. If no matches, the day is skipped with a warning message.\ |
|
84 |
\ |
|
85 |
Arguments:\ |
|
86 |
destdir -- path to local destination directory\ |
|
87 |
tile -- tile identifier (e.g., 'h08v04')\ |
|
88 |
start_doy -- integer start of day range (0-366)\ |
|
89 |
end_doy -- integer end of day range (0-366)\ |
|
90 |
year -- integer year (>=2000)\ |
|
91 |
ver -- MODIS version (4 or 5)\ |
|
92 |
\ |
|
93 |
Returns list of absolute paths to the downloaded HDF files.\ |
|
94 |
"""\cf0 \ |
|
95 |
\cf2 # connect to data pool and change to the right directory\cf0 \ |
|
96 |
ftp = ftplib.FTP(\cf5 'e4ftl01.cr.usgs.gov'\cf0 )\ |
|
97 |
ftp.login(\cf5 'anonymous'\cf0 , \cf5 ''\cf0 )\ |
|
98 |
ftp.cwd(\cf5 'MOLT/MOD11A1.%03d'\cf0 % ver)\ |
|
99 |
\cf2 # make list of daily directories to traverse\cf0 \ |
|
100 |
available_days = list_contents(ftp)\ |
|
101 |
desired_days = [yj_to_ymd(year, x) \cf3 for\cf0 x \cf3 in\cf0 range(start_doy, end_doy+\cf5 1\cf0 )]\ |
|
102 |
days_to_get = filter(\cf3 lambda\cf0 day: day \cf3 in\cf0 desired_days, available_days)\ |
|
103 |
\cf3 if\cf0 len(days_to_get) < len(desired_days):\ |
|
104 |
missing_days = [day \cf3 for\cf0 day \cf3 in\cf0 desired_days \cf3 if\cf0 day \cf3 not\cf0 \cf3 in\cf0 days_to_get]\ |
|
105 |
\cf3 print\cf0 \cf5 'skipping %d day(s) not found on server:'\cf0 % len(missing_days)\ |
|
106 |
\cf3 print\cf0 \cf5 '\\n'\cf0 .join(missing_days)\ |
|
107 |
\cf2 # get each tile in turn\cf0 \ |
|
108 |
hdfs = list()\ |
|
109 |
\cf3 for\cf0 day \cf3 in\cf0 days_to_get:\ |
|
110 |
ftp.cwd(day)\ |
|
111 |
files_to_get = [file \cf3 for\cf0 file \cf3 in\cf0 list_contents(ftp)\ |
|
112 |
\cf3 if\cf0 tile \cf3 in\cf0 file \cf3 and\cf0 file[-\cf5 3\cf0 :]==\cf5 'hdf'\cf0 ]\ |
|
113 |
\cf3 if\cf0 len(files_to_get)==\cf5 0\cf0 :\ |
|
114 |
\cf2 # no file found -- print message and move on\cf0 \ |
|
115 |
\cf3 print\cf0 \cf5 'no hdf found on server for tile'\cf0 , tile, \cf5 'on'\cf0 , day\ |
|
116 |
ftp.cwd(\cf5 '..'\cf0 ) \cf2 #added by Benoit\cf0 \ |
|
117 |
\cf3 continue\cf0 \ |
|
118 |
\cf3 elif\cf0 \cf5 1\cf0 <len(files_to_get):\ |
|
119 |
\cf2 # multiple files found! -- just use the first...\cf0 \ |
|
120 |
\cf3 print\cf0 \cf5 'multiple hdfs found on server for tile'\cf0 , tile, \cf5 'on'\cf0 , day\ |
|
121 |
file = files_to_get[\cf5 0\cf0 ]\ |
|
122 |
local_file = os.path.join(destdir, file)\ |
|
123 |
ftp.retrbinary(\cf5 'RETR %s'\cf0 % file, open(local_file, \cf5 'wb'\cf0 ).write)\ |
|
124 |
hdfs.append(os.path.abspath(local_file))\ |
|
125 |
ftp.cwd(\cf5 '..'\cf0 )\ |
|
126 |
\cf2 # politely disconnect\cf0 \ |
|
127 |
ftp.quit()\ |
|
128 |
\cf2 # return list of downloaded paths\cf0 \ |
|
129 |
\cf3 return\cf0 hdfs\ |
|
130 |
\ |
|
131 |
\cf3 def\cf0 get_hdf_paths(hdfdir, tile, start_doy, end_doy, year):\ |
|
132 |
\cf4 """Look in supplied directory to find the MODIS 11A1 HDF files\ |
|
133 |
matching a given tile, year, and day range. If for some unexpected\ |
|
134 |
reason there are multiple matches for a given day, only the first is\ |
|
135 |
used. If no matches, the day is skipped with a warning message.\ |
|
136 |
\ |
|
137 |
Arguments:\ |
|
138 |
hdfdir -- path to directory containing the desired HDFs\ |
|
139 |
tile -- tile identifier (e.g., 'h08v04')\ |
|
140 |
start_doy -- integer start of day range (0-366)\ |
|
141 |
end_doy -- integer end of day range (0-366)\ |
|
142 |
year -- integer year (>=2000)\ |
|
143 |
\ |
|
144 |
Returns list of absolute paths to the located HDF files.\ |
|
145 |
"""\cf0 \ |
|
146 |
hdfs = list()\ |
|
147 |
\cf3 for\cf0 doy \cf3 in\cf0 range(start_doy, end_doy+\cf5 1\cf0 ):\ |
|
148 |
fileglob = \cf5 'MOD11A1.A%d%03d.%s*hdf'\cf0 % (year, doy, tile)\ |
|
149 |
pathglob = os.path.join(hdfdir, fileglob)\ |
|
150 |
files = glob.glob(pathglob)\ |
|
151 |
\cf3 if\cf0 len(files)==\cf5 0\cf0 :\ |
|
152 |
\cf2 # no file found -- print message and move on\cf0 \ |
|
153 |
\cf3 print\cf0 \cf5 'cannot access %s: no such file'\cf0 % pathglob\ |
|
154 |
\cf3 continue\cf0 \ |
|
155 |
\cf3 elif\cf0 \cf5 1\cf0 <len(files):\ |
|
156 |
\cf2 # multiple files found! -- just use the first...\cf0 \ |
|
157 |
\cf3 print\cf0 \cf5 'multiple hdfs found for tile'\cf0 , tile, \cf5 'on'\cf0 , day\ |
|
158 |
hdfs.append(os.path.abspath(files[\cf5 0\cf0 ]))\ |
|
159 |
\cf3 return\cf0 hdfs\ |
|
160 |
\ |
|
161 |
\cf3 def\cf0 calc_clim(maplist, name, overwrite=False):\ |
|
162 |
\cf4 """Generate some climatalogies in GRASS based on the input list of\ |
|
163 |
maps. As usual, current GRASS region settings apply. Produces the\ |
|
164 |
following output rasters:\ |
|
165 |
* nobs: count of number of (non-null) values over the input maps\ |
|
166 |
* mean: arithmetic mean of (non-null) values over the input maps\ |
|
167 |
\ |
|
168 |
Arguments:\ |
|
169 |
maplist -- list of names of GRASS maps to be aggregated\ |
|
170 |
name -- template (suffix) for GRASS output map names\ |
|
171 |
\ |
|
172 |
Returns list of names of the output maps created in GRASS.\ |
|
173 |
"""\cf0 \ |
|
174 |
denominator = \cf5 '(%s)'\cf0 % \cf5 '+'\cf0 .join([\cf5 'if(!isnull(%s))'\cf0 % m\ |
|
175 |
\cf3 for\cf0 m \cf3 in\cf0 maplist])\ |
|
176 |
gs.mapcalc(\cf5 'nobs_%s = %s'\cf0 % (name, denominator), overwrite=overwrite)\ |
|
177 |
numerator = \cf5 '(%s)'\cf0 % \cf5 '+'\cf0 .join([\cf5 'if(isnull(%s), 0, %s)'\cf0 % (m, m)\ |
|
178 |
\cf3 for\cf0 m \cf3 in\cf0 maplist])\ |
|
179 |
gs.mapcalc(\cf5 'mean_%s = round(float(%s)/nobs_%s)'\cf0 % (name, numerator, name),\ |
|
180 |
overwrite=overwrite)\ |
|
181 |
\cf3 return\cf0 [\cf5 '%s_%s'\cf0 % (s, name) \cf3 for\cf0 s \cf3 in\cf0 [\cf5 'nobs'\cf0 , \cf5 'mean'\cf0 ]]\ |
|
182 |
\ |
|
183 |
\cf3 def\cf0 load_qc_adjusted_lst(hdf):\ |
|
184 |
\cf4 """Load LST_Day_1km into GRASS from the specified hdf file, and\ |
|
185 |
nullify any cells for which QA flags indicate anything other than\ |
|
186 |
high quality.\ |
|
187 |
\ |
|
188 |
Argument:\ |
|
189 |
hdf -- local path to the 11A1 HDF file\ |
|
190 |
\ |
|
191 |
Returns the name of the QC-adjusted LST map created in GRASS.\ |
|
192 |
"""\cf0 \ |
|
193 |
lstname = \cf5 'LST_'\cf0 + \cf5 '_'\cf0 .join(os.path.basename(hdf).split(\cf5 '.'\cf0 )[\cf5 1\cf0 :\cf5 3\cf0 ]) \ |
|
194 |
\cf2 # read in lst grid\cf0 \ |
|
195 |
gs.run_command(\cf5 'r.in.gdal'\cf0 ,\ |
|
196 |
input=\cf5 'HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:LST_Day_1km'\cf0 % hdf,\ |
|
197 |
output=lstname) \cf2 # No new location created since it has already been done with r.in.gdal before!!!\cf0 \ |
|
198 |
\ |
|
199 |
\cf2 # read in qc grid\cf0 \ |
|
200 |
gs.run_command(\cf5 'r.in.gdal'\cf0 ,\ |
|
201 |
input = \cf5 'HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:QC_Day'\cf0 % hdf,\ |
|
202 |
output = \cf5 'qc_this_day'\cf0 )\ |
|
203 |
gs.run_command(\cf5 'g.region'\cf0 , rast=lstname)\ |
|
204 |
\cf2 # null out any LST cells for which QC is not high [00] #THIS WHERE WE MIGHT NEED TO CHANGE...\cf0 \ |
|
205 |
gs.mapcalc(\cf5 '$\{lst\} = if((qc_this_day & 0x03)==0, $\{lst\}, null())'\cf0 ,\ |
|
206 |
lst=lstname, overwrite=True)\ |
|
207 |
\cf2 # clean up\cf0 \ |
|
208 |
gs.run_command(\cf5 'g.remove'\cf0 , rast=\cf5 'qc_this_day'\cf0 )\ |
|
209 |
\cf2 # return name of qc-adjusted LST raster in GRASS\cf0 \ |
|
210 |
\cf3 return\cf0 lstname\ |
|
211 |
\ |
|
212 |
\ |
|
213 |
\cf3 def\cf0 main():\ |
|
214 |
\cf2 #--------------------------------------------\cf0 \ |
|
215 |
\cf2 # Download and Calculate monthly climatology for daily LST time series for a specific area\cf0 \ |
|
216 |
\cf2 #--------------------------------------------\cf0 \ |
|
217 |
\ |
|
218 |
\cf2 # TODO: set up a (temporary?) GRASS database to use for processing? code\cf0 \ |
|
219 |
\cf2 # currently assumes it's being run within an existing GRASS session\cf0 \ |
|
220 |
\cf2 # using an appropriately configured LOCATION...\cf0 \ |
|
221 |
\cf2 #\cf0 \ |
|
222 |
\cf2 # note the following trick to fix datum for modis sinu;\cf0 \ |
|
223 |
\cf2 # TODO: check that this works properly...compare to MRT output?\cf0 \ |
|
224 |
\cf2 # gs.run_command('g.proj', flags='c',\cf0 \ |
|
225 |
\cf2 # proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere')\cf0 \ |
|
226 |
\cf2 ## proj4='+proj=sinu +R=6371007.181 +nadgrids=@null +wktext')\cf0 \ |
|
227 |
\cf5 \ |
|
228 |
\ |
|
229 |
\cf6 #### Added by Benoit on 18 October 2012\ |
|
230 |
\ |
|
231 |
#Parameters\cf5 \ |
|
232 |
\cf0 tiles = ['h11v08','h11v07','h12v08']\ |
|
233 |
start_year = 2001\ |
|
234 |
end_year = 2010\ |
|
235 |
end_month=12\ |
|
236 |
start_month=1\ |
|
237 |
\ |
|
238 |
\cf2 ################## First step: download data ###################\cf0 \ |
|
239 |
\cf2 # Added tile loop \cf0 \ |
|
240 |
\cf3 for\cf0 tile \cf3 in\cf0 tiles: \ |
|
241 |
\cf3 for\cf0 year \cf3 in\cf0 range(start_year, end_year+\cf5 1\cf0 ):\ |
|
242 |
if year==end_year:\ |
|
243 |
start_doy, end_doy = get_doy_range(year, end_month)\ |
|
244 |
start_doy=1 #Fix problem later\ |
|
245 |
hdfs = download_mod11a1(hdfdir, tile, start_doy, end_doy, year)\ |
|
246 |
\ |
|
247 |
\cf3 elif\cf0 year==start_year:\ |
|
248 |
start_doy, end_doy = get_doy_range(year, start_month)\ |
|
249 |
end_doy=365\ |
|
250 |
if calendar.isleap(year):\ |
|
251 |
end_doy=366\ |
|
252 |
hdfs = download_mod11a1(hdfdir, tile, start_doy, end_doy, year)\ |
|
253 |
\ |
|
254 |
\cf3 elif\cf0 (year!=start_year and year!=end_year):\ |
|
255 |
start_doy=1\ |
|
256 |
end_doy=365\ |
|
257 |
if calendar.isleap(year):\ |
|
258 |
end_doy=366\ |
|
259 |
hdfs = download_mod11a1(hdfdir, tile, start_doy, end_doy, year)\ |
|
260 |
\ |
|
261 |
\cf2 # modify loop to take into account "hdfs", list of hdf files needed in the next steps\'85 \cf0 \ |
|
262 |
\ |
|
263 |
\cf2 ################# Second step: compute climatology #################\ |
|
264 |
\cf0 \ |
|
265 |
\ |
|
266 |
\cf2 # Set up a temporary GRASS data base \ |
|
1 |
# -*- coding: utf-8 -*- |
|
2 |
""" |
|
3 |
Created on Fri Oct 26 23:48:57 2012 |
|
4 |
|
|
5 |
@author: benoitparmentier |
|
6 |
""" |
|
7 |
|
|
8 |
#!/usr/bin/python |
|
9 |
# |
|
10 |
# Early version of assorted Python code to download MODIS 11A1 (daily) |
|
11 |
# data associated with specified dates and tiles, then interface with |
|
12 |
# GRASS to calculate temporally averaged LST in a way that accounts for |
|
13 |
# QC flags. |
|
14 |
# |
|
15 |
# TODO: |
|
16 |
# - functionalize to encapsulate high level procedural steps |
|
17 |
# - create a 'main' function to allow script to be run as a command |
|
18 |
# that can accept arguments? |
|
19 |
# - put downloaded data somewhere other than working directory? |
|
20 |
# - write code to set up and take down a temporary GRASS location/mapset |
|
21 |
# - write code to export climatology rasters to file (geotiff? compressed?) |
|
22 |
# - calculate other aggregate stats? (stddev, ...) |
|
23 |
# - deal with nightly LST too? |
|
24 |
# - deal with more than just first two bits of QC flags? |
|
25 |
# e.g.: r.mapcalc 'qc_level2 = (qc>>2) & 0x03' |
|
26 |
# 0x03?? that means 0000_0011 and >>2 means shit twice on the right for |
|
27 |
# - record all downloads to a log file? |
|
28 |
# |
|
29 |
# Jim Regetz |
|
30 |
# NCEAS |
|
31 |
# Created on 16-May-2012 |
|
32 |
# Modified by Benoit Parmentier |
|
33 |
# NCEAS on 18-October-2012 |
|
34 |
# |
|
35 |
|
|
36 |
import os, glob |
|
37 |
import datetime, calendar |
|
38 |
import ftplib |
|
39 |
import grass.script as gs |
|
40 |
#from datetime import date |
|
41 |
#from datetime import datetime |
|
42 |
|
|
43 |
|
|
44 |
#------------------ |
|
45 |
# helper functions |
|
46 |
#------------------ |
|
47 |
|
|
48 |
def yj_to_ymd(year, doy): |
|
49 |
"""Return date as e.g. '2000.03.05' based on specified year and |
|
50 |
numeric day-of-year (doy) """ |
|
51 |
date = datetime.datetime.strptime('%d%03d' % (year, doy), '%Y%j') |
|
52 |
return date.strftime('%Y.%m.%d') |
|
53 |
|
|
54 |
def get_doy_range(year, month): |
|
55 |
"""Determine starting and ending numeric day-of-year (doy) |
|
56 |
asscociated with the specified month and year. |
|
57 |
|
|
58 |
Arguments: |
|
59 |
year -- four-digit integer year |
|
60 |
month -- integer month (1-12) |
|
61 |
|
|
62 |
Returns tuple of start and end doy for that month/year. |
|
63 |
""" |
|
64 |
last_day_of_month = calendar.monthrange(year, month)[1] |
|
65 |
start_doy = int(datetime.datetime.strptime('%d.%02d.%02d' % (year, |
|
66 |
month, 1), '%Y.%m.%d').strftime('%j')) |
|
67 |
end_doy = int(datetime.datetime.strptime('%d.%02d.%02d' % (year, |
|
68 |
month, last_day_of_month), '%Y.%m.%d').strftime('%j')) |
|
69 |
return (int(start_doy), int(end_doy)) |
|
70 |
|
|
71 |
#added by Benoit, to create a list of raster images for particular month (2001-2010) |
|
72 |
|
|
73 |
def date_sequence(year_start,month_start,day_start,year_end,month_end,day_end): |
|
74 |
"""Create a date sequence for a start date and end date defined by year, month, date |
|
75 |
Arguments: |
|
76 |
year_start -- starting year |
|
77 |
|
|
78 |
Returns list of absolute paths to the downloaded HDF files.Note the format yyyydoy |
|
79 |
""" |
|
80 |
from datetime import date |
|
81 |
from dateutil.rrule import rrule, DAILY |
|
82 |
|
|
83 |
a = date(year_start, month_start,day_start) |
|
84 |
b = date(year_end, month_end, day_end) |
|
85 |
date_seq=list() |
|
86 |
for dt in rrule(DAILY, dtstart=a, until=b): |
|
87 |
print dt.strftime("%Y%j") |
|
88 |
date_seq.append(dt.strftime("%Y%j")) |
|
89 |
|
|
90 |
return(date_seq) |
|
91 |
|
|
92 |
#Added on January 16, 2012 by Benoit NCEAS |
|
93 |
def list_raster_per_month(listmaps): |
|
94 |
"""Determine starting and ending numeric day-of-year (doy) |
|
95 |
asscociated with the specified month and year. |
|
96 |
|
|
97 |
Arguments: |
|
98 |
maplist -- list of raster grass-digit integer year |
|
99 |
month -- integer month (1-12) |
|
100 |
year_range -- list of years |
|
101 |
|
|
102 |
Returns tuple of start and end doy for that month/year. |
|
103 |
""" |
|
104 |
list_maps_month=list() |
|
105 |
nb_month=12 |
|
106 |
for i in range(1,nb_month+1): |
|
107 |
#i=i+1 |
|
108 |
#filename = listfiles_wd2[i] #list the files of interest |
|
109 |
#monthlist[:]=[] #empty the list |
|
110 |
monthlist=list() |
|
111 |
#monthlist2[:]=[] #empty the list |
|
112 |
j=0 #counter |
|
113 |
for mapname in listmaps: |
|
114 |
#tmp_str=LST[0] |
|
115 |
date_map=mapname.split('_')[1][1:8] |
|
116 |
#date = filename[filename.rfind('_MOD')-8:filename.rfind('_MOD')] |
|
117 |
d=datetime.datetime.strptime(date_map,'%Y%j') #This produces a date object from a string extracted from the file name. |
|
118 |
month=d.strftime('%m') #Extract the month of the current map |
|
119 |
if int(month)==i: |
|
120 |
#monthlist.append(listmaps[j]) |
|
121 |
monthlist.append(mapname) |
|
122 |
#monthlist2.append(listmaps[j]) |
|
123 |
j=j+1 |
|
124 |
list_maps_month.append(monthlist) #This is a list of a list containing a list for each month |
|
125 |
#list_files2.append(monthlist2) |
|
126 |
return(list_maps_month) |
|
267 | 127 |
|
268 |
\fs22 \cf0 \ |
|
269 |
|
|
270 |
\fs28 tmp_location = \cf5 'tmp'\cf0 + str(os.getpid()) \cf2 # Name of the temporary GRASS data base \cf0 \ |
|
271 |
orig_location = gs.gisenv()[\cf5 'LOCATION_NAME'\cf0 ]\ |
|
272 |
orig_mapset = gs.gisenv()[\cf5 'MAPSET'\cf0 ] |
|
273 |
\fs22 \ |
|
274 |
|
|
275 |
\fs28 gs.os.environ[\cf5 'GRASS_OVERWRITE'\cf0 ] = \cf5 '1' \cf2 # Allow for overwrite?? GRASS data base\cf5 \cf0 \ |
|
276 |
\cf2 \ |
|
277 |
\ |
|
278 |
#load a file that will define the region of interest and location at the same time\cf0 \ |
|
279 |
name_of_file = \cf5 'MOD11A1.A2003364.h12v08.005.2008163211649.hdf'\cf0 \ |
|
280 |
\cf2 # first load daytime LST, by accessing specific layers in the hdd??\ |
|
281 |
#SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MOD11A1.A2003364.h12v08.005.2008163211649.hdf":MODIS_Grid_Daily_1km_LST:LST_Day_1km\ |
|
282 |
\cf0 path_file = os.getcwd() \ |
|
283 |
lst1day = \cf5 'HDF4_EOS:EOS_GRID:%s/%s:%s'\cf0 % (\ |
|
284 |
path_file,\ |
|
285 |
name_of_file,\ |
|
286 |
\cf5 '\cf2 MODIS_Grid_Daily_1km_LST:LST_Day_1km\cf5 '\cf0 )\ |
|
287 |
gs.run_command(\cf5 'r.in.gdal'\cf0 , input=lst1day, output=\cf5 'LST_1day_h12v08'\cf0 ,\ |
|
288 |
location=tmp_location)\ |
|
289 |
\ |
|
290 |
\cf2 # Now that the new location has been create, switch to new location\cf0 \ |
|
291 |
gs.run_command(\cf5 'g.gisenv'\cf0 , set=\cf5 'LOCATION_NAME=%s'\cf0 % tmp_location)\ |
|
292 |
gs.run_command(\cf5 'g.gisenv'\cf0 , set=\cf5 'MAPSET=PERMANENT'\cf0 )\ |
|
293 |
\ |
|
294 |
\cf2 # set GRASS the projection system and GRASS region to match the extent of the file\cf0 \ |
|
295 |
gs.run_command(\cf5 'g.proj'\cf0 , flags=\cf5 'c'\cf0 ,\ |
|
296 |
proj4=\cf5 '+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere'\cf0 )\ |
|
297 |
gs.run_command(\cf5 'g.region'\cf0 , rast=\cf5 'LST_1day_h12v08'\cf0 )\ |
|
298 |
\cf2 \ |
|
299 |
\ |
|
300 |
#generate monthly pixelwise mean & count of high-quality daytime LST\ |
|
301 |
\cf0 \ |
|
302 |
gs.os.environ[\cf5 'GRASS_OVERWRITE'\cf0 ] = \cf5 '1'\ |
|
303 |
\ |
|
304 |
\cf2 #Provide a list of file "hdfs" to process\'85\ |
|
305 |
\cf0 tile= 'h12v08'\cf2 \ |
|
306 |
\cf0 fileglob = \cf5 'MOD11A1.A*.%s*hdf'\cf0 % (tile)\ |
|
307 |
pathglob = os.path.join(path_file, fileglob)\ |
|
308 |
files = glob.glob(pathglob)\ |
|
309 |
hdfs = files\ |
|
310 |
\cf2 ### LST values in images must be masked out if they have low quality: this is done using QC flags from the layer subset in hdf files\cf0 \ |
|
311 |
LST = [load_qc_adjusted_lst(hdf) \cf3 for\cf0 hdf \cf3 in\cf0 hdfs]\ |
|
312 |
\cf2 ### LST is a list that contains the name of the new lst files in the GRASS format. The files are stored in the temporary GRASS database.\cf0 \ |
|
313 |
clims = calc_clim(LST, \cf5 'LST_%s_%d_%02d'\cf0 % (tile, year, month))\ |
|
314 |
\cf2 # clean up\cf0 \ |
|
315 |
gs.run_command(\cf5 'g.remove'\cf0 , rast=\cf5 ','\cf0 .join(LST))\ |
|
316 |
gs.os.environ[\cf5 'GRASS_OVERWRITE'\cf0 ] = \cf5 '0'\cf0 \ |
|
317 |
\ |
|
318 |
\cf2 # clean up\cf0 \ |
|
319 |
gs.run_command(\cf5 'g.gisenv'\cf0 , set=\cf5 'LOCATION_NAME=%s'\cf0 % orig_location)\ |
|
320 |
gs.run_command(\cf5 'g.gisenv'\cf0 , set=\cf5 'MAPSET=%s'\cf0 % orig_mapset)\ |
|
321 |
shutil.rmtree(os.path.join(gs.gisenv()[\cf5 'GISDBASE'\cf0 ], tmp_location)) |
|
322 |
\fs22 \ |
|
323 |
|
|
324 |
\fs28 \ |
|
325 |
\cf3 return\cf0 \cf3 None\cf0 \ |
|
326 |
} |
|
128 |
# quick function to return list of dirs in wd |
|
129 |
def list_contents(ftp): |
|
130 |
"""Parse ftp directory listing into list of names of the files |
|
131 |
and/or directories contained therein. May or may not be robust in |
|
132 |
general, but seems to work fine for LP DAAP ftp server.""" |
|
133 |
listing = [] |
|
134 |
ftp.dir(listing.append) |
|
135 |
contents = [item.split()[-1] for item in listing[1:]] |
|
136 |
return contents |
|
137 |
|
|
138 |
def download_mod11a1(destdir, tile, start_doy, end_doy, year, ver=5): |
|
139 |
"""Download into destination directory the MODIS 11A1 HDF files |
|
140 |
matching a given tile, year, and day range. If for some unexpected |
|
141 |
reason there are multiple matches for a given day, only the first is |
|
142 |
used. If no matches, the day is skipped with a warning message. |
|
143 |
|
|
144 |
Arguments: |
|
145 |
destdir -- path to local destination directory |
|
146 |
tile -- tile identifier (e.g., 'h08v04') |
|
147 |
start_doy -- integer start of day range (0-366) |
|
148 |
end_doy -- integer end of day range (0-366) |
|
149 |
year -- integer year (>=2000) |
|
150 |
ver -- MODIS version (4 or 5) |
|
151 |
|
|
152 |
Returns list of absolute paths to the downloaded HDF files. |
|
153 |
""" |
|
154 |
# connect to data pool and change to the right directory |
|
155 |
ftp = ftplib.FTP('e4ftl01.cr.usgs.gov') |
|
156 |
ftp.login('anonymous', '') |
|
157 |
ftp.cwd('MOLT/MOD11A1.%03d' % ver) |
|
158 |
# make list of daily directories to traverse |
|
159 |
available_days = list_contents(ftp) |
|
160 |
desired_days = [yj_to_ymd(year, x) for x in range(start_doy, end_doy+1)] |
|
161 |
days_to_get = filter(lambda day: day in desired_days, available_days) |
|
162 |
if len(days_to_get) < len(desired_days): |
|
163 |
missing_days = [day for day in desired_days if day not in days_to_get] |
|
164 |
print 'skipping %d day(s) not found on server:' % len(missing_days) |
|
165 |
print '\n'.join(missing_days) |
|
166 |
# get each tile in turn |
|
167 |
hdfs = list() |
|
168 |
for day in days_to_get: |
|
169 |
ftp.cwd(day) |
|
170 |
files_to_get = [file for file in list_contents(ftp) |
|
171 |
if tile in file and file[-3:]=='hdf'] |
|
172 |
if len(files_to_get)==0: |
|
173 |
# no file found -- print message and move on |
|
174 |
print 'no hdf found on server for tile', tile, 'on', day |
|
175 |
ftp.cwd('..') #added by Benoit |
|
176 |
continue |
|
177 |
elif 1<len(files_to_get): |
|
178 |
# multiple files found! -- just use the first... |
|
179 |
print 'multiple hdfs found on server for tile', tile, 'on', day |
|
180 |
file = files_to_get[0] |
|
181 |
local_file = os.path.join(destdir, file) |
|
182 |
ftp.retrbinary('RETR %s' % file, open(local_file, 'wb').write) |
|
183 |
hdfs.append(os.path.abspath(local_file)) |
|
184 |
ftp.cwd('..') |
|
185 |
# politely disconnect |
|
186 |
ftp.quit() |
|
187 |
# return list of downloaded paths |
|
188 |
return hdfs |
|
189 |
|
|
190 |
def get_hdf_paths(hdfdir, tile, start_doy, end_doy, year): |
|
191 |
"""Look in supplied directory to find the MODIS 11A1 HDF files |
|
192 |
matching a given tile, year, and day range. If for some unexpected |
|
193 |
reason there are multiple matches for a given day, only the first is |
|
194 |
used. If no matches, the day is skipped with a warning message. |
|
195 |
|
|
196 |
Arguments: |
|
197 |
hdfdir -- path to directory containing the desired HDFs |
|
198 |
tile -- tile identifier (e.g., 'h08v04') |
|
199 |
start_doy -- integer start of day range (0-366) |
|
200 |
end_doy -- integer end of day range (0-366) |
|
201 |
year -- integer year (>=2000) |
|
202 |
|
|
203 |
Returns list of absolute paths to the located HDF files. |
|
204 |
""" |
|
205 |
hdfs = list() |
|
206 |
for doy in range(start_doy, end_doy+1): |
|
207 |
fileglob = 'MOD11A1.A%d%03d.%s*hdf' % (year, doy, tile) |
|
208 |
pathglob = os.path.join(hdfdir, fileglob) |
|
209 |
files = glob.glob(pathglob) |
|
210 |
if len(files)==0: |
|
211 |
# no file found -- print message and move on |
|
212 |
print 'cannot access %s: no such file' % pathglob |
|
213 |
continue |
|
214 |
elif 1<len(files): |
|
215 |
# multiple files found! -- just use the first... |
|
216 |
print 'multiple hdfs found for tile', tile, 'on', day |
|
217 |
hdfs.append(os.path.abspath(files[0])) |
|
218 |
return hdfs |
|
219 |
|
|
220 |
def calc_clim(maplist, name, overwrite=False): |
|
221 |
"""Generate some climatalogies in GRASS based on the input list of |
|
222 |
maps. As usual, current GRASS region settings apply. Produces the |
|
223 |
following output rasters: |
|
224 |
* nobs: count of number of (non-null) values over the input maps |
|
225 |
* mean: arithmetic mean of (non-null) values over the input maps |
|
226 |
|
|
227 |
Arguments: |
|
228 |
maplist -- list of names of GRASS maps to be aggregated |
|
229 |
name -- template (suffix) for GRASS output map names |
|
230 |
|
|
231 |
Returns list of names of the output maps created in GRASS. |
|
232 |
""" |
|
233 |
denominator = '(%s)' % '+'.join(['if(!isnull(%s))' % m |
|
234 |
for m in maplist]) |
|
235 |
gs.mapcalc('nobs_%s = %s' % (name, denominator), overwrite=overwrite) |
|
236 |
numerator = '(%s)' % '+'.join(['if(isnull(%s), 0, %s)' % (m, m) |
|
237 |
for m in maplist]) |
|
238 |
gs.mapcalc('mean_%s = round(float(%s)/nobs_%s)' % (name, numerator, name), |
|
239 |
overwrite=overwrite) |
|
240 |
return ['%s_%s' % (s, name) for s in ['nobs', 'mean']] |
|
241 |
|
|
242 |
def load_qc_adjusted_lst(hdf): |
|
243 |
"""Load LST_Day_1km into GRASS from the specified hdf file, and |
|
244 |
nullify any cells for which QA flags indicate anything other than |
|
245 |
high quality. |
|
246 |
|
|
247 |
Argument: |
|
248 |
hdf -- local path to the 11A1 HDF file |
|
249 |
|
|
250 |
Returns the name of the QC-adjusted LST map created in GRASS. |
|
251 |
""" |
|
252 |
lstname = 'LST_' + '_'.join(os.path.basename(hdf).split('.')[1:3]) |
|
253 |
# read in lst grid |
|
254 |
gs.run_command('r.in.gdal', |
|
255 |
input='HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:LST_Day_1km' % hdf, |
|
256 |
output=lstname) # No new location created since it has already been done with r.in.gdal before!!! |
|
257 |
|
|
258 |
# read in qc grid |
|
259 |
gs.run_command('r.in.gdal', |
|
260 |
input = 'HDF4_EOS:EOS_GRID:%s:MODIS_Grid_Daily_1km_LST:QC_Day' % hdf, |
|
261 |
output = 'qc_this_day') |
|
262 |
gs.run_command('g.region', rast=lstname) |
|
263 |
# null out any LST cells for which QC is not high [00] #THIS WHERE WE MIGHT NEED TO CHANGE... |
|
264 |
gs.mapcalc('${lst} = if((qc_this_day & 0x03)==0, ${lst}, null())', |
|
265 |
lst=lstname, overwrite=True) |
|
266 |
# clean up |
|
267 |
gs.run_command('g.remove', rast='qc_this_day') |
|
268 |
# return name of qc-adjusted LST raster in GRASS |
|
269 |
return lstname |
|
270 |
|
|
271 |
|
|
272 |
def main(): |
|
273 |
#-------------------------------------------- |
|
274 |
# Download and Calculate monthly climatology for daily LST time series for a specific area |
|
275 |
#-------------------------------------------- |
|
276 |
|
|
277 |
# TODO: set up a (temporary?) GRASS database to use for processing? code |
|
278 |
# currently assumes it's being run within an existing GRASS session |
|
279 |
# using an appropriately configured LOCATION... |
|
280 |
# |
|
281 |
# note the following trick to fix datum for modis sinu; |
|
282 |
# TODO: check that this works properly...compare to MRT output? |
|
283 |
# gs.run_command('g.proj', flags='c', |
|
284 |
# proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere') |
|
285 |
## proj4='+proj=sinu +R=6371007.181 +nadgrids=@null +wktext') |
|
286 |
|
|
287 |
|
|
288 |
#### Added by Benoit on 18 October 2012 |
|
289 |
|
|
290 |
#Parameters |
|
291 |
tiles = ['h11v08','h11v07','h12v08'] #These tiles correspond to Venezuela. |
|
292 |
start_year = 2001 |
|
293 |
end_year = 2010 |
|
294 |
end_month=12 |
|
295 |
start_month=1 |
|
296 |
hdfdir = '/home/parmentier/Data/benoit_test' #destination file where hdf files are stored locally after download. |
|
297 |
################## First step: download data ################### |
|
298 |
# Added tile loop |
|
299 |
year_list=range(start_year,end_year+1) #list of year to loop through |
|
300 |
year_list=[2001,2002] |
|
301 |
for tile in tiles: |
|
302 |
for year in year_list: |
|
303 |
if year==end_year: |
|
304 |
start_doy, end_doy = get_doy_range(year, end_month) |
|
305 |
start_doy=1 #Fix problem later |
|
306 |
hdfs = download_mod11a1(hdfdir, tile, start_doy, end_doy, year) |
|
307 |
|
|
308 |
elif year==start_year: |
|
309 |
start_doy, end_doy = get_doy_range(year, start_month) |
|
310 |
end_doy=365 |
|
311 |
if calendar.isleap(year): |
|
312 |
end_doy=366 |
|
313 |
hdfs = download_mod11a1(hdfdir, tile, start_doy, end_doy, year) |
|
314 |
|
|
315 |
elif (year!=start_year and year!=end_year): |
|
316 |
start_doy=1 |
|
317 |
end_doy=365 |
|
318 |
if calendar.isleap(year): |
|
319 |
end_doy=366 |
|
320 |
hdfs = download_mod11a1(hdfdir, tile, start_doy, end_doy, year) |
|
321 |
|
|
322 |
# modify loop to take into account "hdfs", list of hdf files needed in the next steps… |
|
323 |
|
|
324 |
################# Second step: compute climatology ################# |
|
325 |
|
|
326 |
# Set up a temporary GRASS data base |
|
327 |
|
|
328 |
tmp_location = 'tmp' + str(os.getpid()) # Name of the temporary GRASS data base |
|
329 |
orig_location = gs.gisenv()['LOCATION_NAME'] |
|
330 |
orig_mapset = gs.gisenv()['MAPSET'] |
|
331 |
gs.os.environ['GRASS_OVERWRITE'] = '1' # Allow for overwrite?? GRASS data base |
|
332 |
|
|
333 |
##NEED TO LOOP OVER DIFFERENT TILES!!! |
|
334 |
##start loop for tile |
|
335 |
#load a file that will define the region of interest and location at the same time |
|
336 |
|
|
337 |
tile= 'h12v08' |
|
338 |
tile=tiles[2] |
|
339 |
#tile= 'h12v08' |
|
340 |
|
|
341 |
path_file = '/home/parmentier/Data/benoit_test' |
|
342 |
os.chdir(path_file) #set working directory |
|
343 |
os.getcwd() #get current working directory |
|
344 |
listfiles_wd2 = glob.glob('*'+tile+'*'+'.hdf') #list the all the LST files of interest |
|
345 |
name_of_file = listfiles_wd2[0] |
|
346 |
#name_of_file = 'MOD11A1.A2003364.h12v08.005.2008163211649.hdf' |
|
347 |
# first load daytime LST, by accessing specific layers in the hdd?? |
|
348 |
#SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MOD11A1.A2003364.h12v08.005.2008163211649.hdf":MODIS_Grid_Daily_1km_LST:LST_Day_1km |
|
349 |
#path_file = os.getcwd() |
|
350 |
#Create the name of the layer to extract from the hdf file used for definition of the database |
|
351 |
lst1day = 'HDF4_EOS:EOS_GRID:%s/%s:%s' % ( |
|
352 |
path_file, |
|
353 |
name_of_file, |
|
354 |
'MODIS_Grid_Daily_1km_LST:LST_Day_1km') |
|
355 |
|
|
356 |
#now import one image and set the location, mapset and database !!create name per tile |
|
357 |
gs.run_command('r.in.gdal', input=lst1day, output='LST_1day_h12v08', |
|
358 |
location=tmp_location) |
|
359 |
|
|
360 |
# Now that the new location has been create, switch to new location |
|
361 |
gs.run_command('g.gisenv', set='LOCATION_NAME=%s' % tmp_location) |
|
362 |
gs.run_command('g.gisenv', set='MAPSET=PERMANENT') |
|
363 |
|
|
364 |
# set GRASS the projection system and GRASS region to match the extent of the file |
|
365 |
gs.run_command('g.proj', flags='c', |
|
366 |
proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere') |
|
367 |
gs.run_command('g.region', rast='LST_1day_h12v08') |
|
368 |
|
|
369 |
|
|
370 |
#generate monthly pixelwise mean & count of high-quality daytime LST |
|
371 |
|
|
372 |
gs.os.environ['GRASS_OVERWRITE'] = '1' |
|
373 |
|
|
374 |
#Provide a list of file "hdfs" to process… |
|
375 |
#tile= 'h12v08' |
|
376 |
fileglob = 'MOD11A1.A*.%s*hdf' % (tile) |
|
377 |
pathglob = os.path.join(path_file, fileglob) |
|
378 |
files = glob.glob(pathglob) |
|
379 |
hdfs = files |
|
380 |
### LST values in images must be masked out if they have low quality: this is done using QC flags from the layer subset in hdf files |
|
381 |
LST = [load_qc_adjusted_lst(hdf) for hdf in hdfs] |
|
382 |
### LST is a list that contains the name of the new lst files in the GRASS format. The files are stored in the temporary GRASS database. |
|
383 |
list_maps_month=list_raster_per_month(LST) |
|
384 |
list_maps_name=['jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec'] |
|
385 |
##Now calculate clim per month, do a test for year1 |
|
386 |
nb_month=12 |
|
387 |
for i in range(1,nb_month+1): |
|
388 |
#clims = calc_clim(LST, 'LST_%s_%d_%02d' % (tile, year, month)) |
|
389 |
clims = calc_clim(list_maps_month[i-1],list_maps_name[i-1]+'_'+str(i-1)) |
|
390 |
|
|
391 |
#clims = calc_clim(LST, 'LST_%s_%d_%02d' % (tile, year, month)) |
|
392 |
|
|
393 |
# clean up if necessary |
|
394 |
gs.run_command('g.remove', rast=','.join(LST)) |
|
395 |
gs.os.environ['GRASS_OVERWRITE'] = '0' |
|
396 |
|
|
397 |
# clean up |
|
398 |
gs.run_command('g.gisenv', set='LOCATION_NAME=%s' % orig_location) |
|
399 |
gs.run_command('g.gisenv', set='MAPSET=%s' % orig_mapset) |
|
400 |
shutil.rmtree(os.path.join(gs.gisenv()['GISDBASE'], tmp_location)) |
|
401 |
|
|
402 |
return None |
Also available in: Unified diff
LST climatology, added function to create list of files per month and other modifications