source: issm/oecreview/Archive/16554-17801/ISSM-17696-17697.diff@ 17802

../trunk-jpl/src/m/interp/interp.py

@@ -1 +1 @@
 # module for inperpolating/smoothing data
 import numpy as npy
 from scipy.interpolate import CloughTocher2DInterpolator, Rbf
+from scipy.spatial import cKDTree
 try:
         import matplotlib.pyplot as plt
 except ImportError:
         print 'could not import matplotlib, no plotting functions enabled.\
                         Set plotonly=False in function call'
 
-def MeshSplineToMesh2d(x,y,data,xi,yi,tol=1e-6,**kwargs):#{{{
+def MeshSplineToMesh2d(x,y,data,xi,yi,tol=1e-6,fill_nans=False,**kwargs):#{{{
         '''
         Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.
         The interpolant is guaranteed to be continuously differentiable,
@@ -20 +21 @@
         data:                   data to be interpolated (same length as x,y)
         xi,yi:          coordintes to interpolate data onto
         tol:                    tolerance for gradient estimation (default 1e-6)
+        fill_nans:      fill nan's (holes) in data using the spline fit? 
         **kwargs:       optional keywork arguments:
                                         maxiter: maximum iterations in gradient estimation
         
@@ -56 +58 @@
         # mask out any nan's in the data and corresponding coordinate points
         mask=npy.isnan(subdata)
         ind=npy.nonzero(mask)[0]
-        if len(ind):
-                print "WARNING: input data grid contains nan values. Spline interpolation may be questionable."
+        if len(ind) and fill_nans:
+                print "WARNING: filling nans using spline fit through good data points,\
+                                which may or may not be appropriate. Check results carefully."
         subdata=npy.delete(subdata,ind)
         points=npy.delete(points,ind,axis=0)
 
@@ -67 +70 @@
                 spline=CloughTocher2DInterpolator(points,subdata,tol)
 
         interpdata=spline(xi,yi)
+        
+        if not fill_nans:
+                # identify nan's in xi,yi using nearest neighbors
+                xyinterp=npy.dstack([xi,yi])[0]
+                xg,yg=npy.meshgrid(subx,suby)
+                xydata=npy.dstack([xg.ravel(),yg.ravel()])[0]
+                tree=cKDTree(xydata)
+                nearest=tree.query(xyinterp)[1]
+                pos=npy.nonzero(npy.isnan(subdata[nearest]))
+                interpdata[pos]=subdata[nearest][pos]
 
         return interpdata
 #}}}
-def GridSplineToMesh2d(x,y,data,xi,yi,default_value=npy.nan,plotonly=False):#{{{
+def GridSplineToMesh2d(x,y,data,xi,yi,default_value=npy.nan,plotonly=False,fill_nans=False):#{{{
         '''
         python analog to InterpFromGridToMesh.  This routine uses
         scipy.interpolate.CloughTocher2dInterpolator to create a bivariate spline
@@ -89 +102 @@
         default_value:  default value if points lie outside the convex hull of input
                                                 points (defaults to nan if not specified)
         plotonly:               plot the data to be interpolated using imshow (useful for
-                                                identifying holes in data and problems with interpolation)
+        fill_nans:              fill nan's (holes) in data using the spline fit? 
 
         Usage:
-                interpdata=GridToMesh(x,y,data,xi,yi,default_value=npy.nan,plotonly=False)
+                interpdata=GridToMesh(x,y,data,xi,yi,default_value=npy.nan,plotonly=False,fill_nans=False)
 
         Examples:
                 interpdata=GridToMesh(x_m,y_m,data,md.mesh.x,md.mesh.y,0)
@@ -129 +142 @@
         # mask out any nan's in the data and corresponding coordinate points
         mask=npy.isnan(flatsubdata)
         ind=npy.nonzero(mask)[0]
-        if len(ind):
-                print "WARNING: input data grid contains nan values. Spline interpolation may be questionable."
-        flatsubdata=npy.delete(flatsubdata,ind)
-        points=npy.delete(points,ind,axis=0)
+        if len(ind) and fill_nans:
+                print "WARNING: filling nans using spline fit through good data points,\
+                                which may or may not be appropriate. Check results carefully."
+        goodsubdata=npy.delete(flatsubdata,ind)
+        goodpoints=npy.delete(points,ind,axis=0)
 
         # create spline and index spline at mesh points
-        spline=CloughTocher2DInterpolator(points,flatsubdata)
+        spline=CloughTocher2DInterpolator(goodpoints,goodsubdata)
         interpdata=spline(xi,yi)
 
+        if not fill_nans:
+                # identify nan's in xi,yi using nearest neighbors
+                xyinterp=npy.dstack([xi,yi])[0]
+                xydata=npy.dstack([xg.ravel(),yg.ravel()])[0]
+                tree=cKDTree(xydata)
+                nearest=tree.query(xyinterp)[1]
+                pos=npy.nonzero(npy.isnan(flatsubdata[nearest]))
+                interpdata[pos]=flatsubdata[nearest][pos]
+
         return interpdata
 #}}}
 def RadialInterp(x,y,data,**kwargs):#{{{