Project

General

Profile

« Previous | Next » 

Revision 195bbbb0

Added by Benoit Parmentier over 11 years ago

LST climatology, changes to output location and adding compression to tif output

View differences:

climate/research/oregon/modis-lst/LST_climatology.py
129 129 
    # return name of qc-adjusted LST raster in GRASS
130 130 
    return lstname
131 131 

  		




  
  132
  
    
def create_dir_and_check_existence(path):


  		




  
  133
  
    
    try:


  		




  
  134
  
    
        os.makedirs(path)


  		




  
  135
  
    
    except OSError as exception:


  		




  
  136
  
    
        if exception.errno != errno.EEXIST:


  		




  
  137
  
    
            raise


  		




  
  138
  
    
            


  		




  132
  139
  
    
def main():


  		




  133
  140
  
    
    #--------------------------------------------


  		




  134
  141
  
    
    # Download and Calculate monthly climatology for daily LST time series for a specific area


  		




  ......




  150
  157
  
    
    ### INPUT Parameters


  		




  151
  158
  
    
    #Inputs from R?? there are 9 parameters


  		




  152
  159
  
    
   #tiles = ['h11v08','h11v07','h12v07','h12v08','h10v07','h10v08'] #These tiles correspond to Venezuela.


  		




  153
  
  
    
    #tiles= ['h08v04','h09v04']    


  		




  154
  
  
    
    #start_year = 2001


  		




  155
  
  
    
    #end_year = 2010


  		




  156
  
  
    
    #start_month=1


  		




  157
  
  
    
    #end_month=12


  		




  
  160
  
    
    tiles= ['h08v04','h09v04']    


  		




  
  161
  
    
    start_year = 2001


  		




  
  162
  
    
    end_year = 2010


  		




  
  163
  
    
    start_month=1


  		




  
  164
  
    
    end_month=12


  		




  158
  165
  
    
    #hdfdir =  '/home/layers/commons/modis/MOD11A1_tiles' #destination file where hdf files are stored locally after download.


  		




  159
  
  
    
    #hdfdir =  '/home/parmentier/Data/IPLANT_project/MOD11A1_tiles' #destination file where hdf files are stored locally after download.


  		




  160
  
  
    
    #night=1    # if 1 then produce night climatology


  		




  161
  
  
    
    #out_suffix="_03192013"


  		




  162
  
  
    
    #download=1  # if 1 then download files


  		




  
  166
  
    
    hdfdir =  '/home/parmentier/Data/IPLANT_project/MOD11A1_tiles' #destination file where hdf files are stored locally after download.


  		




  
  167
  
    
    night=1    # if 1 then produce night climatology


  		




  
  168
  
    
    out_suffix="_03192013"


  		




  
  169
  
    
    download=0  # if 1 then download files


  		




  163
  170
  
    
   


  		




  164
  171
  
    
    #Passing arguments from the shell...using positional assignment


  		




  165
  172
  
    
    parser = argparse.ArgumentParser()


  		




  ......




  201
  208
  
    
    if night==0:


  		




  202
  209
  
    
        lst_var = var_name[0]


  		




  203
  210
  
    
        


  		




  
  211
  
    
    DirLST="LST_averages"


  		




  
  212
  
    
    outDir = os.path.join(hdfdir,DirLST) #create output filesfor LST averages


  		




  
  213
  
    
    create_dir_and_check_existence(outDir)


  		




  204
  214
  
    
    #start = time.clock()    


  		




  205
  215
  
    
    for tile in tiles:        


  		




  206
  216
  
    
        # Set up a temporary GRASS data base   


  		




  ......




  255
  265
  
    
        ##Now calculate clim per month, do a test for year1


  		




  256
  266
  
    
        nb_month = 12


  		




  257
  267
  
    
        for i in range(1,nb_month+1):


  		




  258
  
  
    
            clims = calc_clim(list_maps_month[i-1],lst_var+'_'+tile+'_'+list_maps_name[i-1]+'_'+str(i-1))


  		




  
  268
  
    
            clims = calc_clim(list_maps_month[i-1],lst_var+'_'+tile+'_'+list_maps_name[i-1]+'_'+str(i-1)) #length 2, list contains avg and nobs


  		




  259
  269
  
    
            for j in range(1, len(clims)+1):


  		




  260
  
  
    
                if j-1 ==0:


  		




  261
  
  
    
                    gs.run_command('r.out.gdal', input= clims[j-1], output=clims[j-1]+ out_suffix +'.tif', type='Float64')


  		




  262
  
  
    
                if j-1 ==1:


  		




  
  270
  
    
                if j-1 ==0:  #if image is nobs (number of observation, then output as such)


  		




  
  271
  
    
                    #gs.run_command('r.out.gdal', input= clims[j-1], output=clims[j-1]+ out_suffix +'.tif', type='Float64')


  		




  
  272
  
    
                    gs.run_command('r.out.gdal', input= clims[j-1], output=os.path.join(outDir,clims[j-1])+out_suffix+'.tif', type='Float64',createopt='COMPRESS=LZW')


  		




  
  273
  
    
                if j-1 ==1: #if image is avg clim then rescaled in celsius degrees


  		




  263
  274
  
    
                    gs.mapcalc(' clim_rescaled = ('+ clims[j-1 ]+ ' * 0.02) -273.15')  


  		




  264
  
  
    
                    gs.run_command('r.out.gdal', input= 'clim_rescaled', output=clims[j-1]+ out_suffix+'.tif', type='Float64')


  		




  
  275
  
    
                    #gs.run_command('r.out.gdal', input= 'clim_rescaled', output=clims[j-1]+ out_suffix+'.tif', type='Float64')


  		




  
  276
  
    
                    gs.run_command('r.out.gdal', input= 'clim_rescaled', output=os.path.join(outDir,clims[j-1])+out_suffix+'.tif', type='Float64',createopt='COMPRESS=LZW')


  		




  265
  277
  
    
        #clims = calc_clim(LST, 'LST_%s_%d_%02d' % (tile, year, month))


  		




  266
  278
  
    
        


  		




  267
  279
  
    
        # clean up GRASS DATABASE: ok working


  		












Also available in:  Unified diff