Added vertical profile extractor from 3D dfsu files

[abja-dhi]

extract_vertical function added as a new method for dataset and dataarray objects.
The base input arguments are xs, ys, and mode. xs and ys are the coordinates of the transect polyline vertices. There are two modes for calling this function: discrete or interpolate. According to the selected mode, the following inputs should be provided:

discrete:

• layer_min
• layer_max

interpolate:

• n_horizontal
• n_vertical
• z_min
• z_max

[abja-dhi]

@otzi5300 Could you please review this pull request? @ecomodeller I cannot assign a reviewer to my pull request.

[ecomodeller]

@abja-dhi FYI, we appreciate your contribution, but it will take some time until we can do a proper review of this.

[abja-dhi]

@ecomodeller thank you! I totally understand it requires multiple tests to make sure it works as expected in all different scenarios. Please let me know if you found any issues in the script and I'll do my best to modify the function.

[ecomodeller]

Thanks for this — I've now gone through it properly, and the core of it is good. I checked the discrete mode against tests/testdata/oresund_vertical_slice.dfsu, which turns out to be MIKE Zero's own vertical transect of tests/testdata/oresund_sigma_z.dfsu — same items, same timesteps. Feeding your extract_vertical just the two endpoints of that transect reproduces MIKE's element table, node coordinates and all 441 element values with zero error, and 545 of 550 node z-values bit-exact. Getting the node ordering, element node order and surface-aligned layer numbering right is the hard part, and you have.

A few things to fix before this can go in.

1. Bed nodes get z = 0 at bathymetry steps (sigma-z). In _extract_discrete, each output node's z is resolved from the right-hand column, falling back to the left only when that column is outside the domain. It also needs to fall back when the chosen column simply does not contain that global layer — which happens constantly in sigma-z, where n_layers_per_column varies (4–9 in oresund). When it doesn't, the node keeps node_w = 0 and ends up at z = 0.0. On the transect above that is 5 of 550 nodes, off by 13–23 m, and because the static node_coordinates are taken from zn_out[0] the geometry itself is wrong, not just the plot. A node is only ever created by a column that has the layer, so preferring right and falling back left is sufficient — with that change all 550 nodes match MIKE to float32 precision. A node that still ends up with zero weights should raise rather than silently produce z = 0.

2. Layer numbering is already defined in MIKE IO, and this doesn't follow it. get_layer_elements() documents "layer between 0 (bottom) and n_layers-1 (top), can also be negative counting from -1 at the top layer", and sel(layers=...) / read(layers=...) repeat the same thing plus "top"/"bottom". layer_min/layer_max here are 1-based instead, so layer_min=-1 means "clamped to the bed" — the opposite of what sel(layers=-1) means. They are also silently clamped rather than raising (layer_min=-5, layer_max=99 quietly becomes 1–9, where get_layer_elements would raise). Please reuse the existing convention.

Related: the docstring says "Lowest layer (1 = bed)", but layer ids are surface-aligned, so in a sigma-z file the lowest layer is the deepest possible layer, not each column's own bed. On the transect above, layer_min=1, layer_max=3 returns one element rather than a three-layer band, because only one column is deep enough to have those layers. That is consistent with the rest of MIKE IO, but it needs saying explicitly.

3. ruff format fails on both new files (trailing whitespace, missing final newline), so CI will be red. uv run ruff format .

4. Every interpolate call emits two RuntimeWarning: divide by zero — np.where(denom != 0, (z - zc_lo) / denom, 0.0) evaluates both branches. np.divide(..., out=..., where=denom != 0) avoids it.

5. The tests need to pin some numbers. As they stand they assert types and shapes only, so issue 1 above passes all of them. The reference file makes this cheap: extract along that transect and assert the element count, the element values and zn against oresund_vertical_slice.dfsu directly. Worth adding alongside: the pure-sigma case (basin_3d.dfsu), a layer subset with a pinned element count, a multi-vertex polyline (all the current tests use a straight two-point line), and a dfsu write/re-read round-trip of the discrete result (see tests/test_dfsu2dv.py for that pattern).

Two smaller things while you're in there: GeometryFMVerticalColumn passes the isinstance(geometry, GeometryFM3D) guard (it subclasses it) and returns an all-NaN result instead of raising, and _find_col catches bare Exception where it means OutsideModelDomainError.

Separately, mode="interpolate" only interpolates vertically — values come straight from the containing column, so the horizontal direction is piecewise constant while the layer interfaces are smoothly interpolated. Is that what the MIKE tool does? If so it just needs saying in the docstring.

One API question I'll come back to separately, which doesn't need to hold up the fixes above: the x-axis units for geographic meshes. _resample_polyline measures arc length in raw coordinate units, so for a LONG/LAT mesh the Grid2D dx ends up in degrees, while GeometryFMVerticalProfile.relative_element_distance is geodesic metres and the profile plotter labels its axis "relative distance [m]". The same transect therefore gets two different x axes depending on mode. relative_cumulative_distance(..., is_geo=True) already exists for this.