Changeset 17598

Timestamp:

03/28/14 12:18:11 (11 years ago)

Author:

cborstad

Message:

CHG: subset the input data for better efficiency, mask out nan's since they corrupt the spline calculation

File:

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

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

@@ -36 +36 @@
         assert len(kwargs)==0, 'error, unexpected or misspelled kwargs'
 
-        points=npy.array([x,y]).T
+        # create sub-vectors that just cover the limits of xi and yi
+        dx=x[1]-x[0]
+        dy=y[1]-y[0]
+        xlim=[min(xi)-dx,max(xi)+dx]
+        ylim=[min(yi)-dy,max(yi)+dy]
+        xflag=npy.logical_and(x>xlim[0],x<xlim[1])
+        yflag=npy.logical_and(y>ylim[0],y<ylim[1])
+        bothind=npy.nonzero(npy.logical_and(xflag,yflag))
+        subdata=data[bothind]
+        subx=x[bothind]
+        suby=y[bothind]
+        points=npy.array([subx,suby]).T
+
+        # 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."
+        subdata=npy.delete(subdata,ind)
+        points=npy.delete(points,ind,axis=0)
 
         if maxiter:
-                curvy=CloughTocher2DInterpolator(points,data,tol,maxiter=maxiter)
+                spline=CloughTocher2DInterpolator(points,subdata,tol,maxiter=maxiter)
         else:
-                curvy=CloughTocher2DInterpolator(points,data,tol)
+                spline=CloughTocher2DInterpolator(points,subdata,tol)
 
-        interpdata=curvy(xi,yi)
+        interpdata=spline(xi,yi)
 
         return interpdata
@@ -99 +118 @@
         mask=npy.isnan(flatsubdata)
         ind=npy.nonzero(mask)[0]
-        if float(len(ind))/float(len(flatsubdata)) > 0.1:
-                print "WARNING: over 10% of data grid contains nan values. Spline interpolation may be questionable."
+        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)
 
         # create spline and index spline at mesh points
-        spl=CloughTocher2DInterpolator(points,flatsubdata)
-        interpdata=spl(xi,yi)
+        spline=CloughTocher2DInterpolator(points,flatsubdata)
+        interpdata=spline(xi,yi)
 
         return interpdata