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.005  # 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                 print 'cols:',len(lons_interp)
  77                 lats_interp = np.arange(bounds[2],bounds[3],res)
  78                 print 'rows:',len(lats_interp)
  79                 lons_interp,lats_interp = np.meshgrid(lons_interp,lats_interp)
  80                 glons = np.reshape(lons_interp,np.prod(lons_interp.shape))
  81                 glats = np.reshape(lats_interp,np.prod(lats_interp.shape))
  82
  83                 # initialize fire_interp
  84                 if not whole:
  85                         fires_interp = np.zeros(np.prod(lons_interp.shape))
  86
  87                 # creating geotransform array with necessary geolocation information
  88                 geotransform = [bounds[0]-res*.5,res,rot,bounds[2]-res*.5,rot,res]
  89
  90                 # approximation of the radius for the tree
  91                 dlon = abs(data[gran].lon[1,1]-data[gran].lon[0,0])
  92                 dlat = abs(data[gran].lat[1,1]-data[gran].lat[0,0])
  93                 if opt_rad == 'max':
  94                         radius = max(dlon,dlat)
  95                 elif opt_rad == 'mean':
  96                         radius = (dlon+dlat)*.5
  97                 elif opt_rad == 'min':
  98                         radius = min(dlon,dlat)
  99                 print 'radius =',radius
 100
 101                 # flatten granule data into 1d arrays
 102                 lons = np.reshape(data[gran].lon,np.prod(data[gran].lon.shape))
 103                 lats = np.reshape(data[gran].lat,np.prod(data[gran].lat.shape))
 104                 fires = np.reshape(data[gran].fire,np.prod(data[gran].fire.shape)).astype(np.int8)
 105                 if not whole:
 106                         M = (lons>fbounds[0])*(lons<fbounds[1])*(lats>fbounds[2])*(lats<fbounds[3])
 107                         Mg = (glons>fbounds[0])*(glons<fbounds[1])*(glats>fbounds[2])*(glats<fbounds[3])
 108                         lons = lons[M]
 109                         lats = lats[M]
 110                         fires = fires[M]
 111                         glons = glons[Mg]
 112                         glats = glats[Mg]
 113
 114                 if len(lons) > 0:
 115                         # Making tree
 116                         print '> Making the tree...'
 117                         t1 = time()
 118                         tree = spatial.cKDTree(np.column_stack((lons,lats)))
 119                         t2 = time()
 120                         print 'elapsed time: ',abs(t2-t1),'s'
 121                         print '> Interpolating the data...'
 122                         t1 = time()
 123                         indexes = np.array(tree.query_ball_point(np.column_stack((glons,glats)),radius,return_sorted=True))
 124                         t2 = time()
 125                         print 'elapsed time: ',abs(t2-t1),'s'
 126                         filtered_indexes = np.array([index[0] if len(index) > 0 else np.nan for index in indexes])
 127                         fire1d = [fires[int(ii)] if not np.isnan(ii) else 0 for ii in filtered_indexes]
 128                         if whole:
 129                                 fires_interp = np.reshape(fire1d,lons_interp.shape)
 130                         else:
 131                                 fires_interp[Mg] = fire1d
 132
 133                 fires_interp = np.reshape(fires_interp,lons_interp.shape)
 134
 135                 plot = False
 136                 if plot:
 137                         print '> Plotting the data...'
 138                         #number of elements in arrays to plot
 139                         nums = 100
 140                         nums1 = 50
 141                         fig1, (ax1,ax2) = plt.subplots(nrows = 2, ncols =1)
 142                         ax1.pcolor(lons_interp[0::nums1],lats_interp[0::nums1],fires_interp[0::nums1])
 143                         ax2.scatter(lons[0::nums],lats[0::nums],c=fires[0::nums],edgecolors='face')
 144                         plt.show()
 145
 146                 print '> Saving the data...'
 147                 t1 = time()
 148                 result = {'data': fires_interp.astype(np.int8),'geotransform':geotransform}
 149                 sio.savemat(file_name, mdict = result)
 150                 t2 = time()
 151                 print 'elapsed time: ',abs(t2-t1),'s'
 152
 153                 print '> File saved as ',file_name
 154         else:
 155                 print '> File already exists ',file_name
 156
 157         tg2 = time()
 158         print '> Elapsed time for the granule: ',abs(tg2-tg1),'s'
 159         print ''
 160         sys.stdout.flush()
 161
 162 t_final = time()
 163 print '>> Total elapsed time: ',abs(t_final-t_init),'s <<'
 164
 165
 166