Revision b8506b64
Added by Benoit Parmentier over 11 years ago
| climate/research/oregon/modis-lst/LST_downloading_and_climatology.py | ||
|---|---|---|
| 37 | 37 |
import datetime, calendar |
| 38 | 38 |
import ftplib |
| 39 | 39 |
import grass.script as gs |
| 40 |
|
|
| 40 | 41 |
#from datetime import date |
| 41 | 42 |
#from datetime import datetime |
| 42 | 43 |
|
| ... | ... | |
| 285 | 286 |
## proj4='+proj=sinu +R=6371007.181 +nadgrids=@null +wktext') |
| 286 | 287 |
|
| 287 | 288 |
|
| 288 |
#### MOdified by Benoit in October 2012 and January 2013
|
|
| 289 |
#### MOdified by Benoit in March 2013
|
|
| 289 | 290 |
|
| 290 | 291 |
#Parameters |
| 292 |
#Inputs from R?? |
|
| 291 | 293 |
tiles = ['h11v08','h11v07','h12v07','h12v08','h10v07','h10v08'] #These tiles correspond to Venezuela. |
| 292 | 294 |
start_year = 2001 |
| 293 | 295 |
end_year = 2010 |
| 294 | 296 |
end_month=12 |
| 295 | 297 |
start_month=1 |
| 296 | 298 |
hdfdir = '/home/parmentier/Data/benoit_test' #destination file where hdf files are stored locally after download. |
| 297 |
night=0 # if 1 then produce night climatology
|
|
| 298 |
out_suffix="_01252013"
|
|
| 299 |
night=1 # if 1 then produce night climatology
|
|
| 300 |
out_suffix="_03192013"
|
|
| 299 | 301 |
download=1 # if 1 then download files |
| 300 | 302 |
|
| 301 | 303 |
################## First step: download data ################### |
| ... | ... | |
| 340 | 342 |
os.getcwd() #get current working directory |
| 341 | 343 |
listfiles_wd2 = glob.glob('*'+tile+'*'+'.hdf') #list the all the LST files of interest
|
| 342 | 344 |
name_of_file = listfiles_wd2[0] |
| 343 |
#Create the name of the layer to extract from the hdf file used for definition of the database |
|
| 345 |
#Create the name of the layer to extract from the hdf file used for definition of the database
|
|
| 344 | 346 |
lst1day = 'HDF4_EOS:EOS_GRID:%s/%s:%s' % ( |
| 345 | 347 |
path_file, |
| 346 | 348 |
name_of_file, |
| 347 |
'MODIS_Grid_Daily_1km_LST:'+ lst_var)
|
|
| 348 |
#'LST_Night_1km'
|
|
| 349 |
#change the above line later to allow night LST for tmin
|
|
| 349 |
'MODIS_Grid_Daily_1km_LST:'+ lst_var) |
|
| 350 |
#'LST_Night_1km' |
|
| 351 |
#change the above line later to allow night LST for tmin |
|
| 350 | 352 |
|
| 351 | 353 |
#now import one image and set the location, mapset and database !!create name per tile |
| 352 | 354 |
gs.run_command('r.in.gdal', input=lst1day, output='LST_1day_'+tile,
|
| ... | ... | |
| 358 | 360 |
|
| 359 | 361 |
# set GRASS the projection system and GRASS region to match the extent of the file |
| 360 | 362 |
gs.run_command('g.proj', flags='c',
|
| 361 |
proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere') |
|
| 363 |
proj4='+proj=sinu +a=6371007.181 +b=6371007.181 +ellps=sphere') #This should be set in the parameters
|
|
| 362 | 364 |
gs.run_command('g.region', rast='LST_1day_'+tile)
|
| 363 | 365 |
|
| 364 | 366 |
#generate monthly pixelwise mean & count of high-quality daytime LST |
| ... | ... | |
| 378 | 380 |
##Now calculate clim per month, do a test for year1 |
| 379 | 381 |
nb_month = 12 |
| 380 | 382 |
for i in range(1,nb_month+1): |
| 381 |
clims = calc_clim(list_maps_month[i-1],tile+'_'+list_maps_name[i-1]+'_'+str(i-1)) |
|
| 383 |
clims = calc_clim(list_maps_month[i-1],lst_var+'_'+tile+'_'+list_maps_name[i-1]+'_'+str(i-1))
|
|
| 382 | 384 |
for j in range(1, len(clims)+1): |
| 383 | 385 |
if j-1 ==0: |
| 384 | 386 |
gs.run_command('r.out.gdal', input= clims[j-1], output=clims[j-1]+ out_suffix +'.tif', type='Float32')
|
Also available in: Unified diff
Climatology LST, modifications to allow LST night calculation for TMIN predictions