test if p.lon/p.lat are between min/max of grid

Author: delandmeterp

parcels/field.py

@@ -403,9 +403,13 @@ def search_indices_rectilinear(self, x, y, z, ti=-1, time=-1, search2D=False):
         grid = self.grid
         xi = yi = -1
 
-        if grid.mesh is not 'spherical':
-            if x < grid.lon[0] or x > grid.lon[-1]:
+        if not grid.zonal_periodic:
+            if x < grid.lonlat_minmax[0] or x > grid.lonlat_minmax[1]:
                 raise FieldSamplingError(x, y, z, field=self)
+        if y < grid.lonlat_minmax[2] or y > grid.lonlat_minmax[3]:
+            raise FieldSamplingError(x, y, z, field=self)
+
+        if grid.mesh is not 'spherical':
             lon_index = grid.lon < x
             if lon_index.all():
                 xi = len(grid.lon) - 2
@@ -425,11 +429,6 @@ def search_indices_rectilinear(self, x, y, z, ti=-1, time=-1, search2D=False):
                 lon_fixed[indices.argmin():] += 360
             if x < lon_fixed[0]:
                 lon_fixed -= 360
-            if not grid.zonal_periodic:
-                if (grid.lon[0] < grid.lon[-1]) and (x < grid.lon[0] or x > grid.lon[-1]):
-                    raise FieldSamplingError(x, y, z, field=self)
-                elif (grid.lon[0] >= grid.lon[-1]) and (x < grid.lon[0] and x > grid.lon[-1]):
-                    raise FieldSamplingError(x, y, z, field=self)
 
             lon_index = lon_fixed < x
             if lon_index.all():
@@ -444,8 +443,6 @@ def search_indices_rectilinear(self, x, y, z, ti=-1, time=-1, search2D=False):
                 xi += 1
                 xsi = (x-lon_fixed[xi]) / (lon_fixed[xi+1]-lon_fixed[xi])
 
-        if y < grid.lat[0] or y > grid.lat[-1]:
-            raise FieldSamplingError(x, y, z, field=self)
         lat_index = grid.lat < y
         if lat_index.all():
             yi = len(grid.lat) - 2
@@ -485,12 +482,13 @@ def search_indices_curvilinear(self, x, y, z, xi, yi, ti=-1, time=-1, search2D=F
                          [1, -1, 1, -1]])
         maxIterSearch = 1e6
         it = 0
-        if (not grid.zonal_periodic) or grid.mesh == 'flat':
-            if (grid.lon[0, 0] < grid.lon[0, -1]) and (x < grid.lon[0, 0] or x > grid.lon[0, -1]):
-                raise FieldSamplingError(x, y, z, field=self)
-            elif (grid.lon[0, 0] >= grid.lon[0, -1]) and (x < grid.lon[0, 0] and x > grid.lon[0, -1]):
-                raise FieldSamplingError(x, y, z, field=self)
-        if y < np.min(grid.lat) or y > np.max(grid.lat):
+        if not grid.zonal_periodic:
+            if x < grid.lonlat_minmax[0] or x > grid.lonlat_minmax[1]:
+                if grid.lon[0, 0] < grid.lon[0, -1]:
+                    raise FieldSamplingError(x, y, z, field=self)
+                elif x < grid.lon[0, 0] and x > grid.lon[0, -1]:  # This prevents from crashing in [160, -160]
+                    raise FieldSamplingError(x, y, z, field=self)
+        if y < grid.lonlat_minmax[2] or y > grid.lonlat_minmax[3]:
             raise FieldSamplingError(x, y, z, field=self)
 
         while xsi < 0 or xsi > 1 or eta < 0 or eta > 1:

parcels/grid.py

@@ -52,7 +52,7 @@ def __init__(self, lon, lat, time, time_origin, mesh):
         self.meridional_halo = 0
         self.lat_flipped = False
         self.defer_load = False
-        self.lonlat_minmax = np.array([np.min(lon), np.max(lon), np.min(lat), np.max(lat)])
+        self.lonlat_minmax = np.array([np.min(lon), np.max(lon), np.min(lat), np.max(lat)], dtype=np.float32)
 
     @property
     def ctypes_struct(self):
@@ -203,6 +203,7 @@ def add_periodic_halo(self, zonal, meridional, halosize=5):
                                       self.lat, self.lat[0:halosize] + latshift))
             self.ydim = self.lat.size
             self.meridional_halo = halosize
+        self.lonlat_minmax = np.array([np.min(self.lon), np.max(self.lon), np.min(self.lat), np.max(self.lat)], dtype=np.float32)
 
 
 class RectilinearZGrid(RectilinearGrid):

parcels/include/index_search.h

@@ -129,21 +129,24 @@ static inline ErrorCode search_indices_rectilinear(float x, float y, float z, CS
   float *xvals = grid->lon;
   float *yvals = grid->lat;
   float *zvals = grid->depth;
+  float *xy_minmax = grid->lonlat_minmax;
   int sphere_mesh = grid->sphere_mesh;
   int zonal_periodic = grid->zonal_periodic;
   int z4d = grid->z4d;
 
+  if (zonal_periodic == 0){
+    if ((x < xy_minmax[0]) || (x > xy_minmax[1]))
+      return ERROR_OUT_OF_BOUNDS;
+  }
+  if ((y < xy_minmax[2]) || (y > xy_minmax[3]))
+    return ERROR_OUT_OF_BOUNDS;
+
   if (sphere_mesh == 0){
-    if (x < xvals[0] || x > xvals[xdim-1]) {return ERROR_OUT_OF_BOUNDS;}
     while (*xi < xdim-1 && x > xvals[*xi+1]) ++(*xi);
     while (*xi > 0 && x < xvals[*xi]) --(*xi);
     *xsi = (x - xvals[*xi]) / (xvals[*xi+1] - xvals[*xi]);
   }
   else{
-    if (zonal_periodic == 0){
-      if ((xvals[0] < xvals[xdim-1]) && (x < xvals[0] || x > xvals[xdim-1])) {return ERROR_OUT_OF_BOUNDS;}
-      else if ((xvals[0] >= xvals[xdim-1]) && (x < xvals[0] && x > xvals[xdim-1])) {return ERROR_OUT_OF_BOUNDS;}
-    }
 
     float xvalsi = xvals[*xi];
     // TODO: this will fail if longitude is e.g. only [-180, 180] (so length 2)
@@ -176,7 +179,6 @@ static inline ErrorCode search_indices_rectilinear(float x, float y, float z, CS
     *xsi = (x - xvalsi) / (xvalsi1 - xvalsi);
   }
 
-  if (y < yvals[0] || y > yvals[ydim-1]) {return ERROR_OUT_OF_BOUNDS;}
   while (*yi < ydim-1 && y > yvals[*yi+1]) ++(*yi);
   while (*yi > 0 && y < yvals[*yi]) --(*yi);
 
@@ -220,6 +222,7 @@ static inline ErrorCode search_indices_curvilinear(float x, float y, float z, CS
   float *xvals = grid->lon;
   float *yvals = grid->lat;
   float *zvals = grid->depth;
+  float *xy_minmax = grid->lonlat_minmax;
   int sphere_mesh = grid->sphere_mesh;
   int zonal_periodic = grid->zonal_periodic;
   int z4d = grid->z4d;
@@ -228,10 +231,14 @@ static inline ErrorCode search_indices_curvilinear(float x, float y, float z, CS
   float (* xgrid)[xdim] = (float (*)[xdim]) xvals;
   float (* ygrid)[xdim] = (float (*)[xdim]) yvals;
 
-  if (zonal_periodic == 0 || sphere_mesh == 0) {
-    if ((xgrid[0][0] < xgrid[0][xdim-1]) && (x < xgrid[0][0] || x > xgrid[0][xdim-1])) {return ERROR_OUT_OF_BOUNDS;}
-    else if ((xgrid[0][0] >= xgrid[0][xdim-1]) && (x < xgrid[0][0] && x > xgrid[0][xdim-1])) {return ERROR_OUT_OF_BOUNDS;}
+  if (zonal_periodic == 0){
+    if ((x < xy_minmax[0]) || (x > xy_minmax[1])){
+      if (xgrid[0][0] < xgrid[0][xdim-1]) {return ERROR_OUT_OF_BOUNDS;}
+      else if (x < xgrid[0][0] && x > xgrid[0][xdim-1]) {return ERROR_OUT_OF_BOUNDS;}
+    }
   }
+  if ((y < xy_minmax[2]) || (y > xy_minmax[3]))
+    return ERROR_OUT_OF_BOUNDS;
 
   double a[4], b[4];