Skip to content

Commit 0e6ade3

Browse files
delandmeterpdelandmeterp
authored
Merge pull request #440 from OceanParcels/lonlat_minmax_values
Lonlat minmax values
2 parents e006cb2 + f14e3ba commit 0e6ade3

5 files changed

Lines changed: 89 additions & 61 deletions

File tree

‎parcels/field.py‎

Lines changed: 13 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -403,9 +403,13 @@ def search_indices_rectilinear(self, x, y, z, ti=-1, time=-1, search2D=False):
403403
grid = self.grid
404404
xi = yi = -1
405405

406-
if grid.mesh is not 'spherical':
407-
if x < grid.lon[0] or x > grid.lon[-1]:
406+
if not grid.zonal_periodic:
407+
if x < grid.lonlat_minmax[0] or x > grid.lonlat_minmax[1]:
408408
raise FieldSamplingError(x, y, z, field=self)
409+
if y < grid.lonlat_minmax[2] or y > grid.lonlat_minmax[3]:
410+
raise FieldSamplingError(x, y, z, field=self)
411+
412+
if grid.mesh is not 'spherical':
409413
lon_index = grid.lon < x
410414
if lon_index.all():
411415
xi = len(grid.lon) - 2
@@ -425,11 +429,6 @@ def search_indices_rectilinear(self, x, y, z, ti=-1, time=-1, search2D=False):
425429
lon_fixed[indices.argmin():] += 360
426430
if x < lon_fixed[0]:
427431
lon_fixed -= 360
428-
if not grid.zonal_periodic:
429-
if (grid.lon[0] < grid.lon[-1]) and (x < grid.lon[0] or x > grid.lon[-1]):
430-
raise FieldSamplingError(x, y, z, field=self)
431-
elif (grid.lon[0] >= grid.lon[-1]) and (x < grid.lon[0] and x > grid.lon[-1]):
432-
raise FieldSamplingError(x, y, z, field=self)
433432

434433
lon_index = lon_fixed < x
435434
if lon_index.all():
@@ -444,8 +443,6 @@ def search_indices_rectilinear(self, x, y, z, ti=-1, time=-1, search2D=False):
444443
xi += 1
445444
xsi = (x-lon_fixed[xi]) / (lon_fixed[xi+1]-lon_fixed[xi])
446445

447-
if y < grid.lat[0] or y > grid.lat[-1]:
448-
raise FieldSamplingError(x, y, z, field=self)
449446
lat_index = grid.lat < y
450447
if lat_index.all():
451448
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
485482
[1, -1, 1, -1]])
486483
maxIterSearch = 1e6
487484
it = 0
488-
if (not grid.zonal_periodic) or grid.mesh == 'flat':
489-
if (grid.lon[0, 0] < grid.lon[0, -1]) and (x < grid.lon[0, 0] or x > grid.lon[0, -1]):
490-
raise FieldSamplingError(x, y, z, field=self)
491-
elif (grid.lon[0, 0] >= grid.lon[0, -1]) and (x < grid.lon[0, 0] and x > grid.lon[0, -1]):
492-
raise FieldSamplingError(x, y, z, field=self)
493-
if y < np.min(grid.lat) or y > np.max(grid.lat):
485+
if not grid.zonal_periodic:
486+
if x < grid.lonlat_minmax[0] or x > grid.lonlat_minmax[1]:
487+
if grid.lon[0, 0] < grid.lon[0, -1]:
488+
raise FieldSamplingError(x, y, z, field=self)
489+
elif x < grid.lon[0, 0] and x > grid.lon[0, -1]: # This prevents from crashing in [160, -160]
490+
raise FieldSamplingError(x, y, z, field=self)
491+
if y < grid.lonlat_minmax[2] or y > grid.lonlat_minmax[3]:
494492
raise FieldSamplingError(x, y, z, field=self)
495493

496494
while xsi < 0 or xsi > 1 or eta < 0 or eta > 1:

‎parcels/grid.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,7 @@ def __init__(self, lon, lat, time, time_origin, mesh):
5252
self.meridional_halo = 0
5353
self.lat_flipped = False
5454
self.defer_load = False
55+
self.lonlat_minmax = np.array([np.nanmin(lon), np.nanmax(lon), np.nanmin(lat), np.nanmax(lat)], dtype=np.float32)
5556

5657
@property
5758
def ctypes_struct(self):
@@ -70,6 +71,7 @@ class CStructuredGrid(Structure):
7071
_fields_ = [('xdim', c_int), ('ydim', c_int), ('zdim', c_int),
7172
('tdim', c_int), ('z4d', c_int),
7273
('mesh_spherical', c_int), ('zonal_periodic', c_int),
74+
('lonlat_minmax', POINTER(c_float)),
7375
('lon', POINTER(c_float)), ('lat', POINTER(c_float)),
7476
('depth', POINTER(c_float)), ('time', POINTER(c_double))
7577
]
@@ -79,6 +81,7 @@ class CStructuredGrid(Structure):
7981
self.cstruct = CStructuredGrid(self.xdim, self.ydim, self.zdim,
8082
self.tdim, self.z4d,
8183
self.mesh == 'spherical', self.zonal_periodic,
84+
self.lonlat_minmax.ctypes.data_as(POINTER(c_float)),
8285
self.lon.ctypes.data_as(POINTER(c_float)),
8386
self.lat.ctypes.data_as(POINTER(c_float)),
8487
self.depth.ctypes.data_as(POINTER(c_float)),
@@ -200,6 +203,7 @@ def add_periodic_halo(self, zonal, meridional, halosize=5):
200203
self.lat, self.lat[0:halosize] + latshift))
201204
self.ydim = self.lat.size
202205
self.meridional_halo = halosize
206+
self.lonlat_minmax = np.array([np.nanmin(self.lon), np.nanmax(self.lon), np.nanmin(self.lat), np.nanmax(self.lat)], dtype=np.float32)
203207

204208

205209
class RectilinearZGrid(RectilinearGrid):

‎parcels/include/index_search.h‎

Lines changed: 67 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,22 @@ extern "C" {
1010

1111
#define CHECKERROR(res) do {if (res != SUCCESS) return res;} while (0)
1212

13+
typedef struct
14+
{
15+
int gtype;
16+
void *grid;
17+
} CGrid;
18+
19+
typedef struct
20+
{
21+
int xdim, ydim, zdim, tdim, z4d;
22+
int sphere_mesh, zonal_periodic;
23+
float *lonlat_minmax;
24+
float *lon, *lat, *depth;
25+
double *time;
26+
} CStructuredGrid;
27+
28+
1329
typedef enum
1430
{
1531
SUCCESS=0, REPEAT=1, DELETE=2, ERROR=3, ERROR_OUT_OF_BOUNDS=4, ERROR_TIME_EXTRAPOLATION =5
@@ -102,22 +118,35 @@ static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, i
102118
}
103119

104120

105-
static inline ErrorCode search_indices_rectilinear(float x, float y, float z, int xdim, int ydim, int zdim,
106-
float *xvals, float *yvals, float *zvals, int sphere_mesh, int zonal_periodic, GridCode gcode,
107-
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
108-
int z4d, int ti, int tdim, double time, double t0, double t1)
121+
static inline ErrorCode search_indices_rectilinear(float x, float y, float z, CStructuredGrid *grid, GridCode gcode,
122+
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
123+
int ti, double time, double t0, double t1)
109124
{
125+
int xdim = grid->xdim;
126+
int ydim = grid->ydim;
127+
int zdim = grid->zdim;
128+
int tdim = grid->tdim;
129+
float *xvals = grid->lon;
130+
float *yvals = grid->lat;
131+
float *zvals = grid->depth;
132+
float *xy_minmax = grid->lonlat_minmax;
133+
int sphere_mesh = grid->sphere_mesh;
134+
int zonal_periodic = grid->zonal_periodic;
135+
int z4d = grid->z4d;
136+
137+
if (zonal_periodic == 0){
138+
if ((x < xy_minmax[0]) || (x > xy_minmax[1]))
139+
return ERROR_OUT_OF_BOUNDS;
140+
}
141+
if ((y < xy_minmax[2]) || (y > xy_minmax[3]))
142+
return ERROR_OUT_OF_BOUNDS;
143+
110144
if (sphere_mesh == 0){
111-
if (x < xvals[0] || x > xvals[xdim-1]) {return ERROR_OUT_OF_BOUNDS;}
112145
while (*xi < xdim-1 && x > xvals[*xi+1]) ++(*xi);
113146
while (*xi > 0 && x < xvals[*xi]) --(*xi);
114147
*xsi = (x - xvals[*xi]) / (xvals[*xi+1] - xvals[*xi]);
115148
}
116149
else{
117-
if (zonal_periodic == 0){
118-
if ((xvals[0] < xvals[xdim-1]) && (x < xvals[0] || x > xvals[xdim-1])) {return ERROR_OUT_OF_BOUNDS;}
119-
else if ((xvals[0] >= xvals[xdim-1]) && (x < xvals[0] && x > xvals[xdim-1])) {return ERROR_OUT_OF_BOUNDS;}
120-
}
121150

122151
float xvalsi = xvals[*xi];
123152
// TODO: this will fail if longitude is e.g. only [-180, 180] (so length 2)
@@ -150,7 +179,6 @@ static inline ErrorCode search_indices_rectilinear(float x, float y, float z, in
150179
*xsi = (x - xvalsi) / (xvalsi1 - xvalsi);
151180
}
152181

153-
if (y < yvals[0] || y > yvals[ydim-1]) {return ERROR_OUT_OF_BOUNDS;}
154182
while (*yi < ydim-1 && y > yvals[*yi+1]) ++(*yi);
155183
while (*yi > 0 && y < yvals[*yi]) --(*yi);
156184

@@ -183,19 +211,34 @@ static inline ErrorCode search_indices_rectilinear(float x, float y, float z, in
183211
}
184212

185213

186-
static inline ErrorCode search_indices_curvilinear(float x, float y, float z, int xdim, int ydim, int zdim,
187-
float *xvals, float *yvals, float *zvals, int sphere_mesh, int zonal_periodic, GridCode gcode,
188-
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
189-
int z4d, int ti, int tdim, double time, double t0, double t1)
214+
static inline ErrorCode search_indices_curvilinear(float x, float y, float z, CStructuredGrid *grid, GridCode gcode,
215+
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
216+
int ti, double time, double t0, double t1)
190217
{
218+
int xdim = grid->xdim;
219+
int ydim = grid->ydim;
220+
int zdim = grid->zdim;
221+
int tdim = grid->tdim;
222+
float *xvals = grid->lon;
223+
float *yvals = grid->lat;
224+
float *zvals = grid->depth;
225+
float *xy_minmax = grid->lonlat_minmax;
226+
int sphere_mesh = grid->sphere_mesh;
227+
int zonal_periodic = grid->zonal_periodic;
228+
int z4d = grid->z4d;
229+
191230
// NEMO convention
192231
float (* xgrid)[xdim] = (float (*)[xdim]) xvals;
193232
float (* ygrid)[xdim] = (float (*)[xdim]) yvals;
194233

195-
if (zonal_periodic == 0 || sphere_mesh == 0) {
196-
if ((xgrid[0][0] < xgrid[0][xdim-1]) && (x < xgrid[0][0] || x > xgrid[0][xdim-1])) {return ERROR_OUT_OF_BOUNDS;}
197-
else if ((xgrid[0][0] >= xgrid[0][xdim-1]) && (x < xgrid[0][0] && x > xgrid[0][xdim-1])) {return ERROR_OUT_OF_BOUNDS;}
234+
if (zonal_periodic == 0){
235+
if ((x < xy_minmax[0]) || (x > xy_minmax[1])){
236+
if (xgrid[0][0] < xgrid[0][xdim-1]) {return ERROR_OUT_OF_BOUNDS;}
237+
else if (x < xgrid[0][0] && x > xgrid[0][xdim-1]) {return ERROR_OUT_OF_BOUNDS;}
238+
}
198239
}
240+
if ((y < xy_minmax[2]) || (y > xy_minmax[3]))
241+
return ERROR_OUT_OF_BOUNDS;
199242

200243
double a[4], b[4];
201244

@@ -292,23 +335,20 @@ static inline ErrorCode search_indices_curvilinear(float x, float y, float z, in
292335
/* Local linear search to update grid index
293336
* params ti, sizeT, time. t0, t1 are only used for 4D S grids
294337
* */
295-
static inline ErrorCode search_indices(float x, float y, float z, int xdim, int ydim, int zdim,
296-
float *xvals, float *yvals, float *zvals,
297-
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
298-
int sphere_mesh, int zonal_periodic,
299-
GridCode gcode, int z4d,
300-
int ti, int tdim, double time, double t0, double t1)
338+
static inline ErrorCode search_indices(float x, float y, float z, CStructuredGrid *grid,
339+
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
340+
GridCode gcode, int ti, double time, double t0, double t1)
301341
{
302342
switch(gcode){
303343
case RECTILINEAR_Z_GRID:
304344
case RECTILINEAR_S_GRID:
305-
return search_indices_rectilinear(x, y, z, xdim, ydim, zdim, xvals, yvals, zvals, sphere_mesh, zonal_periodic, gcode, xi, yi, zi, xsi, eta, zeta,
306-
z4d, ti, tdim, time, t0, t1);
345+
return search_indices_rectilinear(x, y, z, grid, gcode, xi, yi, zi, xsi, eta, zeta,
346+
ti, time, t0, t1);
307347
break;
308348
case CURVILINEAR_Z_GRID:
309349
case CURVILINEAR_S_GRID:
310-
return search_indices_curvilinear(x, y, z, xdim, ydim, zdim, xvals, yvals, zvals, sphere_mesh, zonal_periodic, gcode, xi, yi, zi, xsi, eta, zeta,
311-
z4d, ti, tdim, time, t0, t1);
350+
return search_indices_curvilinear(x, y, z, grid, gcode, xi, yi, zi, xsi, eta, zeta,
351+
ti, time, t0, t1);
312352
break;
313353
default:
314354
printf("Only RECTILINEAR_Z_GRID, RECTILINEAR_S_GRID, CURVILINEAR_Z_GRID and CURVILINEAR_S_GRID grids are currently implemented\n");

‎parcels/include/parcels.h‎

Lines changed: 4 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -18,20 +18,6 @@ typedef enum
1818
LINEAR=0, NEAREST=1, CGRID_LINEAR=2
1919
} InterpCode;
2020

21-
typedef struct
22-
{
23-
int gtype;
24-
void *grid;
25-
} CGrid;
26-
27-
typedef struct
28-
{
29-
int xdim, ydim, zdim, tdim, z4d;
30-
int sphere_mesh, zonal_periodic;
31-
float *lon, *lat, *depth;
32-
double *time;
33-
} CStructuredGrid;
34-
3521
typedef struct
3622
{
3723
int xdim, ydim, zdim, tdim, igrid, allow_time_extrapolation, time_periodic;
@@ -119,7 +105,7 @@ static inline ErrorCode temporal_interpolation_structured_grid(float x, float y,
119105
float f0, f1;
120106
double t0 = grid->time[ti[igrid]]; double t1 = grid->time[ti[igrid]+1];
121107
/* Identify grid cell to sample through local linear search */
122-
err = search_indices(x, y, z, grid->xdim, grid->ydim, grid->zdim, grid->lon, grid->lat, grid->depth, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, grid->sphere_mesh, grid->zonal_periodic, gcode, grid->z4d, ti[igrid], grid->tdim, time, t0, t1); CHECKERROR(err);
108+
err = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gcode, ti[igrid], time, t0, t1); CHECKERROR(err);
123109
if (interp_method == CGRID_LINEAR){
124110
xsi = 0.;
125111
eta = 0.;
@@ -152,7 +138,7 @@ static inline ErrorCode temporal_interpolation_structured_grid(float x, float y,
152138
return SUCCESS;
153139
} else {
154140
double t0 = grid->time[ti[igrid]];
155-
err = search_indices(x, y, z, grid->xdim, grid->ydim, grid->zdim, grid->lon, grid->lat, grid->depth, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, grid->sphere_mesh, grid->zonal_periodic, gcode, grid->z4d, ti[igrid], grid->tdim, t0, t0, t0+1); CHECKERROR(err);
141+
err = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gcode, ti[igrid], t0, t0, t0+1); CHECKERROR(err);
156142
if (interp_method == CGRID_LINEAR){
157143
xsi = 0.;
158144
eta = 0.;
@@ -301,7 +287,7 @@ static inline ErrorCode temporal_interpolation_UV_c_grid(float x, float y, float
301287
float u0, u1, v0, v1;
302288
double t0 = grid->time[ti[igrid]]; double t1 = grid->time[ti[igrid]+1];
303289
/* Identify grid cell to sample through local linear search */
304-
err = search_indices(x, y, z, grid->xdim, grid->ydim, grid->zdim, grid->lon, grid->lat, grid->depth, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, grid->sphere_mesh, grid->zonal_periodic, gcode, grid->z4d, ti[igrid], grid->tdim, time, t0, t1); CHECKERROR(err);
290+
err = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gcode, ti[igrid], time, t0, t1); CHECKERROR(err);
305291
if (grid->zdim==1){
306292
err = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gcode, (float**)(dataU[ti[igrid]]) , (float**)(dataV[ti[igrid]]), &u0, &v0);
307293
err = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gcode, (float**)(dataU[ti[igrid]+1]), (float**)(dataV[ti[igrid]+1]), &u1, &v1);
@@ -314,7 +300,7 @@ static inline ErrorCode temporal_interpolation_UV_c_grid(float x, float y, float
314300
return SUCCESS;
315301
} else {
316302
double t0 = grid->time[ti[igrid]];
317-
err = search_indices(x, y, z, grid->xdim, grid->ydim, grid->zdim, grid->lon, grid->lat, grid->depth, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, grid->sphere_mesh, grid->zonal_periodic, gcode, grid->z4d, ti[igrid], grid->tdim, t0, t0, t0+1); CHECKERROR(err);
303+
err = search_indices(x, y, z, grid, &xi[igrid], &yi[igrid], &zi[igrid], &xsi, &eta, &zeta, gcode, ti[igrid], t0, t0, t0+1); CHECKERROR(err);
318304
if (grid->zdim==1){
319305
err = spatial_interpolation_UV_c_grid(xsi, eta, xi[igrid], yi[igrid], grid, gcode, (float**)(dataU[ti[igrid]][zi[igrid]]) , (float**)(dataV[ti[igrid]]), u, v);
320306
}

‎tests/test_grids.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -327,7 +327,7 @@ def sampleSpeed(particle, fieldset, time, dt):
327327
class MyParticle(ptype[mode]):
328328
speed = Variable('speed', dtype=np.float32, initial=0.)
329329

330-
pset = ParticleSet.from_list(field_set, MyParticle, lon=[400], lat=[600])
330+
pset = ParticleSet.from_list(field_set, MyParticle, lon=[400, -200], lat=[600, 600])
331331
pset.execute(pset.Kernel(sampleSpeed), runtime=0, dt=0)
332332
assert(np.allclose(pset[0].speed, 1000))
333333

0 commit comments

Comments
 (0)