Skip to content

Reordering dimensions to [time, depth, lat, lon] in Grid class #1813

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
erikvansebille wants to merge 16 commits into from
Closed

Reordering dimensions to [time, depth, lat, lon] in Grid class #1813

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
0118b65
Updating variable order in parcels.h interpolation
erikvansebille Jan 2, 2025
8a1daf0
Updating variable order in spatial_interpolation methods
erikvansebille Jan 2, 2025
8d4e8c5
Renaming zet to zeta throughout
erikvansebille Jan 2, 2025
26c467f
Remvoing non-used phi3D_lin function
erikvansebille Jan 2, 2025
998a047
Reordering variable order in phi2D_in and associated functions
erikvansebille Jan 2, 2025
c72e7e7
Reordering variables in search_indices and getCell methods
erikvansebille Jan 2, 2025
4b697be
Reordering variables in field.py search_indices
erikvansebille Jan 2, 2025
371fda9
Updating unitconverters to z,y,x order
erikvansebille Jan 2, 2025
a7ea578
Updating Error codes to use depth,lat,lon order
erikvansebille Jan 2, 2025
09533f1
Removing unused Interpolation Error
erikvansebille Jan 2, 2025
54ee8e8
Update variable order in reconnect_bnd_indices
erikvansebille Jan 2, 2025
eb01ba4
Update variable order in dist() function
erikvansebille Jan 2, 2025
15a54cf
Updating variable order in jacobian calculations
erikvansebille Jan 2, 2025
9ad4246
Removing unused length2d_lin_edge() method
erikvansebille Jan 2, 2025
e3b00cd
Reordering dimensions to time,depth,lat,lon in Grid class
erikvansebille Jan 3, 2025
0a7ab33
fixing mypy typing errors
erikvansebille Jan 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 7 additions & 7 deletions parcels/application_kernels/advection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,8 +184,8 @@ def AdvectionAnalytical(particle, fieldset, time):
time_i = np.linspace(0, fieldset.U.grid.time[ti + 1] - fieldset.U.grid.time[ti], I_s)
ds_t = min(ds_t, time_i[np.where(time - fieldset.U.grid.time[ti] < time_i)[0][0]])

xsi, eta, zeta, xi, yi, zi = fieldset.U._search_indices(
particle.lon, particle.lat, particle.depth, particle=particle
zeta, eta, xsi, zi, yi, xi = fieldset.U._search_indices(
-1, particle.depth, particle.lat, particle.lon, particle=particle
)
if withW:
if abs(xsi - 1) < tol:
Expand Down Expand Up @@ -232,14 +232,14 @@ def AdvectionAnalytical(particle, fieldset, time):
else:
dz = 1.0

c1 = fieldset.UV.dist(px[0], px[1], py[0], py[1], grid.mesh, np.dot(i_u.phi2D_lin(xsi, 0.0), py))
c2 = fieldset.UV.dist(px[1], px[2], py[1], py[2], grid.mesh, np.dot(i_u.phi2D_lin(1.0, eta), py))
c3 = fieldset.UV.dist(px[2], px[3], py[2], py[3], grid.mesh, np.dot(i_u.phi2D_lin(xsi, 1.0), py))
c4 = fieldset.UV.dist(px[3], px[0], py[3], py[0], grid.mesh, np.dot(i_u.phi2D_lin(0.0, eta), py))
c1 = fieldset.UV.dist(py[0], py[1], px[0], px[1], grid.mesh, np.dot(i_u.phi2D_lin(0.0, xsi), py))
c2 = fieldset.UV.dist(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
c3 = fieldset.UV.dist(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
c4 = fieldset.UV.dist(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
rad = np.pi / 180.0
deg2m = 1852 * 60.0
meshJac = (deg2m * deg2m * math.cos(rad * particle.lat)) if grid.mesh == "spherical" else 1
dxdy = fieldset.UV.jacobian(xsi, eta, px, py) * meshJac
dxdy = fieldset.UV.jacobian(py, px, eta, xsi) * meshJac

if withW:
U0 = direction * fieldset.U.data[ti, zi + 1, yi + 1, xi] * c4 * dz
Expand Down
4 changes: 2 additions & 2 deletions parcels/compilation/codegenerator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -834,7 +834,7 @@ def visit_FieldEvalNode(self, node):
statements_croco = [
c.Statement(f"float cs_w[] = {*Cs_w, }".replace("(", "{").replace(")", "}")),
c.Statement(
f"{node.var} = croco_from_z_to_sigma(U, H, Zeta, {args[3]}, {args[2]}, {args[1]}, time, &particles->xi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->ti[pnum*ngrid], hc, &cs_w)"
f"{node.var} = croco_from_z_to_sigma(time, {args[1]}, {args[2]}, {args[3]}, U, H, Zeta, &particles->ti[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->xi[pnum*ngrid], hc, &cs_w)"
),
]
args = (args[0], node.var, args[2], args[3])
Expand Down Expand Up @@ -863,7 +863,7 @@ def visit_VectorFieldEvalNode(self, node):
statements_croco = [
c.Statement(f"float cs_w[] = {*Cs_w, }".replace("(", "{").replace(")", "}")),
c.Statement(
f"{node.var4} = croco_from_z_to_sigma(U, H, Zeta, {args[3]}, {args[2]}, {args[1]}, time, &particles->xi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->ti[pnum*ngrid], hc, &cs_w)"
f"{node.var4} = croco_from_z_to_sigma(time, {args[1]}, {args[2]}, {args[3]}, U, H, Zeta, &particles->ti[pnum*ngrid], &particles->zi[pnum*ngrid], &particles->yi[pnum*ngrid], &particles->xi[pnum*ngrid], hc, &cs_w)"
),
]
args = (args[0], node.var4, args[2], args[3])
Expand Down
194 changes: 100 additions & 94 deletions parcels/field.py
View file
Open in desktop

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion parcels/fieldset.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ def from_data(
time = np.array([time_origin.reltime(t) for t in time])
else:
time_origin = kwargs.pop("time_origin", TimeConverter(0))
grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
grid = Grid.create_grid(time, depth, lat, lon, time_origin=time_origin, mesh=mesh)
if "creation_log" not in kwargs.keys():
kwargs["creation_log"] = "from_data"

Expand Down
67 changes: 35 additions & 32 deletions parcels/grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,9 +44,10 @@

def __init__(
self,
lon: npt.NDArray,
lat: npt.NDArray,
time: npt.NDArray | None,
depth: npt.NDArray | None,
lat: npt.NDArray,
lon: npt.NDArray,
time_origin: TimeConverter | None,
mesh: Mesh,
):
Expand All @@ -71,6 +72,7 @@

self._lon = lon
self._lat = lat
self._depth = depth
self.time = time
self.time_full = self.time # needed for deferred_loaded Fields
self._time_origin = TimeConverter() if time_origin is None else time_origin
Expand Down Expand Up @@ -98,7 +100,7 @@
with np.printoptions(threshold=5, suppress=True, linewidth=120, formatter={"float": "{: 0.2f}".format}):
return (
f"{type(self).__name__}("
f"lon={self.lon!r}, lat={self.lat!r}, time={self.time!r}, "
f"lon={self.lon!r}, lat={self.lat!r}, time={self.time!r}, depth={self.depth!r}, "

Check warning on line 103 in parcels/grid.py

View check run for this annotation

Codecov / codecov/patch

parcels/grid.py#L103 

Added line #L103 was not covered by tests
f"time_origin={self.time_origin!r}, mesh={self.mesh!r})"
)

Expand Down Expand Up @@ -195,10 +197,10 @@

@staticmethod
def create_grid(
lon: npt.ArrayLike,
lat: npt.ArrayLike,
depth,
time,
time: npt.NDArray | None,
depth: npt.NDArray | None,
lat: npt.NDArray,
lon: npt.NDArray,
time_origin,
mesh: Mesh,
**kwargs,
Expand All @@ -211,14 +213,14 @@

if len(lon.shape) <= 1:
if depth is None or len(depth.shape) <= 1:
return RectilinearZGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh, **kwargs)
return RectilinearZGrid(time, depth, lat, lon, time_origin=time_origin, mesh=mesh, **kwargs)
else:
return RectilinearSGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh, **kwargs)
return RectilinearSGrid(time, depth, lat, lon, time_origin=time_origin, mesh=mesh, **kwargs)
else:
if depth is None or len(depth.shape) <= 1:
return CurvilinearZGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh, **kwargs)
return CurvilinearZGrid(time, depth, lat, lon, time_origin=time_origin, mesh=mesh, **kwargs)
else:
return CurvilinearSGrid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh, **kwargs)
return CurvilinearSGrid(time, depth, lat, lon, time_origin=time_origin, mesh=mesh, **kwargs)

@property
def ctypes_struct(self):
Expand Down Expand Up @@ -461,14 +463,14 @@

"""

def __init__(self, lon, lat, time, time_origin, mesh: Mesh):
def __init__(self, time, depth, lat, lon, time_origin, mesh: Mesh):
assert isinstance(lon, np.ndarray) and len(lon.shape) <= 1, "lon is not a numpy vector"
assert isinstance(lat, np.ndarray) and len(lat.shape) <= 1, "lat is not a numpy vector"
assert isinstance(time, np.ndarray) or not time, "time is not a numpy array"
if isinstance(time, np.ndarray):
assert len(time.shape) == 1, "time is not a vector"

super().__init__(lon, lat, time, time_origin, mesh)
super().__init__(time, depth, lat, lon, time_origin, mesh)
self.tdim = self.time.size

if self.ydim > 1 and self.lat[-1] < self.lat[0]:
Expand Down Expand Up @@ -559,8 +561,8 @@
2. flat: No conversion, lat/lon are assumed to be in m.
"""

def __init__(self, lon, lat, depth=None, time=None, time_origin=None, mesh: Mesh = "flat"):
super().__init__(lon, lat, time, time_origin, mesh)
def __init__(self, time, depth, lat, lon, time_origin=None, mesh: Mesh = "flat"):
super().__init__(time, depth, lat, lon, time_origin, mesh)
if isinstance(depth, np.ndarray):
assert len(depth.shape) <= 1, "depth is not a vector"

Expand Down Expand Up @@ -592,7 +594,7 @@
which are s-coordinates.
s-coordinates can be terrain-following (sigma) or iso-density (rho) layers,
or any generalised vertical discretisation.
The depth of each node depends then on the horizontal position (lon, lat),
The depth of each node depends then on the horizontal position (lat, lon),
the number of the layer and the time is depth is a 4D array.
depth array is either a 4D array[xdim][ydim][zdim][tdim] or a 3D array[xdim][ydim[zdim].
time :
Expand All @@ -610,14 +612,14 @@

def __init__(
self,
lon: npt.NDArray,
time: npt.NDArray | None,
depth: npt.NDArray | None,
lat: npt.NDArray,
depth: npt.NDArray,
time: npt.NDArray | None = None,
lon: npt.NDArray,
time_origin: TimeConverter | None = None,
mesh: Mesh = "flat",
):
super().__init__(lon, lat, time, time_origin, mesh)
super().__init__(time, depth, lat, lon, time_origin, mesh)
assert isinstance(depth, np.ndarray) and len(depth.shape) in [3, 4], "depth is not a 3D or 4D numpy array"

self._gtype = GridType.RectilinearSGrid
Expand Down Expand Up @@ -656,9 +658,10 @@
class CurvilinearGrid(Grid):
def __init__(
self,
lon: npt.NDArray,
time: npt.NDArray | None,
depth: npt.NDArray | None,
lat: npt.NDArray,
time: npt.NDArray | None = None,
lon: npt.NDArray,
time_origin: TimeConverter | None = None,
mesh: Mesh = "flat",
):
Expand All @@ -670,7 +673,7 @@

lon = lon.squeeze()
lat = lat.squeeze()
super().__init__(lon, lat, time, time_origin, mesh)
super().__init__(time, None, lat, lon, time_origin, mesh)
self.tdim = self.time.size

@property
Expand Down Expand Up @@ -770,14 +773,14 @@

def __init__(
self,
lon: npt.NDArray,
time: npt.NDArray | None,
depth: npt.NDArray | None,
lat: npt.NDArray,
depth: npt.NDArray | None = None,
time: npt.NDArray | None = None,
lon: npt.NDArray,
time_origin: TimeConverter | None = None,
mesh: Mesh = "flat",
):
super().__init__(lon, lat, time, time_origin, mesh)
super().__init__(time, depth, lat, lon, time_origin, mesh)
if isinstance(depth, np.ndarray):
assert len(depth.shape) == 1, "depth is not a vector"

Expand Down Expand Up @@ -808,7 +811,7 @@
which are s-coordinates.
s-coordinates can be terrain-following (sigma) or iso-density (rho) layers,
or any generalised vertical discretisation.
The depth of each node depends then on the horizontal position (lon, lat),
The depth of each node depends then on the horizontal position (lat, lon),
the number of the layer and the time is depth is a 4D array.
depth array is either a 4D array[xdim][ydim][zdim][tdim] or a 3D array[xdim][ydim[zdim].
time :
Expand All @@ -826,14 +829,14 @@

def __init__(
self,
lon: npt.NDArray,
time: npt.NDArray | None,
depth: npt.NDArray | None,
lat: npt.NDArray,
depth: npt.NDArray,
time: npt.NDArray | None = None,
lon: npt.NDArray,
time_origin: TimeConverter | None = None,
mesh: Mesh = "flat",
):
super().__init__(lon, lat, time, time_origin, mesh)
super().__init__(time, depth, lat, lon, time_origin, mesh)
assert isinstance(depth, np.ndarray) and len(depth.shape) in [3, 4], "depth is not a 4D numpy array"

self._gtype = GridType.CurvilinearSGrid
Expand Down
2 changes: 1 addition & 1 deletion parcels/gridset.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ def add_grid(self, field):
field.igrid = self.grids.index(field.grid)

def dimrange(self, dim):
"""Returns maximum value of a dimension (lon, lat, depth or time)
"""Returns maximum value of a dimension (time, depth, lat or lon)
on 'left' side and minimum value on 'right' side for all grids
in a gridset. Useful for finding e.g. longitude range that
overlaps on all grids in a gridset.
Expand Down
62 changes: 34 additions & 28 deletions parcels/include/index_search.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,9 +127,10 @@ static inline StatusCode search_indices_vertical_z(type_coord z, int zdim, float
return SUCCESS;
}

static inline StatusCode search_indices_vertical_s(type_coord z, int xdim, int ydim, int zdim, float *zvals,
int xi, int yi, int *zi, double xsi, double eta, double *zeta,
int z4d, int ti, int tdim, double time, double t0, double t1, int interp_method)
static inline StatusCode search_indices_vertical_s(double time, type_coord z,
int tdim, int zdim, int ydim, int xdim, float *zvals,
int ti, int *zi, int yi, int xi, double *zeta, double eta, double xsi,
int z4d, double t0, double t1, int interp_method)
{
if (interp_method == BGRID_VELOCITY || interp_method == BGRID_W_VELOCITY || interp_method == BGRID_TRACER){
xsi = 1;
Expand Down Expand Up @@ -184,7 +185,7 @@ static inline StatusCode search_indices_vertical_s(type_coord z, int xdim, int y
return SUCCESS;
}

static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, int onlyX, int sphere_mesh)
static inline void reconnect_bnd_indices(int *yi, int *xi, int ydim, int xdim, int onlyX, int sphere_mesh)
{
if (*xi < 0){
if (sphere_mesh)
Expand All @@ -211,10 +212,11 @@ static inline void reconnect_bnd_indices(int *xi, int *yi, int xdim, int ydim, i
}


static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridType gtype,
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
int ti, double time, double t0, double t1, int interp_method,
int gridindexingtype)
static inline StatusCode search_indices_rectilinear(double time, type_coord z, type_coord y, type_coord x,
CStructuredGrid *grid, GridType gtype,
int ti, int *zi, int *yi, int *xi, double *zeta, double *eta, double *xsi,
double t0, double t1, int interp_method,
int gridindexingtype)
{
int xdim = grid->xdim;
int ydim = grid->ydim;
Expand Down Expand Up @@ -261,7 +263,7 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y,
++(*xi);
else if (xvalsi > x)
--(*xi);
reconnect_bnd_indices(xi, yi, xdim, ydim, 1, 1);
reconnect_bnd_indices(yi, xi, ydim, xdim, 1, 1);
xvalsi = xvals[*xi];
if (xvalsi < x - 225) xvalsi += 360;
if (xvalsi > x + 225) xvalsi -= 360;
Expand Down Expand Up @@ -294,9 +296,9 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y,
status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype);
break;
case RECTILINEAR_S_GRID:
status = search_indices_vertical_s(z, xdim, ydim, zdim, zvals,
*xi, *yi, zi, *xsi, *eta, zeta,
z4d, ti, tdim, time, t0, t1, interp_method);
status = search_indices_vertical_s(time, z, tdim, zdim, ydim, xdim, zvals,
ti, zi, *yi, *xi, zeta, *eta, *xsi,
z4d, t0, t1, interp_method);
break;
default:
status = ERRORINTERPOLATION;
Expand All @@ -321,10 +323,12 @@ static inline StatusCode search_indices_rectilinear(type_coord x, type_coord y,
}


static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid, GridType gtype,
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
int ti, double time, double t0, double t1, int interp_method,
int gridindexingtype)
static inline StatusCode search_indices_curvilinear(double time, type_coord z, type_coord y, type_coord x,
CStructuredGrid *grid, GridType gtype,
int ti, int *zi, int *yi, int *xi,
double *zeta, double *eta, double *xsi,
double t0, double t1, int interp_method,
int gridindexingtype)
{
int xi_old = *xi;
int yi_old = *yi;
Expand Down Expand Up @@ -407,7 +411,7 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y,
(*yi)--;
if (*eta > 1+tol)
(*yi)++;
reconnect_bnd_indices(xi, yi, xdim, ydim, 0, sphere_mesh);
reconnect_bnd_indices(yi, xi, ydim, xdim, 0, sphere_mesh);
it++;
if ( it > maxIterSearch){
printf("Correct cell not found for (%f, %f) after %d iterations\n", x, y, maxIterSearch);
Expand Down Expand Up @@ -438,9 +442,9 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y,
status = search_indices_vertical_z(z, zdim, zvals, zi, zeta, gridindexingtype);
break;
case CURVILINEAR_S_GRID:
status = search_indices_vertical_s(z, xdim, ydim, zdim, zvals,
*xi, *yi, zi, *xsi, *eta, zeta,
z4d, ti, tdim, time, t0, t1, interp_method);
status = search_indices_vertical_s(time, z, tdim, ydim, xdim, zdim, zvals,
ti, zi, *yi, *xi, zeta, *eta, *xsi,
z4d, t0, t1, interp_method);
break;
default:
status = ERRORINTERPOLATION;
Expand All @@ -460,21 +464,23 @@ static inline StatusCode search_indices_curvilinear(type_coord x, type_coord y,
/* Local linear search to update grid index
* params ti, sizeT, time. t0, t1 are only used for 4D S grids
* */
static inline StatusCode search_indices(type_coord x, type_coord y, type_coord z, CStructuredGrid *grid,
int *xi, int *yi, int *zi, double *xsi, double *eta, double *zeta,
GridType gtype, int ti, double time, double t0, double t1, int interp_method,
int gridindexingtype)
static inline StatusCode search_indices(double time, type_coord z, type_coord y, type_coord x,
CStructuredGrid *grid,
int ti, int *zi, int *yi, int *xi,
double *zeta, double *eta, double *xsi,
GridType gtype, double t0, double t1, int interp_method,
int gridindexingtype)
{
switch(gtype){
case RECTILINEAR_Z_GRID:
case RECTILINEAR_S_GRID:
return search_indices_rectilinear(x, y, z, grid, gtype, xi, yi, zi, xsi, eta, zeta,
ti, time, t0, t1, interp_method, gridindexingtype);
return search_indices_rectilinear(time, z, y, x, grid, gtype, ti, zi, yi, xi, zeta, eta, xsi,
t0, t1, interp_method, gridindexingtype);
break;
case CURVILINEAR_Z_GRID:
case CURVILINEAR_S_GRID:
return search_indices_curvilinear(x, y, z, grid, gtype, xi, yi, zi, xsi, eta, zeta,
ti, time, t0, t1, interp_method, gridindexingtype);
return search_indices_curvilinear(time, z, y, x, grid, gtype, ti, zi, yi, xi, zeta, eta, xsi,
t0, t1, interp_method, gridindexingtype);
break;
default:
printf("Only RECTILINEAR_Z_GRID, RECTILINEAR_S_GRID, CURVILINEAR_Z_GRID and CURVILINEAR_S_GRID grids are currently implemented\n");
Expand Down
Loading
Loading