Index: /issm/trunk-jpl/src/m/interp/interp.py
===================================================================
--- /issm/trunk-jpl/src/m/interp/interp.py	(revision 17596)
+++ /issm/trunk-jpl/src/m/interp/interp.py	(revision 17596)
@@ -0,0 +1,169 @@
+# module for inperpolating/smoothing data
+import numpy as npy
+from scipy.interpolate import CloughTocher2DInterpolator, Rbf
+
+def MeshSplineToMesh2d(x,y,data,xi,yi,tol=1e-6,**kwargs):#{{{
+	'''
+	Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.
+	The interpolant is guaranteed to be continuously differentiable,
+	and the gradients are chosen such that the curvature of the interpolant
+	is approximately minimized.
+
+	Uses scipy.interpolate.CloughTocher2DInterpolator
+
+	x,y:			data point coordinates
+	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)
+	**kwargs:	optional keywork arguments:
+					maxiter: maximum iterations in gradient estimation
+	
+	Returns interpolated data at given x,y coordinates.
+
+	Usage:
+		interpdata=CloughToucher2d(x,y,data)
+
+	Examples:
+		interpdata=CloughToucher2d(md.mesh.x,md.mesh.y,data)
+		interpdata=CloughToucher2d(md.mesh.x,md.mesh.y,data,tol=1e-3,maxiter=100)
+	'''
+
+	# unpack kwargs
+	maxiter=kwargs.pop('maxiter',None)
+	if 'maxiter' in kwargs: del kwargs['maxiter']
+	if maxiter:
+		assert type(maxiter)==int, 'error, maxiter should be an integer'
+	assert len(kwargs)==0, 'error, unexpected or misspelled kwargs'
+
+	points=npy.array([x,y]).T
+
+	if maxiter:
+		curvy=CloughTocher2DInterpolator(points,data,tol,maxiter=maxiter)
+	else:
+		curvy=CloughTocher2DInterpolator(points,data,tol)
+
+	interpdata=curvy(xi,yi)
+
+	return interpdata
+#}}}
+def GridSplineToMesh2d(x,y,data,xi,yi,default_value=npy.nan):#{{{
+	'''
+	python analog to InterpFromGridToMesh.  This routine uses
+	scipy.interpolate.CloughTocher2dInterpolator to create a bivariate spline
+	interpolation of the input data and then return values of the spline
+	on the x,y coordinates of the model mesh.  The interpolant is piece-wise
+	cubic, C1 smooth (continuously differentiable) and has approximately 
+	minimized curvature.  See "help(scipy.interpolate.CloughTocher2dInterpolator)"
+	for more information on the routine.
+
+	NOTE: this routine will not be appropriate if there are large holes (nan's) in 
+	the input data.  A non-spline interpolation scheme should be used in that case.
+
+	x,y:				vectors defining the coordinates of the input data
+	data:				2D array of input data
+	xi,yi:			x and y coordinates to be interpolated onto
+	default_value:	default value if points lie outside the convex hull of input
+						points (defaults to nan if not specified)
+
+	Usage:
+		interpdata=GridToMesh(x,y,data,xi,yi,default_value=npy.nan)
+
+	Examples:
+		interpdata=GridToMesh(x_m,y_m,data,md.mesh.x,md.mesh.y,0)
+	'''
+
+	if npy.ndim(x)==2:
+		x=x.reshape(-1,)
+	if npy.ndim(y)==2:
+		y=y.reshape(-1,)
+	if len(x) > data.shape[1]+1:
+		raise ValueError('x should have same length as ncols(data)+1')
+	if len(y) > data.shape[0]+1:
+		raise ValueError('y should have same length as nrows(data)+1')
+	
+	# create sub-grid that just covers 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]
+	xind=npy.nonzero(npy.logical_and(x>xlim[0],x<xlim[1]))[0]
+	yind=npy.nonzero(npy.logical_and(y>ylim[0],y<ylim[1]))[0]
+	subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1]
+
+	# create points array and flattened data array
+	xg,yg=npy.meshgrid(x[xind],y[yind])
+	points=npy.array([xg.ravel(),yg.ravel()]).T
+	flatsubdata=subdata.ravel()
+
+	# mask out any nan's in the data and corresponding coordinate points
+	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."
+	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)
+
+	return interpdata
+#}}}
+def RadialInterp(x,y,data,**kwargs):#{{{
+	'''
+	Interpolation using a radial basis function in 2 or 3 dimensions.
+	Useful for smoothing input data after interpolation.
+
+	Uses scipy.interpolate.Rbf
+
+	x,y:			data point coordinates
+	data:			data to be interpolated (same length as x,y)
+
+		function: form of radial basis function for interpolation:
+			'multiquadric': sqrt((r/self.epsilon)**2 + 1) (default)
+			'inverse': 1.0/sqrt((r/self.epsilon)**2 + 1)
+			'gaussian': exp(-(r/self.epsilon)**2)
+			'linear': r
+			'cubic': r**3
+			'quintic': r**5
+			'thin_plate': r**2 * log(r)
+		epsilon: adjustable constant for scaling radial distance.  Defaults to 
+					approximate average distance between nodes.
+		smooth: float>0, adjusts the amount of smoothing applied.  Defaults to 0,
+					such that the function always passes through nodal points.
+		z:	coordinate array if interpolating in 3 dimensions
+
+	Usage:
+		interpdata=RadialInterp(x,y,data,**kwargs)
+
+	Examples:
+		interpdata=RadialInterp(md.mesh.x,md.mesh.y,data)
+		interpdata=RadialInterp(md.mesh.x,md.mesh.y,data,function='gaussian',epsilon=100,smooth=1)
+	'''
+
+	# unpack kwargs
+	function=kwargs.pop('function','multiquadric')
+	if 'function' in kwargs: del kwargs['function']
+	epsilon=kwargs.pop('epsilon',None)
+	if 'epsilon' in kwargs: del kwargs['epsilon']
+	smooth=kwargs.pop('smooth',0)
+	if 'smooth' in kwargs: del kwargs['smooth']
+	z=kwargs.pop('z',None)
+	if 'z' in kwargs: del kwargs['z']
+	assert len(kwargs)==0, 'error, unexpected or misspelled kwargs'
+
+	if z:
+		if epsilon:
+			rbfi=Rbf(x,y,z,data,function=function,smooth=smooth,epsilon=epsilon)
+		else:
+			rbfi=Rbf(x,y,z,data,function=function,smooth=smooth)
+		interpdata=rbfi(x,y,z)
+	else:
+		if epsilon:
+			rbfi=Rbf(x,y,data,function=function,smooth=smooth,epsilon=epsilon)
+		else:
+			rbfi=Rbf(x,y,data,function=function,smooth=smooth)
+		interpdata=rbfi(x,y)
+	
+	return interpdata
+#}}}