detection/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 print('Loading data...')
   8 data,fxlon,fxlat,time_num=sl.load('data')
   9 prefixes={'MOD': 'MOD14', 'MYD': 'MYD14', 'VNP': 'NPP_VAF_L2'}
  10
  11 # constants for geotransform
  12 res = 0.01   # resolution
  13 rot = 0.0    # rotation (not currently supported)
  14
  15 # radius for the tree
  16 radius = 0.05
  17
  18 for gran in list(data):
  19         splitted=gran.split('_')
  20         prefix=splitted[0]
  21         date=splitted[1][3:8]+splitted[2]+'00'
  22         file_name = prefixes[prefix]+'.'+date+'.tif.mat'
  23         #print(name)
  24         #G=data[gran]
  25
  26         bounds = [data[gran].lon.min(),data[gran].lon.max(),data[gran].lat.min(),data[gran].lat.max()]
  27         lons_interp = np.arange(bounds[0],bounds[1],res)
  28         lats_interp = np.arange(bounds[2],bounds[3],res)
  29         lons_interp,lats_interp = np.meshgrid(lons_interp,lats_interp)
  30
  31         geotransform = [bounds[0],res,rot,bounds[3],rot,res]
  32         #print(geotransform)
  33
  34         # flatten to 1d the arrays
  35         lons = np.reshape(data[gran].lon,np.prod(data[gran].lon.shape))
  36         lats = np.reshape(data[gran].lat,np.prod(data[gran].lat.shape))
  37         fires = np.reshape(data[gran].fire,np.prod(data[gran].fire.shape)).astype(np.int8)
  38
  39         # making tree
  40
  41         tree = spatial.cKDTree(np.column_stack((lons,lats)))
  42         glons = np.reshape(lons_interp,np.prod(lons_interp.shape))
  43         glats = np.reshape(lats_interp,np.prod(lats_interp.shape))
  44         indexes = np.array(tree.query_ball_point(np.column_stack((glons,glats)),radius))
  45         filtered_indexes = np.array([index[0] if len(index) > 0 else np.nan for index in indexes])
  46         fire1d = [fires[int(ii)] if not np.isnan(ii) else 0 for ii in filtered_indexes]
  47         fires_interp = np.reshape(fire1d,lons_interp.shape)
  48
  49         #coordinates = np.array((lons,lats)).T
  50         #print(coordinates.shape)
  51         #print(coordinates)
  52         #print(lons_interp.shape)
  53         '''
  54         print(lons_interp)
  55         print(lats_interp.shape)
  56         print(lats_interp)
  57         '''
  58         #fires_interp = griddata(coordinates,fires,(lons_interp,lats_interp),method='nearest',fill_value = np.nan)
  59         print('Done interpolating, now plotting')
  60
  61         print('Nans in fire: ',np.isnan(fires_interp).sum())
  62
  63         plot = False
  64         if plot:
  65                 #number of elements in arrays to plot
  66                 nums = 100
  67                 nums1 = 50
  68                 fig1, (ax1,ax2) = plt.subplots(nrows = 2, ncols =1)
  69                 ax1.pcolor(lons_interp[0::nums1],lats_interp[0::nums1],fires_interp[0::nums1])
  70                 #plt.colorbar()
  71                 ax2.scatter(lons[0::nums],lats[0::nums],c=fires[0::nums],edgecolors='face')
  72                 plt.show()
  73         result = {'data': fires_interp.astype(np.int8),'geotransform':geotransform}
  74         sio.savemat(file_name,mdict = result)
  75
  76         print('File saved as ',file_name)
  77
  78
  79
  80
  81
  82