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

Index: ../trunk-jpl/src/m/interp/interp.py
===================================================================
--- ../trunk-jpl/src/m/interp/interp.py (revision 17696)
+++ ../trunk-jpl/src/m/interp/interp.py (revision 17697)
@@ -1,13 +1,14 @@
 # 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,6 +21,7 @@
        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,8 +58,9 @@
        # 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,10 +70,20 @@
                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,10 +102,10 @@
        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,15 +142,25 @@
        # 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):#{{{