ParticleSet default initial z not set correctly (does not respect direction of depth axis)

[erikvansebille]

Parcels version

v4

Description

The current code in main/particleset.py for calculating the initial z of the ParticleSet does not work as expected:

Parcels/src/parcels/_core/particleset.py

Lines 80 to 85 in df9a313

	if z is None:
	minz = 0
	for field in self.fieldset.fields.values():
	if field.grid.depth is not None:
	minz = min(minz, field.grid.depth[0])
	z = np.ones(x.size) * minz

This means that if the Grid.depth is e.g. [1, 5, 10] m, the particles will be initiated at z=0 m, which immediately throws an out of bound error.

Better to initialise zmin= np.inf; but then this still wouldn't work for depths that are negative (e.g. [-1, -5, -10] m). So a proper fix would be to make this code block aware of the direction of the grid.

An alternative (simpler) solution would be to force users to set the initial z explicitly; lik ewe already do for x and y. But this would be an API-change from v3, where depth would be inferred from he grid when users don't specify it.

[wyatt-fluidnumerics]

Is it not the case that parcels requires users to input data with strictly increasing depth dims? So the directionality here is predetermined and we can simply go the zmin=np.inf route?

Parcels/src/parcels/_core/xgrid.py

Lines 190 to 191 in df9a313

	if "Z" in grid.axes:
	assert_valid_depth(ds["depth"])

Parcels/src/parcels/_core/xgrid.py

Lines 581 to 586 in df9a313

	def assert_valid_depth(da_depth):
	if not np.all(np.diff(da_depth.values) > 0):
	raise ValueError(
	f"Depth DataArray {da_depth.name!r} with dims {da_depth.dims} must be strictly increasing. "
	f'HINT: you may be able to use ds.reindex to flip depth - e.g., ds = ds.reindex({da_depth.name}=ds["{da_depth.name}"][::-1])'
	)

[erikvansebille]

We can't simply change to minz = np.inf, because we want the level closest to zero; not per se the maximum. While the depth acts must be strongly increasing, its values can be positive or negative. Most modern ocean models (e.g. MITgcm) have negative depths, but for example NEMO still uses positive depths. And we also want to support atmospheric models, where depth/elevation is positive. Note that this is also related to #2307

I recently even found that there is a vertical_axis option in the copernicusmarine toolbox to give output as either depth (positive increasing downward) or elevation (negative and increasing upward): https://toolbox-docs.marine.copernicus.eu/en/v2.0.0/python-interface.html#copernicusmarine.open_dataset

So, in short, the v3 API returned the depth closest to zero that was present in all fields (except for those with only one depth level).

Thinking about it, pephaps we can use a similar method as FieldSet.time_interval, and then pick the value closest to zero

Parcels/src/parcels/_core/fieldset.py

Lines 123 to 134 in 22e203d

	def time_interval(self):
	"""Returns the valid executable time interval of the FieldSet,
	which is the intersection of the time intervals of all fields
	in the FieldSet.
	"""
	time_intervals = (m.time_interval for m in self.models)
	# Filter out Nones from constant Fields
	time_intervals = [t for t in time_intervals if t is not None]
	if len(time_intervals) == 0: # All fields are constant fields
	return None
	return functools.reduce(lambda x, y: x.intersection(y), time_intervals)