reader_netCDF_CF_generic selects wrong vertical layers when z-coordinate is ascending (altitude data)

[BenSpin]

When using reader_netCDF_CF_generic with a NetCDF file where the z-coordinate is ascending (altitude convention), the reader silently selects incorrect vertical layers. This causes particles to sample wind data from the wrong altitudes.

Problematic Code is reader_netCDF_CF_generic.py method get_variables

# NB: may need to flip if self.z is ascending
indices = np.searchsorted(-self.z, [-z.min(), -z.max()])

np.searchsorted requires the input array to be ascending.

• If self.z is descending (ocean depth convention), then -self.z becomes ascending and works correctly.
• If self.z is ascending (altitude convention), then -self.z becomes descending and returns incorrect indices.

The comment acknowledges this issue, but as far as I can tell no correction is implemented.

This results in a silent failure (no warning or error) and particles receive wind data from incorrect altitude layers. What happens is the simulations appear valid but are physically incorrect.

To get aroudn this you can store the z-coordinate in descending order in the NetCDF file:

[25800, 25500, ..., 0]

A potential correction would be to detect the intended order from the dataset within the reader

if self.z[0] > self.z[-1]:  # descending
    indices = np.searchsorted(-self.z, [-z.min(), -z.max()])
else:  # ascending 
    indices = np.searchsorted(self.z, [z.min(), z.max()])

Environment

• OpenDrift: 1.14.9
• Python: 3.12

[knutfrode]

Hi,

In fact, whereas OpenDrift is designed in a generic way to support 3D drift (or motion) in the atmosphere just as in the ocean, the only atmospheric module in the repository (WindBlow) is just a triviality ("placeholder") without any vertical motion.

Since I have personally no experience with calculating 3D drift in the atmosphere, I worry a little if there are other parts of OpenDrift that would cause similar errors (Errors that go unnoticed are of course the worst).
Thus it would be good to use this opportunity to do some additional testing, in addition to your suggested correction which seems plausible.

Did you make an atmospheric module, since you have this use case?
If so, could you share this module and/or this specific netCDF file (or the first few timesteps if it is large)?
Eventually some unit tests should be made also for atmospheric motion, of which there is presently noting.

Btw, did your own suggested fix provide output as expected?

[BenSpin]

Hi,

Thank you for reply! Yes, the fix I mentioned did resolve the issue, the simulation now runs and produces outputs close to what I would expect.

I did build a custom atmospheric module on top of WindBlow for my specific use case, which is simulating high altitude balloon trajectories. Rather than a general-purpose atmospheric module, its pretty targeted, but maybe its still useful as a test case?

The key elements are:
A PassiveTracer subclass where I added climb rate and target altitude as per-element variables
A WindBlow subclass whose update() moves elements vertically towards their target altitude, then calls update_positions() with horizontal wind components.

I also do a preprocessing step that gets GFS/GEFS GRIB2 data, and then interpolates wind profiles from presure levels onto a fixed altitude grid. This then writes a NetCDF file for reader_netCDF_CF_generic

The fix itself was that the z-axis needed to be stored high to low to match OpenDrifts internal convention, and I compare my results to sondehub for sanity.

I am happy to share the code or a sample netCDF if its useful for adding unit tests.

[knutfrode]

Ok, if your output seems reasonable, it should indicate that the ascending z-coordinate is apparently not a problem elsewhere in the code.

Thus I now added your correction:
e350fb6
We can eventually come back to making some unit test for the atmospheric transport case.