Index: /issm/trunk-jpl/src/m/interp/interp.py
===================================================================
--- /issm/trunk-jpl/src/m/interp/interp.py	(revision 17708)
+++ /issm/trunk-jpl/src/m/interp/interp.py	(revision 17709)
@@ -116,8 +116,8 @@
 	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')
+	if len(x) != data.shape[1]+1 and len(x) != data.shape[1]:
+		raise ValueError('x should have same length as ncols(data) or ncols(data)+1')
+	if len(y) != data.shape[0]+1 and len(y) != data.shape[0]:
+		raise ValueError('y should have same length as nrows(data) or nrows(data)+1')
 	
 	# create sub-grid that just covers the limits of xi and yi
@@ -126,10 +126,25 @@
 	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]
-
+
+	# TODO create grid differently depending on whether data is defined at x,y
+	# or at the center of a grid cell with corner coordinates defined by xi,yi
 	# create points array and flattened data array
-	xg,yg=npy.meshgrid(x[xind],y[yind])
+	if len(x)==data.shape[1] and len(y)==data.shape[0]:
+		print 'x,y taken to define the center of data grid cells'
+		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]
+		xg,yg=npy.meshgrid(x[xind],y[yind])
+		subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1]
+	elif len(x)==data.shape[1]+1 and len(y)==data.shape[0]+1:
+		print 'x,y taken to define the corners of data grid cells'
+		xcenter=npy.fromiter(((x[i]+x[i+1])/2 for i in range(len(x)-1)),npy.float)
+		ycenter=npy.fromiter(((y[i]+y[i+1])/2 for i in range(len(y)-1)),npy.float)
+		xind=npy.nonzero(npy.logical_and(xcenter>xlim[0],xcenter<xlim[1]))[0]
+		yind=npy.nonzero(npy.logical_and(ycenter>ylim[0],ycenter<ylim[1]))[0]
+		xg,yg=npy.meshgrid(xcenter[xind],ycenter[yind])
+		subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1]
+	else:
+		raise ValueError('x and y have inconsistent sizes: both should have length ncols(data)/nrows(data) or ncols(data)+1/nrows(data)+1')
+
 	points=npy.array([xg.ravel(),yg.ravel()]).T
 	flatsubdata=subdata.ravel()