Index: /issm/trunk-jpl/src/m/interp/interp.py
===================================================================
--- /issm/trunk-jpl/src/m/interp/interp.py	(revision 17696)
+++ /issm/trunk-jpl/src/m/interp/interp.py	(revision 17697)
@@ -2,4 +2,5 @@
 import numpy as npy
 from scipy.interpolate import CloughTocher2DInterpolator, Rbf
+from scipy.spatial import cKDTree
 try:
 	import matplotlib.pyplot as plt
@@ -8,5 +9,5 @@
 			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,4 +22,5 @@
 	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,6 +59,7 @@
 	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,8 +71,18 @@
 
 	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,8 +103,8 @@
 						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,12 +143,22 @@
 	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