Changeset 17697

Timestamp:

04/09/14 17:04:50 (11 years ago)

Author:

cborstad

Message:

CHG: by default, do not fill nan's using the spline fit, rather mark them using nearest neighbors

File:

• issm/trunk-jpl/src/m/interp/interp.py (modified) (7 diffs)

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

@@ -2 +2 @@
 import numpy as npy
 from scipy.interpolate import CloughTocher2DInterpolator, Rbf
+from scipy.spatial import cKDTree
 try:
         import matplotlib.pyplot as plt
@@ -8 +9 @@
                         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.
@@ -21 +22 @@
         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
@@ -57 +59 @@
         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)
@@ -68 +71 @@
 
         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
@@ -90 +103 @@
                                                 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:
@@ -130 +143 @@
         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