Reading Delft3D Flow NetCDF output into parcels

[VeckoTheGecko]

Currently parcels doesn't have the capability to read Delf3D Flow simulation NetCDF flow output. Writing a function like FieldSet.from_nemo() that parses the NetCDF file from Flow into a parcels fieldset would be a useful addition to the codebase for those working with Delft3D.

---

Discussed in https://github.com/OceanParcels/parcels/discussions/1204

^{Originally posted by **VeckoTheGecko** July 27, 2022} I've been working with [Delft3D Flow](https://content.oss.deltares.nl/delft3d/manuals/Delft3D-FLOW_User_Manual.pdf) simulations which have an option to export to NetCDF files. The format of the output NetCDF of `Delf3D Flow` simulations and the input of `parcels` flowfields are different. Are there existing functions (within the codebase, or external tools) to wrangle flow data from `Delf3D Flow` simulations into the `parcels` format? If not, I may develop these functions. Would this be a suitable fit to have somewhere in the `parcels` codebase?

[VeckoTheGecko]

Happy to be assigned this. Will develop it and hopefully have it up within a couple weeks

[VeckoTheGecko]

@erikvansebille I've been working on this feature, and have been making quite a bit of progress. I have also been encountering a roadblock, which I feel the only way to overcome would be through knowing more about the internal workings and design decisions of parcels.

Consideration of "land" cells within parcels

In Delft3D Flow, you specify an arbitrary mesh which represents the water cells for the simulation. Of course, this mesh is not necessarily rectangular, and can stagger towards the shape of a coastline. There is no concept of land cells (whatever is outside the mesh is either land or just outside the domain of the simulation).

The NetCDF output from Delft3D Flow provides the hydrodynamical data, as well as the coordinates of the gridpoints for the mesh. As the NetCDF output is array-like, the shape of the mesh data is (M_MAX, N_MAX) which correspond to the maximum width and height of the mesh. For mesh-points that do not exist, Delft3D fills the corresponding value in the array with 0 (which can be replaced for nan).

My question is, how does parcels handle land cells? Options:

1. Does it assume that the mesh also includes gridpoints for the land cells? (but the U and V currents are always 0 for those land cells)
2. Does it allow the gridpoints of land cells to be undefined, ie. to be nans? (indicating that the gridpoint is not part of the mesh)

If option 1, I would need to generate (interpolate/extrapolate) these missing gridcells before passing it to parcels, which would need to fundamentally work off the assumption that the grid is flat or geopolar curvilinear (as arbitrary curvilinear grids won't have an intuitive way to generate the missing gridpoints).

This leads me into my other question.

Arbitrary curvilinear grids in parcels

Delft3D Flow allows users to specify a completely arbitrary curvilinear grid (eg. for modelling a river). From my understanding parcels has support for flat meshes and for geopolar curvilinear meshes. Are there any limitations within parcels that prevent it from dealing with arbitrary curvilinear grids (eg. interpolation schemes or fieldset plotting)?

---

Let me know if you'd like any code or Delft Flow data to assist with this discussion.

[erikvansebille]

Hi @VeckoTheGecko, sorry for the delay in responding; I was on holidays last week.

To answer your questions:

1. Parcels does not directly have a concept of 'land cells'. Depending on which interpolation scheme is used (see the interpolation tutorial), NaN points in the FieldSet are set to zero. This can sometimes lead to particles getting stuck on land, see the Preventing Stuck particles documentation for more info.
2. Yes, Parcels can work with arbitrary curvilinear grids (see the Parcels v2.0 paper and the NEMO curvilinear tutorial), as long as each grid cell has four neighbors in the horizontal and up to two in the vertical. This means that there is no support (yet) for unstructured/triangular meshes.

[VeckoTheGecko]

Hi @erikvansebille

I was on holidays last week.

All good 😁. Hope you had a good break

NaN points in the FieldSet are set to zero

To clarify, I don't mean nans for currents in the Fields U and V (which would lead to "beaching" of particles). Moreso I mean whether parcels can handle nan values for the (lon, lat) gridpoints themselves. Delft doesn't have land cells in the domain of the mesh, so these gridpoints conceptually just don't exist in Delft.

Whipped up a diagram which should illustrate.

Thanks for the links.

[erikvansebille]

Thanks for the sketch @VeckoTheGecko, that clarifies what you meant.

The answer is 'No', Parcels can't handle NaNs in lot and lat, unfortunately. Because we use a binary search to find the indices of the grid cell where a particle is located, NaNs don' play nice and it may be tricky to adapt the code (see here for the Scipy version of the index search, and here for the C/JIT version).

My hunch is that ti will be easier to preprocess the Fields so that they end up as 'clean' rectangular arrays (your bottom-left sketch)

[VeckoTheGecko]

I've managed to get a working implementation of this up and going for my own work with cartesian meshes.

I'd like to refactor and develop tests for this before submitting a PR (possibly also supporting curvilinear meshes), but am distracted other work at the moment.

I'll aim to get this done and cleaned up by late November/early Dec.

@erikvansebille let me know if you want me to submit a (very rough) draft PR for this in the meantime.

[erikvansebille]

Great to hear that you've got a working implementation! Perhaps you can indeed already upload a rough PR, so that others may also be able to work/advance it. That's what open source development is about 😉