make_mat.py

   1 import saveload as sl
   2 import numpy as np
   3 import matplotlib.pyplot as plt
   4 from scipy.interpolate import griddata
   5 from scipy import spatial
   6 import scipy.io as sio
   7 from time import time
   8 import sys
   9 import os
  10
  11 t_init = time()
  12
  13 # use a test granule or all the granules?
  14 test = False
  15 # interpolate the whole granule?
  16 whole = False
  17 # radius options: 'max', 'mean', or 'min'
  18 opt_rad = 'max'
  19
  20 print '>> Options <<'
  21 if test:
  22         print '1> Test granule'
  23 else:
  24         print '1> All the granules'
  25 if whole:
  26         print '2> Interpolating the whole granule'
  27 else:
  28         print '2> Interpolating just the fire mesh bounding box'
  29 print '3> Radius computed using',opt_rad
  30 print ''
  31
  32 print '>> Loading data <<'
  33 # loading the data file
  34 data,fxlon,fxlat,time_num=sl.load('data')
  35 t1 = time()
  36 print 'elapsed time: ',abs(t_init-t1),'s'
  37 print ''
  38
  39 # mapping prefixes from JPSSD to cycling
  40 prefixes={'MOD': 'MOD14', 'MYD': 'MYD14', 'VNP': 'NPP_VAF_L2'}
  41
  42 # constants for geotransform
  43 if whole:
  44         res = 0.05   # resolution
  45 else:
  46         res = 0.01   # resolution
  47         fbounds = [fxlon.min(),fxlon.max(),fxlat.min(),fxlat.max()] # fire mesh bounds
  48 rot = 0.0    # rotation (not currently supported)
  49 print '>> General geotransform constants <<'
  50 print 'resolution =',res
  51 print 'rotation =',rot
  52 print ''
  53
  54 # defining granules to process
  55 if test:
  56         grans = ['VNP_A2018312_2130']
  57 else:
  58         grans = list(data)
  59
  60 NG = len(grans)
  61 print '>> Processing',NG,'granules <<'
  62 # for each granule
  63 for ng,gran in enumerate(grans):
  64         print '>> Processing granule',gran,'-','%d/%d' % (ng+1,NG)
  65         tg1 = time()
  66         # creating file_name
  67         splitted=gran.split('_')
  68         prefix=splitted[0]
  69         date=splitted[1][3:8]+splitted[2]+'00'
  70         file_name = prefixes[prefix]+'.'+date+'.tif.mat'
  71
  72         if not os.path.isfile(file_name):
  73                 # creating meshgrid where to interpolate the data
  74                 bounds = [data[gran].lon.min(),data[gran].lon.max(),data[gran].lat.min(),data[gran].lat.max()]
  75                 lons_interp = np.arange(bounds[0],bounds[1],res)
  76                 lats_interp = np.arange(bounds[2],bounds[3],res)
  77                 lons_interp,lats_interp = np.meshgrid(lons_interp,lats_interp)
  78                 glons = np.reshape(lons_interp,np.prod(lons_interp.shape))
  79                 glats = np.reshape(lats_interp,np.prod(lats_interp.shape))
  80
  81                 # initialize fire_interp
  82                 if not whole:
  83                         fires_interp = np.zeros(np.prod(lons_interp.shape))
  84
  85                 # creating geotransform array with necessary geolocation information
  86                 geotransform = [bounds[0],res,rot,bounds[3],rot,res]
  87
  88                 # approximation of the radius for the tree
  89                 dlon = abs(data[gran].lon[1,1]-data[gran].lon[0,0])
  90                 dlat = abs(data[gran].lat[1,1]-data[gran].lat[0,0])
  91                 if opt_rad == 'max':
  92                         radius = max(dlon,dlat)
  93                 elif opt_rad == 'mean':
  94                         radius = (dlon+dlat)*.5
  95                 elif opt_rad == 'min':
  96                         radius = min(dlon,dlat)
  97                 print 'radius =',radius
  98
  99                 # flatten granule data into 1d arrays
 100                 lons = np.reshape(data[gran].lon,np.prod(data[gran].lon.shape))
 101                 lats = np.reshape(data[gran].lat,np.prod(data[gran].lat.shape))
 102                 fires = np.reshape(data[gran].fire,np.prod(data[gran].fire.shape)).astype(np.int8)
 103                 if not whole:
 104                         M = (lons>fbounds[0])*(lons<fbounds[1])*(lats>fbounds[2])*(lats<fbounds[3])
 105                         Mg = (glons>fbounds[0])*(glons<fbounds[1])*(glats>fbounds[2])*(glats<fbounds[3])
 106                         lons = lons[M]
 107                         lats = lats[M]
 108                         fires = fires[M]
 109                         glons = glons[Mg]
 110                         glats = glats[Mg]
 111
 112                 if len(lons) > 0:
 113                         # Making tree
 114                         print '> Making the tree...'
 115                         t1 = time()
 116                         tree = spatial.cKDTree(np.column_stack((lons,lats)))
 117                         t2 = time()
 118                         print 'elapsed time: ',abs(t2-t1),'s'
 119                         print '> Interpolating the data...'
 120                         t1 = time()
 121                         indexes = np.array(tree.query_ball_point(np.column_stack((glons,glats)),radius,return_sorted=True))
 122                         t2 = time()
 123                         print 'elapsed time: ',abs(t2-t1),'s'
 124                         filtered_indexes = np.array([index[0] if len(index) > 0 else np.nan for index in indexes])
 125                         fire1d = [fires[int(ii)] if not np.isnan(ii) else 0 for ii in filtered_indexes]
 126                         if whole:
 127                                 fires_interp = np.reshape(fire1d,lons_interp.shape)
 128                         else:
 129                                 fires_interp[Mg] = fire1d
 130
 131                 fires_interp = np.reshape(fires_interp,lons_interp.shape)
 132
 133                 plot = False
 134                 if plot:
 135                         print '> Plotting the data...'
 136                         #number of elements in arrays to plot
 137                         nums = 100
 138                         nums1 = 50
 139                         fig1, (ax1,ax2) = plt.subplots(nrows = 2, ncols =1)
 140                         ax1.pcolor(lons_interp[0::nums1],lats_interp[0::nums1],fires_interp[0::nums1])
 141                         ax2.scatter(lons[0::nums],lats[0::nums],c=fires[0::nums],edgecolors='face')
 142                         plt.show()
 143
 144                 print '> Saving the data...'
 145                 t1 = time()
 146                 result = {'data': fires_interp.astype(np.int8),'geotransform':geotransform}
 147                 sio.savemat(file_name, mdict = result)
 148                 t2 = time()
 149                 print 'elapsed time: ',abs(t2-t1),'s'
 150
 151                 print '> File saved as ',file_name
 152         else:
 153                 print '> File already exists ',file_name
 154
 155         tg2 = time()
 156         print '> Elapsed time for the granule: ',abs(tg2-tg1),'s'
 157         print ''
 158         sys.stdout.flush()
 159
 160 t_final = time()
 161 print '>> Total elapsed time: ',abs(t_final-t_init),'s <<'
 162
 163
 164