interpolation.py

   1 import warnings
   2 warnings.filterwarnings("ignore")
   3 import numpy as np
   4 from time import time
   5 from datetime import datetime
   6 from scipy import spatial
   7 import itertools
   8 from random import randint
   9
  10 def sort_dates(data):
  11         """
  12     Sorting a dictionary depending on the time number in seconds from January 1, 1970
  13
  14     :param data: dictionary of granules where each granule has a time_num key
  15     :return: an array of dictionaries ordered by time_num key (seconds from January 1, 1970)
  16
  17     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  18     Angel Farguell (angel.farguell@gmail.com), 2018-09-20
  19     """
  20         return sorted(data.iteritems(), key=lambda x: x[1]['time_num'])
  21
  22 def nearest_euclidean(lon,lat,lons,lats,bounds):
  23         """
  24     Returns the longitude and latitude arrays interpolated using Euclidean distance
  25
  26     :param lon: 2D array of longitudes to look the nearest neighbours
  27     :param lat: 2D array of latitudes to look the nearest neighbours
  28     :param lons: 2D array of longitudes interpolating to
  29     :param lats: 2D array of latitudes interpolating to
  30     :param bounds: array of 4 bounding boundaries where to interpolate: [minlon maxlon minlat maxlat]
  31     :return ret: tuple with a 2D array of longitudes and 2D array of latitudes interpolated from (lon,lat) to (lons,lats)
  32
  33     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  34     Angel Farguell (angel.farguell@gmail.com), 2018-09-20
  35     """
  36         vlon=np.reshape(lon,np.prod(lon.shape))
  37         vlat=np.reshape(lat,np.prod(lat.shape))
  38         rlon=np.zeros(vlon.shape)
  39         rlat=np.zeros(vlat.shape)
  40         for k in range(0,len(vlon)):
  41                 if (vlon[k]>bounds[0]) and (vlon[k]<bounds[1]) and (vlat[k]>bounds[2]) and (vlat[k]<bounds[3]):
  42                         dist=np.square(lons-vlon[k])+np.square(lats-vlat[k])
  43                         ii,jj = np.unravel_index(dist.argmin(),dist.shape)
  44                         rlon[k]=lons[ii,jj]
  45                         rlat[k]=lats[ii,jj]
  46                 else:
  47                         rlon[k]=np.nan;
  48                         rlat[k]=np.nan;
  49         ret=(np.reshape(rlon,lon.shape),np.reshape(rlat,lat.shape))
  50         return ret
  51
  52
  53 def nearest_scipy(lon,lat,stree,bounds):
  54         """
  55     Returns the coordinates interpolated from (lon,lat) to (lons,lats) and the value of the indices where NaN values
  56
  57     :param lon: 2D array of longitudes to look the nearest neighbours
  58     :param lat: 2D array of latitudes to look the nearest neighbours
  59     :param lons: 2D array of longitudes interpolating to
  60     :param lats: 2D array of latitudes interpolating to
  61     :param dub: optional, distance upper bound to look for the nearest neighbours
  62     :return ret: A tuple with a 2D array of coordinates interpolated from (lon,lat) to (lons,lats) and the value of the indices where NaN values
  63
  64     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  65     Angel Farguell (angel.farguell@gmail.com), 2018-09-20
  66     """
  67         vlon=np.reshape(lon,np.prod(lon.shape))
  68         vlat=np.reshape(lat,np.prod(lat.shape))
  69         vlonlat=np.column_stack((vlon,vlat))
  70         M=(vlon>bounds[0])*(vlon<bounds[1])*(vlat>bounds[2])*(vlat<bounds[3])
  71         vlonlat=vlonlat[M]
  72         inds=np.array(stree.query(vlonlat)[1])
  73         ret=(inds,M)
  74         return ret
  75
  76 def distance_upper_bound(dx,dy):
  77         """
  78     Computes the distance upper bound
  79
  80     :param dx: array of two elements with fire mesh grid resolutions
  81     :param dy: array of two elements with satellite grid resolutions
  82     :return d: distance upper bound to look for the nearest neighbours
  83
  84     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  85     Angel Farguell (angel.farguell@gmail.com), 2018-09-20
  86     """
  87         rx=np.sqrt(dx[0]**2+dx[1]**2)
  88         ry=np.sqrt(dy[0]**2+dy[1]**2)
  89         d=max(rx,ry)
  90         return d
  91
  92 def neighbor_indices(indices,shape,d=2):
  93         """
  94     Computes all the neighbor indices from an indice list
  95
  96     :param indices: list of coordinates in a 1D array
  97     :param shape: array of two elements with satellite grid size
  98     :param d: optional, distance of the neighbours
  99     :return ind: Returns a numpy array with the indices and the neighbor indices in 1D array
 100
 101     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 102     Angel Farguell (angel.farguell@gmail.com), 2018-09-20
 103     """
 104     # Width and Length of the 2D grid
 105         w=shape[0]
 106         l=shape[1]
 107         # 2D indices of the 1D indices
 108         I=[[np.mod(ind,w),ind/w] for ind in indices]
 109         # All the combinations (x,y) for all the neighbor points from x-d to x+d and y-d to y+d
 110         N=np.concatenate([np.array(list(itertools.product(range(max(i[0]-d,0),min(i[0]+d+1,w)),range(max(i[1]-d,0),min(i[1]+d+1,l))))) for i in I])
 111         # Recompute the 1D indices of the 2D coordinates inside the 2D domain
 112         ret=np.array([x[0]+w*x[1] for x in N])
 113         # Sort them and take each indice once
 114         ind=sorted(np.unique(ret))
 115         return ind
 116
 117 def neighbor_indices_pixel(lons,lats,lon,lat,scan,track):
 118         """
 119     Computes all the neighbor indices depending on the size of the pixel
 120
 121     :param lons: one dimensional array of the fire mesh longitud coordinates
 122     :param lats: one dimensional array of the fire mesh latitude coordinates
 123     :param lon: one dimensional array of the fire detection longitud coordinates
 124     :param lat: one dimensional array of the fire detection latitude coordinates
 125     :param scan: one dimensional array of the fire detection along-scan sizes
 126     :param track: one dimensional array of the fire detection along-track sizes
 127     :return ll: returns a one dimensional logical array with the indices in the fire mesh including pixel neighbours
 128
 129     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 130     Angel Farguell (angel.farguell@gmail.com), 2018-10-16
 131     """
 132         R=6378   # Earth radius
 133         km_lat=180/(np.pi*R)  # 1km in degrees latitude
 134         km_lon=km_lat/np.cos(lat*np.pi/180)  # 1 km in longitude
 135         # pixel sizes in degrees
 136         sqlat=km_lat*track/2
 137         sqlon=km_lon*scan/2
 138         # creating bounds for each fire detection (pixel vertexs)
 139         bounds=[[lon[k]-sqlon[k],lon[k]+sqlon[k],lat[k]-sqlat[k],lat[k]+sqlat[k]] for k in range(len(lat))]
 140         # creating a logical array of indices in the fire mesh with the intersection of all the cases
 141         ll=np.sum([np.array(np.logical_and(np.logical_and(np.logical_and(lons>b[0],lons<b[1]),lats>b[2]),lats<b[3])) for b in bounds],axis=0).astype(bool)
 142         if np.all(ll==False):
 143                 fll = np.array([False]*len(lons))
 144                 return fll
 145         return ll
 146
 147 def neighbor_indices_ellipse(lons,lats,lon,lat,scan,track,mm=1):
 148         """
 149     Computes all the neighbor indices in an ellipse of radius the size of the pixel
 150
 151     :param lons: one dimensional array of the fire mesh longitud coordinates
 152     :param lats: one dimensional array of the fire mesh latitude coordinates
 153     :param lon: one dimensional array of the fire detection longitud coordinates
 154     :param lat: one dimensional array of the fire detection latitude coordinates
 155     :param scan: one dimensional array of the fire detection along-scan sizes
 156     :param track: one dimensional array of the fire detection along-track sizes
 157     :param mm: constant of magnitude of the ellipse, default mm=1 (ellipse inscribed in the pixel)
 158     :return ll: returns a one dimensional logical array with the indices in the fire mesh including pixel neighbours
 159
 160     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 161     Angel Farguell (angel.farguell@gmail.com), 2018-10-18
 162     """
 163         R=6378   # Earth radius
 164         km_lat=180/(np.pi*R)  # 1km in degrees latitude
 165         km_lon=km_lat/np.cos(lat*np.pi/180)  # 1 km in longitude
 166         # pixel sizes in degrees
 167         sqlat=km_lat*track/2
 168         sqlon=km_lon*scan/2
 169         # creating the coordinates in the center of the ellipse
 170         C=[[(np.array(lons)-lon[k])/sqlon[k],(np.array(lats)-lat[k])/sqlat[k]] for k in range(len(lat))]
 171         # creating a logical array of indices in the fire mesh with the intersection of all the cases
 172         ll=np.sum([np.array((np.square(c[0])+np.square(c[1]))<=mm) for c in C],axis=0).astype(bool)
 173         if np.all(ll==False):
 174                 fll = np.array([False]*len(lons))
 175                 return fll
 176         return ll
 177
 178 def neighbor_indices_ball(tree,indices,shape,d=2):
 179         """
 180     Computes all the neighbor indices from an indice list in a grid of indices defined through a cKDTree
 181
 182     :param indices: list of coordinates in a 1D array
 183     :param shape: array of two elements with satellite grid size
 184     :param d: optional, distance of the neighbours computed as: sqrt(2*d**2)
 185     :return ind: returns a numpy array with the indices and the neighbor indices in 1D array respect to the grid used in the tree
 186
 187     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 188     Angel Farguell (angel.farguell@gmail.com), 2018-09-20
 189     """
 190     # Width of the 2D grid
 191         w=shape[0]
 192         # 2D indices of the 1D indices
 193         I=[[np.mod(ind,w),ind/w] for ind in indices]
 194         # Radius to look for
 195         radius=np.sqrt(2*d**2)
 196         # Request all the neighbors in a radius "radius"
 197         N=tree.query_ball_point(I,r=radius)
 198         # Return an unique and sorted array of 1D indices (indices pointing to the grid used for the tree)
 199         ind=sorted(np.unique(np.concatenate(N)))
 200         return ind
 201
 202 if __name__ == "__main__":
 203         t_init = time()
 204         # Initialization of grids
 205         N=100
 206         (dx1,dx2)=(1,1)
 207         (dy1,dy2)=(3,3)
 208         x=np.arange(0,N,dx1)
 209         lons=np.repeat(x[np.newaxis,:],x.shape[0], axis=0)
 210         x=np.arange(0,N,dx2)
 211         lats=np.repeat(x[np.newaxis,:],x.shape[0], axis=0).T
 212         bounds=[lons.min(),lons.max(),lats.min(),lats.max()]
 213         print 'bounds'
 214         print bounds
 215         print 'dim_mesh=(%d,%d)' % lons.shape
 216         y=np.arange(-N*1.432,2.432*N,dy1)
 217         lon=np.repeat(y[np.newaxis,:],y.shape[0], axis=0)
 218         y=np.arange(-N*1.432,2.432*N,dy2)
 219         lat=np.repeat(y[np.newaxis,:],y.shape[0], axis=0)
 220         print 'dim_data=(%d,%d)' % (lon.shape[0], lat.shape[0])
 221
 222         # Result by Euclidean distance
 223         print '>>Euclidean approax<<'
 224         (rlon,rlat)=nearest_euclidean(lon,lat,lons,lats,bounds)
 225         rlon=np.reshape(rlon,np.prod(lon.shape))
 226         rlat=np.reshape(rlat,np.prod(lat.shape))
 227         vlonlatm=np.column_stack((rlon,rlat))
 228         print vlonlatm
 229         t_final = time()
 230         print 'Elapsed time: %ss.' % str(t_final-t_init)
 231
 232         # Result by scipy.spatial.cKDTree function
 233         vlons=np.reshape(lons,np.prod(lons.shape))
 234         vlats=np.reshape(lats,np.prod(lats.shape))
 235         vlonlats=np.column_stack((vlons,vlats))
 236         print vlonlats
 237         stree=spatial.cKDTree(vlonlats)
 238         (inds,K)=nearest_scipy(lon,lat,stree,bounds)
 239         vlonlatm2=np.empty((np.prod(lon.shape),2))
 240         vlonlatm2[:]=np.nan
 241         vlons=np.reshape(lons,np.prod(lons.shape))
 242         vlats=np.reshape(lats,np.prod(lats.shape))
 243         vlonlats=np.column_stack((vlons,vlats))
 244         vlonlatm2[K]=vlonlats[inds]
 245         print '>>cKDTree<<'
 246         print vlonlatm2
 247         t_ffinal = time()
 248         print 'Elapsed time: %ss.' % str(t_ffinal-t_final)
 249
 250         # Comparison
 251         print 'Same results?'
 252         print (np.isclose(vlonlatm,vlonlatm2) | np.isnan(vlonlatm) | np.isnan(vlonlatm2)).all()
 253
 254         # Testing neighbor indices
 255         N=100
 256         shape=(250,250)
 257         ind=sorted(np.unique([randint(0,np.prod(shape)-1) for i in range(0,N)]))
 258         #ind=[0,3,shape[0]/2+shape[1]/2*shape[0],np.prod(shape)-1]
 259         t_init = time()
 260         ne=neighbor_indices(ind,shape,d=8)
 261         t_final = time()
 262         print '1D neighbors:'
 263         print 'elapsed time: %ss.' % str(t_final-t_init)
 264         grid=np.array(list(itertools.product(np.array(range(0,shape[0])),np.array(range(0,shape[1])))))
 265         tree=spatial.cKDTree(grid)
 266         t_init = time()
 267         kk=neighbor_indices_ball(tree,ind,shape,d=8)
 268         t_final = time()
 269         nse=[x[0]+x[1]*shape[0] for x in grid[kk]]
 270         print '1D neighbors scipy:'
 271         print 'elapsed time: %ss.' % str(t_final-t_init)
 272         print 'Difference'
 273         print np.setdiff1d(ne,nse)
 274         print np.setdiff1d(nse,ne)