1
|
# -*- coding: utf-8 -*-
|
2
|
"""
|
3
|
Created on Mon May 13 10:02:02 2013
|
4
|
|
5
|
@author: parmentier
|
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
|
# - record all downloads to a log file?
|
21
|
#
|
22
|
# Created by Jim Regetz
|
23
|
# NCEAS on 16-May-2012
|
24
|
# Modified by Benoit Parmentier
|
25
|
#NCEAS on 18-October-2012
|
26
|
#NCEAS on 19-January-2013
|
27
|
#NCEAS on 14-May-2013
|
28
|
#
|
29
|
|
30
|
import os, glob
|
31
|
import datetime, calendar
|
32
|
import ftplib
|
33
|
#import grass.script as gs
|
34
|
import argparse
|
35
|
import errno
|
36
|
#from datetime import date
|
37
|
#from datetime import datetime
|
38
|
|
39
|
|
40
|
#------------------
|
41
|
# helper functions
|
42
|
#------------------
|
43
|
|
44
|
def yj_to_ymd(year, doy):
|
45
|
"""Return date as e.g. '2000.03.05' based on specified year and
|
46
|
numeric day-of-year (doy) """
|
47
|
date = datetime.datetime.strptime('%d%03d' % (year, doy), '%Y%j')
|
48
|
return date.strftime('%Y.%m.%d')
|
49
|
|
50
|
def get_doy_range(year, month):
|
51
|
"""Determine starting and ending numeric day-of-year (doy)
|
52
|
asscociated with the specified month and year.
|
53
|
|
54
|
Arguments:
|
55
|
year -- four-digit integer year
|
56
|
month -- integer month (1-12)
|
57
|
|
58
|
Returns tuple of start and end doy for that month/year.
|
59
|
"""
|
60
|
last_day_of_month = calendar.monthrange(year, month)[1]
|
61
|
start_doy = int(datetime.datetime.strptime('%d.%02d.%02d' % (year,
|
62
|
month, 1), '%Y.%m.%d').strftime('%j'))
|
63
|
end_doy = int(datetime.datetime.strptime('%d.%02d.%02d' % (year,
|
64
|
month, last_day_of_month), '%Y.%m.%d').strftime('%j'))
|
65
|
return (int(start_doy), int(end_doy))
|
66
|
|
67
|
#added by Benoit, to create a list of raster images for particular month (2001-2010)
|
68
|
|
69
|
def date_sequence(year_start,month_start,day_start,year_end,month_end,day_end):
|
70
|
"""Create a date sequence for a start date and end date defined by year, month, date
|
71
|
Arguments:
|
72
|
year_start -- starting year
|
73
|
|
74
|
Returns list of absolute paths to the downloaded HDF files.Note the format yyyydoy
|
75
|
"""
|
76
|
from datetime import date
|
77
|
from dateutil.rrule import rrule, DAILY
|
78
|
|
79
|
a = date(year_start, month_start,day_start)
|
80
|
b = date(year_end, month_end, day_end)
|
81
|
date_seq=list()
|
82
|
for dt in rrule(DAILY, dtstart=a, until=b):
|
83
|
print dt.strftime("%Y%j")
|
84
|
date_seq.append(dt.strftime("%Y%j"))
|
85
|
|
86
|
return(date_seq)
|
87
|
|
88
|
# quick function to return list of dirs in wd
|
89
|
def list_contents(ftp):
|
90
|
"""Parse ftp directory listing into list of names of the files
|
91
|
and/or directories contained therein. May or may not be robust in
|
92
|
general, but seems to work fine for LP DAAP ftp server."""
|
93
|
listing = []
|
94
|
ftp.dir(listing.append)
|
95
|
contents = [item.split()[-1] for item in listing[1:]]
|
96
|
return contents
|
97
|
|
98
|
def download_mod11a1(destdir, tile, start_doy, end_doy, year, ver=5):
|
99
|
"""Download into destination directory the MODIS 11A1 HDF files
|
100
|
matching a given tile, year, and day range. If for some unexpected
|
101
|
reason there are multiple matches for a given day, only the first is
|
102
|
used. If no matches, the day is skipped with a warning message.
|
103
|
|
104
|
Arguments:
|
105
|
destdir -- path to local destination directory
|
106
|
tile -- tile identifier (e.g., 'h08v04')
|
107
|
start_doy -- integer start of day range (0-366)
|
108
|
end_doy -- integer end of day range (0-366)
|
109
|
year -- integer year (>=2000)
|
110
|
ver -- MODIS version (4 or 5)
|
111
|
|
112
|
Returns list of absolute paths to the downloaded HDF files.
|
113
|
"""
|
114
|
# connect to data pool and change to the right directory
|
115
|
ftp = ftplib.FTP('e4ftl01.cr.usgs.gov')
|
116
|
ftp.login('anonymous', '')
|
117
|
ftp.cwd('MOLT/MOD11A1.%03d' % ver)
|
118
|
# make list of daily directories to traverse
|
119
|
available_days = list_contents(ftp)
|
120
|
desired_days = [yj_to_ymd(year, x) for x in range(start_doy, end_doy+1)]
|
121
|
days_to_get = filter(lambda day: day in desired_days, available_days)
|
122
|
if len(days_to_get) < len(desired_days):
|
123
|
missing_days = [day for day in desired_days if day not in days_to_get]
|
124
|
print 'skipping %d day(s) not found on server:' % len(missing_days)
|
125
|
print '\n'.join(missing_days)
|
126
|
# get each tile in turn
|
127
|
hdfs = list()
|
128
|
for day in days_to_get:
|
129
|
ftp.cwd(day)
|
130
|
files_to_get = [file for file in list_contents(ftp)
|
131
|
if tile in file and file[-3:]=='hdf']
|
132
|
if len(files_to_get)==0:
|
133
|
# no file found -- print message and move on
|
134
|
print 'no hdf found on server for tile', tile, 'on', day
|
135
|
ftp.cwd('..') #added by Benoit
|
136
|
continue
|
137
|
elif 1<len(files_to_get):
|
138
|
# multiple files found! -- just use the first...
|
139
|
print 'multiple hdfs found on server for tile', tile, 'on', day
|
140
|
file = files_to_get[0]
|
141
|
local_file = os.path.join(destdir, file)
|
142
|
ftp.retrbinary('RETR %s' % file, open(local_file, 'wb').write)
|
143
|
hdfs.append(os.path.abspath(local_file))
|
144
|
ftp.cwd('..')
|
145
|
# politely disconnect
|
146
|
ftp.quit()
|
147
|
# return list of downloaded paths
|
148
|
return hdfs
|
149
|
|
150
|
def get_hdf_paths(hdfdir, tile, start_doy, end_doy, year):
|
151
|
"""Look in supplied directory to find the MODIS 11A1 HDF files
|
152
|
matching a given tile, year, and day range. If for some unexpected
|
153
|
reason there are multiple matches for a given day, only the first is
|
154
|
used. If no matches, the day is skipped with a warning message.
|
155
|
|
156
|
Arguments:
|
157
|
hdfdir -- path to directory containing the desired HDFs
|
158
|
tile -- tile identifier (e.g., 'h08v04')
|
159
|
start_doy -- integer start of day range (0-366)
|
160
|
end_doy -- integer end of day range (0-366)
|
161
|
year -- integer year (>=2000)
|
162
|
|
163
|
Returns list of absolute paths to the located HDF files.
|
164
|
"""
|
165
|
hdfs = list()
|
166
|
for doy in range(start_doy, end_doy+1):
|
167
|
fileglob = 'MOD11A1.A%d%03d.%s*hdf' % (year, doy, tile)
|
168
|
pathglob = os.path.join(hdfdir, fileglob)
|
169
|
files = glob.glob(pathglob)
|
170
|
if len(files)==0:
|
171
|
# no file found -- print message and move on
|
172
|
print 'cannot access %s: no such file' % pathglob
|
173
|
continue
|
174
|
elif 1<len(files):
|
175
|
# multiple files found! -- just use the first...
|
176
|
print 'multiple hdfs found for tile', tile, 'on', day
|
177
|
hdfs.append(os.path.abspath(files[0]))
|
178
|
return hdfs
|
179
|
|
180
|
def create_dir_and_check_existence(path):
|
181
|
try:
|
182
|
os.makedirs(path)
|
183
|
except OSError as exception:
|
184
|
if exception.errno != errno.EEXIST:
|
185
|
raise
|
186
|
|
187
|
def main():
|
188
|
#from sys import argv # example client code
|
189
|
#myargs = getopts(argv)
|
190
|
|
191
|
#import argparse
|
192
|
parser = argparse.ArgumentParser()
|
193
|
parser.add_argument("tiles", type=str, help="list of Modis tiles")
|
194
|
parser.add_argument("start_year", type=int, help="start year")
|
195
|
parser.add_argument("end_year", type=int, help="end year")
|
196
|
parser.add_argument("start_month", type=int, help="start month")
|
197
|
parser.add_argument("end_month", type=int, help="end month")
|
198
|
parser.add_argument("hdfdir", type=str, help="destination/source directory for hdf file")
|
199
|
parser.add_argument("night", type=int, help="night")
|
200
|
parser.add_argument("download", type=int, help="out_suffix")
|
201
|
parser.add_argument("out_suffix", type=str, help="out_suffix")
|
202
|
|
203
|
myargs = parser.parse_args()
|
204
|
|
205
|
tiles = myargs.tiles #These tiles correspond to Venezuela.
|
206
|
start_year = myargs.start_year
|
207
|
end_year = myargs.end_year
|
208
|
end_month = myargs.end_month
|
209
|
start_month= myargs.start_month
|
210
|
hdfdir = myargs.hdfdir
|
211
|
night= myargs.night # if 1 then produce night climatology
|
212
|
out_suffix= myargs.out_suffix #"_03192013"
|
213
|
download=myargs.download# if 1 then download files
|
214
|
|
215
|
################## First step: download data ###################
|
216
|
# Added tile loop
|
217
|
year_list=range(start_year,end_year+1) #list of year to loop through
|
218
|
#print 'Number of arguments:', len(myargs), 'arguments.'
|
219
|
#print 'Argument List:', str(myargs)
|
220
|
|
221
|
tiles =tiles.split(",") #Create a list from string
|
222
|
print 'this is the result',str(year_list)
|
223
|
print(myargs)
|
224
|
print(tiles)
|
225
|
|
226
|
|
227
|
#tiles = ['h11v08','h11v07','h12v07','h12v08','h10v07','h10v08'] #These tiles correspond to Venezuela.
|
228
|
#start_year = 2001
|
229
|
#end_year = 2010
|
230
|
#end_month=12
|
231
|
#start_month=1
|
232
|
#hdfdir = '/home/parmentier/Data/benoit_test' #destination file where hdf files are stored locally after download.
|
233
|
#night=1 # if 1 then produce night climatology
|
234
|
#out_suffix="_03192013"
|
235
|
#download=1 # if 1 then download files
|
236
|
|
237
|
################## First step: download data ###################
|
238
|
# Added tile loop
|
239
|
year_list=range(start_year,end_year+1) #list of year to loop through
|
240
|
#for testing...
|
241
|
#year_list=[2001,2002]
|
242
|
#hdfdir = '/home/parmentier/Data/benoit_test2'
|
243
|
#tiles = ['h10v07','h10v08','h12v07']
|
244
|
if download==1:
|
245
|
for tile in tiles:
|
246
|
outDir = os.path.join(hdfdir,tile)
|
247
|
create_dir_and_check_existence(outDir)
|
248
|
|
249
|
for year in year_list:
|
250
|
start_doy = 1
|
251
|
end_doy = 365
|
252
|
if calendar.isleap(year):
|
253
|
end_doy=366
|
254
|
|
255
|
hdfs = download_mod11a1(outDir, tile, start_doy, end_doy, year)
|
256
|
|
257
|
if __name__ == '__main__':
|
258
|
main()
|
259
|
# from sys import argv # example client code
|
260
|
# myargs = getopts(argv)
|
261
|
# if '-a' in myargs:
|
262
|
# print(myargs['-a'])
|
263
|
# if '-b' in myargs:
|
264
|
# print(myargs['-b'])
|
265
|
# print(myargs)
|
266
|
|
267
|
#i = 1
|
268
|
#while i < arg_length:
|
269
|
# arg = argv[i]
|
270
|
# if arg =='-a':
|
271
|
# a = argv[i]
|
272
|
# elif arg =='-b':
|
273
|
# b = argv[i]
|
274
|
|
275
|
#x = a + b
|
276
|
|
277
|
#print x
|