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         return ll
 143
 144 def neighbor_indices_ellipse(lons,lats,lon,lat,scan,track,mm=1):
 145         """
 146     Computes all the neighbor indices in an ellipse of radius the size of the pixel
 147
 148     :param lons: one dimensional array of the fire mesh longitud coordinates
 149     :param lats: one dimensional array of the fire mesh latitude coordinates
 150     :param lon: one dimensional array of the fire detection longitud coordinates
 151     :param lat: one dimensional array of the fire detection latitude coordinates
 152     :param scan: one dimensional array of the fire detection along-scan sizes
 153     :param track: one dimensional array of the fire detection along-track sizes
 154     :param mm: constant of magnitude of the ellipse, default mm=1 (ellipse inscribed in the pixel)
 155     :return ll: returns a one dimensional logical array with the indices in the fire mesh including pixel neighbours
 156
 157     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 158     Angel Farguell (angel.farguell@gmail.com), 2018-10-18
 159     """
 160         R=6378   # Earth radius
 161         km_lat=180/(np.pi*R)  # 1km in degrees latitude
 162         km_lon=km_lat/np.cos(lat*np.pi/180)  # 1 km in longitude
 163         # pixel sizes in degrees
 164         sqlat=km_lat*track/2
 165         sqlon=km_lon*scan/2
 166         # creating the coordinates in the center of the ellipse
 167         C=[[(np.array(lons)-lon[k])/sqlon[k],(np.array(lats)-lat[k])/sqlat[k]] for k in range(len(lat))]
 168         # creating a logical array of indices in the fire mesh with the intersection of all the cases
 169         ll=np.sum([np.array((np.square(c[0])+np.square(c[1]))<=mm) for c in C],axis=0).astype(bool)
 170         return ll
 171
 172 def neighbor_indices_ball(tree,indices,shape,d=2):
 173         """
 174     Computes all the neighbor indices from an indice list in a grid of indices defined through a cKDTree
 175
 176     :param indices: list of coordinates in a 1D array
 177     :param shape: array of two elements with satellite grid size
 178     :param d: optional, distance of the neighbours computed as: sqrt(2*d**2)
 179     :return ind: returns a numpy array with the indices and the neighbor indices in 1D array respect to the grid used in the tree
 180
 181     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 182     Angel Farguell (angel.farguell@gmail.com), 2018-09-20
 183     """
 184     # Width of the 2D grid
 185         w=shape[0]
 186         # 2D indices of the 1D indices
 187         I=[[np.mod(ind,w),ind/w] for ind in indices]
 188         # Radius to look for
 189         radius=np.sqrt(2*d**2)
 190         # Request all the neighbors in a radius "radius"
 191         N=tree.query_ball_point(I,r=radius)
 192         # Return an unique and sorted array of 1D indices (indices pointing to the grid used for the tree)
 193         ind=sorted(np.unique(np.concatenate(N)))
 194         return ind
 195
 196 if __name__ == "__main__":
 197         t_init = time()
 198         # Initialization of grids
 199         N=100
 200         (dx1,dx2)=(1,1)
 201         (dy1,dy2)=(3,3)
 202         x=np.arange(0,N,dx1)
 203         lons=np.repeat(x[np.newaxis,:],x.shape[0], axis=0)
 204         x=np.arange(0,N,dx2)
 205         lats=np.repeat(x[np.newaxis,:],x.shape[0], axis=0).T
 206         bounds=[lons.min(),lons.max(),lats.min(),lats.max()]
 207         print 'bounds'
 208         print bounds
 209         print 'dim_mesh=(%d,%d)' % lons.shape
 210         y=np.arange(-N*1.432,2.432*N,dy1)
 211         lon=np.repeat(y[np.newaxis,:],y.shape[0], axis=0)
 212         y=np.arange(-N*1.432,2.432*N,dy2)
 213         lat=np.repeat(y[np.newaxis,:],y.shape[0], axis=0)
 214         print 'dim_data=(%d,%d)' % (lon.shape[0], lat.shape[0])
 215
 216         # Result by Euclidean distance
 217         print '>>Euclidean approax<<'
 218         (rlon,rlat)=nearest_euclidean(lon,lat,lons,lats,bounds)
 219         rlon=np.reshape(rlon,np.prod(lon.shape))
 220         rlat=np.reshape(rlat,np.prod(lat.shape))
 221         vlonlatm=np.column_stack((rlon,rlat))
 222         print vlonlatm
 223         t_final = time()
 224         print 'Elapsed time: %ss.' % str(t_final-t_init)
 225
 226         # Result by scipy.spatial.cKDTree function
 227         vlons=np.reshape(lons,np.prod(lons.shape))
 228         vlats=np.reshape(lats,np.prod(lats.shape))
 229         vlonlats=np.column_stack((vlons,vlats))
 230         print vlonlats
 231         stree=spatial.cKDTree(vlonlats)
 232         (inds,K)=nearest_scipy(lon,lat,stree,bounds)
 233         vlonlatm2=np.empty((np.prod(lon.shape),2))
 234         vlonlatm2[:]=np.nan
 235         vlons=np.reshape(lons,np.prod(lons.shape))
 236         vlats=np.reshape(lats,np.prod(lats.shape))
 237         vlonlats=np.column_stack((vlons,vlats))
 238         vlonlatm2[K]=vlonlats[inds]
 239         print '>>cKDTree<<'
 240         print vlonlatm2
 241         t_ffinal = time()
 242         print 'Elapsed time: %ss.' % str(t_ffinal-t_final)
 243
 244         # Comparison
 245         print 'Same results?'
 246         print (np.isclose(vlonlatm,vlonlatm2) | np.isnan(vlonlatm) | np.isnan(vlonlatm2)).all()
 247
 248         # Testing neighbor indices
 249         N=100
 250         shape=(250,250)
 251         ind=sorted(np.unique([randint(0,np.prod(shape)-1) for i in range(0,N)]))
 252         #ind=[0,3,shape[0]/2+shape[1]/2*shape[0],np.prod(shape)-1]
 253         t_init = time()
 254         ne=neighbor_indices(ind,shape,d=8)
 255         t_final = time()
 256         print '1D neighbors:'
 257         print 'elapsed time: %ss.' % str(t_final-t_init)
 258         grid=np.array(list(itertools.product(np.array(range(0,shape[0])),np.array(range(0,shape[1])))))
 259         tree=spatial.cKDTree(grid)
 260         t_init = time()
 261         kk=neighbor_indices_ball(tree,ind,shape,d=8)
 262         t_final = time()
 263         nse=[x[0]+x[1]*shape[0] for x in grid[kk]]
 264         print '1D neighbors scipy:'
 265         print 'elapsed time: %ss.' % str(t_final-t_init)
 266         print 'Difference'
 267         print np.setdiff1d(ne,nse)
 268         print np.setdiff1d(nse,ne)