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):#{{{