Changeset 17709

Timestamp:

04/10/14 15:46:42 (11 years ago)

Author:

cborstad

Message:

BUG: define x,y coords in centers of grid cells if appropriate

File:

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

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

@@ -116 +116 @@
         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 +126 @@
         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()