Vertical Coordinate Interpolators by wyatt-fluidnumerics · Pull Request #3037 · Parcels-code/Parcels

wyatt-fluidnumerics · 2026-11-04T22:07:48Z

Description

This PR addresses the need for more complex interpolation methods for time varying z dimensions, while maintaining flexibility for users and alignment with the current interpolator API.

The key change is the introduction of a VerticalCoordinateInterpolator abstract class that sits at the same level as the Vector and Scalar interpolator classes. At its core, the VerticalCoordinateInterpolator class enables any existing or user madeScalarInterpolator to be used to interpolate the vertical coordinate. This is done through its base method _columns_from_scalar, which takes in a scalar interpolator and uses it to run a horizontal interpolation on every vertical interface at each particle's horizontal position. These interfaces are then used to reconstruct the vertical columns at each particle position, which can be used for the vertical grid search.

For time varying z coordinate datasets, a vertical coordinate interpolator is set during Fieldset creation, and can be overwritten by the user.

This PR also adjusts all UxInterpolators to perform linear in time interpolation for data that falls between two model snapshots. This is required to avoid large temporal jump discontinuities in time varying vertical grids.

Checklist

Tests added
This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

This PR contains AI-generated content.
- I have tested any AI-generated content in my PR.
- I take responsibility for any AI-generated content in my PR.
- Describe how you used it (e.g., by pasting your prompt): I guided claude opus 5.5 through the implementation. I have thoroughly reviewed all code, and made edits where needed.

…polator

for more information, see https://pre-commit.ci

…owed array

…or UxGrid

for more information, see https://pre-commit.ci

…s-code/Parcels into vertical-coordinate-interpolators

for more information, see https://pre-commit.ci

erikvansebille

Looks good! Some comments below

erikvansebille · 2026-11-05T06:16:41Z

+        warnings.warn(
+            "Time-varying (3D) z coordinates are experimental and may cause significant memory overhead that "
+            f"leads to OOM errors. This z coordinate has sizes {dict(z.sizes)} ({z.nbytes / 1e9:.3g} GB). "
+            f"Assumptions: z has dims {z.dims}, i.e. it is defined at the layer interfaces ('zf') on the mesh "
+            "nodes or faces, and is strictly increasing along 'zf'; each particle's z column is interpolated with "
+            "the grid's `z_interp_method` (by default barycentrically from its face's nodes for node z, or "
+            "constant over its face for face z), and linearly in time between z snapshots.",
+            FieldSetWarning,
+            stacklevel=5,


But I thought this warning may not be necessary anymore because the 4D depth fields can use the ChunkCachedArray or WindowedArray backend?

This is true, I had kept it because we currently don't set it to be ChunkCachedArray or WindowedArray explicitly, but we should probably do this in _init_time_varying_z and then remove the warning.

Added in cb13494

erikvansebille · 2026-11-05T06:17:59Z

+                f"A 3D z must have dims ('time', 'zf', 'n_node') or ('time', 'zf', 'n_face'), got {z.dims}"
+            )
+
+        if z.dims[-1] == "n_node":


Do we also need to check that the dimension order of z is such that the last is n_node or n_face? Or is that done somewhere else?

The check above this line verifies the dimension ordering and throws an error if the dimension ordering doesn't have n_node/n_face last. That being said, in retrospect, and comparing to the XGrid side, I don't think we should error here, and instead add a function similar to preprocess_sgrid_model_data, which would rearrange the dims to guarantee this order but allow us to ingest data of any order.

I'll add that to this PR.

erikvansebille · 2026-11-05T06:27:53Z

+            isinstance(value, XVerticalCoordinateConstant)
+            and get_dim_position(self.sgrid_metadata, depth_dims[-1]) != "face"
+        ):
+            raise ValueError(f"XVerticalCoordinateConstant needs depth on the cell centres, got dims {depth_dims}")


Can we be more explicit here what the cell centres are called in this dataset? Would we know that?

erikvansebille · 2026-11-05T06:29:31Z


        if "Z" in self.axes and ds.depth.ndim == 4:
-            zi, zeta = self._search_time_varying_depth(z, yi, xi, ti, tau)
+            zi, zeta = self._search_time_varying_depth(z, yi, eta, xi, xsi, ti, tau)


I only realise this now, but would it not make more sense to use the order ti, tau, yi, eta, xi, xsi for the position here - in order to harmonise dimensionality throughout the codebase?

erikvansebille · 2026-11-11T07:33:52Z

I can't comment on individual lines, so below some comments on the notebook

rename the title to use lowercase `time varying'?

use the hide-input tag for the dataset creation cell, so that users don't get overwhelmed by it? It's not that relevant for most users how the dataset is created?

Type-o in coordiante

Is it expected that the particles don't move horizontally on the right panel of the last animation?

erikvansebille · 2026-11-11T07:36:48Z

        self._ds = model_data
        self.mesh = mesh
        self._spatialhash = None
+        self._z_interp_method = None


The _interp_method is set on the Field. Why is the _z_interp_method then set on the grid? Can different Fields on the same grid not have different interpolation methods?

I guess in my mind the time varying z coordinate is still a product of the grid. If interpolation between variables is to be as physically consistent as possible I don't think different fields on the same grid should use different interpolators for z, since it would mean that different fields think a particle is at different positions in the water column. That would be a fundamentally different thing when compared to using different interpolators for different fields that all receive the same particle position from the search.

Additionally from a structural point of view the z_interp_method is used during particle search which is grid functionality rather than field functionality, and it is not clear to me if there is a clean way to expose field specific z interpolators to the particle search.

Good point that these interpolators are essentially for the searching of the particle locations in the grid - not for interpolating the values. So you're right that these belong on the grid.

But perhaps we should make clearer then (in the naming of the functions? in the docstrings?) that these z_interp_methods are indeed fundamentally different from 'normal' interpolators

erikvansebille · 2026-11-11T07:37:48Z

        self._ds = model_data
        self.mesh = mesh
        self._spatialhash = None
+        self._z_interp_method = None


Should we now rename interp_method to h_interp_method or something similar to make explicit that these are horizontal interpolation methods? Or should/could that be done in another pR?

Potentially, though my hope is that anyone who is looking into the internals enough to run into this will also run into the doc string on the VerticalCoordinateInterpolator which states this explicitly. Maybe we try like this and if the naming becomes a problem we change it in another PR?

Given the comment above, I now think that we shouldn't rename the other interpolators, since they do something fundamentally different (interpolating values onto a location, vs finding a location in a grid)

erikvansebille · 2026-11-11T07:40:43Z

+    def __repr__(self):
+        return f"{self.__class__.__name__}(...)"


Should we move this interpolators.__repr__ to the repr_utils.py file, so they are all in one place

erikvansebille · 2026-11-11T07:43:37Z

+    upper_half = bcoord >= 0.5
+    lower_face = np.where(upper_half, face_index, face_index - 1)
+    centre_bcoord = np.where(upper_half, bcoord - 0.5, bcoord + 0.5)
+    below_first_centre = lower_face < 0
+    return np.where(below_first_centre, 0, lower_face), np.where(below_first_centre, 0.0, centre_bcoord)


Should we use the offsets here? Because some models are indexed by the upper-right corner and others by the lower-left corner. Or is that not relevant here?

…s-code/Parcels into vertical-coordinate-interpolators

wyatt-fluidnumerics added 2 commits November 4, 2026 21:34

initial vertical coordinate interpolator plumbing

75a20b6

implement UxLinearNodeLinearZf and corresponding Vertical coord inter…

4e98ebc

…polator

github-project-automation Bot added this to Parcels development Nov 4, 2026

github-project-automation Bot moved this to Backlog in Parcels development Nov 4, 2026

pre-commit-ci Bot and others added 16 commits November 4, 2026 22:08

[pre-commit.ci] auto fixes from pre-commit.com hooks

adc20eb

for more information, see https://pre-commit.ci

Add xgrid interpolators and testing

b00c9d1

[pre-commit.ci] auto fixes from pre-commit.com hooks

7bebb28

for more information, see https://pre-commit.ci

swap time indexing in uxinterpolators to be more compatable with wind…

a8a5de2

…owed array

make all ux interpolators linear in time

dd5f01d

Move vertical interpolator testing to its own file

c016552

Add default z interpolator and linear in time interpolation testing f…

3cc828e

…or UxGrid

[pre-commit.ci] auto fixes from pre-commit.com hooks

2da08b7

for more information, see https://pre-commit.ci

Remove experimental warning on startup

6f5355e

Merge branch 'vertical-coordinate-interpolators' of github.com:Parcel…

3850d9e

…s-code/Parcels into vertical-coordinate-interpolators

Renamed tutorial and added vertical coord interpolator section

24ecc04

Clear notebook output

012c0fe

minor changes

13aaeb3

add py:obj to tutorial interpolators

87e87c3

Merge branch 'main' into vertical-coordinate-interpolators

71f206a

[pre-commit.ci] auto fixes from pre-commit.com hooks

4a0b6cc

for more information, see https://pre-commit.ci

wyatt-fluidnumerics marked this pull request as ready for review November 10, 2026 20:33

wyatt-fluidnumerics requested review from erikvansebille and fluidnumericsJoe November 10, 2026 20:33

erikvansebille reviewed Nov 11, 2026

View reviewed changes

wyatt-fluidnumerics added 2 commits November 13, 2026 16:17

Added chunked cache array default for dask backed time varying z

cb13494

Merge branch 'vertical-coordinate-interpolators' of github.com:Parcel…

494f49f

…s-code/Parcels into vertical-coordinate-interpolators

Conversation

wyatt-fluidnumerics commented Nov 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Checklist

AI Disclosure

Uh oh!

erikvansebille left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wyatt-fluidnumerics Nov 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wyatt-fluidnumerics Nov 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

wyatt-fluidnumerics commented Nov 4, 2026 •

edited

Loading

wyatt-fluidnumerics Nov 13, 2026 •

edited

Loading

wyatt-fluidnumerics Nov 12, 2026 •

edited

Loading