vertical aggregation#984
Conversation
There was a problem hiding this comment.
Pull request overview
This PR adds vertical (layer-collapsing) aggregation support to the mikeio DataArray/Dataset APIs by introducing axis="z" for mean, min, and max on layered FM 3D geometries, returning a 2D result. It also extends FM 3D geometry reduction to explicitly support reducing over "z", and adds regression tests covering sigma-z behavior, dynamic zn, NaN/inactive layers, and layer subsetting.
Changes:
- Add
axis="z"handling forDataArray.mean/min/maxandDataset.mean/min/max(collapsing layered 3D to 2D; mean is thickness-weighted). - Extend
GeometryFM3D.reduce()to support reducing over"z"and combinations with"element". - Add tests for dynamic sigma stretching, NaN/inactive layers, and layer selection/aggregation.
Reviewed changes
Copilot reviewed 5 out of 8 changed files in this pull request and generated 10 comments.
Show a summary per file
| File | Description |
|---|---|
| tests/test_dataset.py | Adds dataset-level vertical aggregation tests (min/max/mean over "z"), including NaN/inactive layer coverage. |
| tests/test_dataarray.py | Adds dataarray-level tests for thickness-weighted vertical mean with dynamic zn, NaN masking, and layer selection behavior. |
| src/mikeio/spatial/_FM_geometry_layered.py | Adds GeometryFM3D.reduce() support for reducing over "z" (to 2D) and "z"+"element" (to 0D). |
| src/mikeio/dataset/_dataset.py | Adds axis="z" dispatch for dataset min/max/mean to call the corresponding dataarray vertical operations. |
| src/mikeio/dataset/_dataarray.py | Implements vertical mean/min/max for layered FM geometries via _mean_z() and _extremum_z(). |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
… into aggregation_in_vertical
Move the per-element dynamic thickness calculation out of the column loop into a vectorized precomputation step (_compute_dynamic_dz). Previously, _compute_column_dz_dynamic was called for each of the column, resulting in large amount small numpy operations. Now dz is computed for all elements at once in timestep chunks (100 at a time) to limit memory usage, and the column loop just slices into the precomputed array — matching the performance pattern of max(axis="z")
Replace the Python column loop in _extremum_z with a single np.fmin.reduceat / np.fmax.reduceat call. Since column elements are contiguous in memory, reduceat reduces each segment without any loop, padding, or intermediate arrays. Also pass np.fmin/fmax directly from the call sites instead of mapping from np.nanmin/nanmax. speedup x3
|
@miab-dhi please take a look at this. |
|
Hi Henrik,
I saw this and reviewed the PR. It seems like it was already completed, tests passed. Did you want me to test on some unseen data or something?
|
Just wanted to get your opinion if this is useful for you, and make sure that the syntax is optimal, before releasing this into the wild😉. |
|
@tile-timley any comments/suggestions? |
|
I think this is a very useful contribution and the syntax is clear. Looking forward to start getting rid of my own very time and memory limited solutions. Does it work to work when sigma-layers crossing border into z-layers (if sigma depth isn't deep enough for ones purposes)?. Great stuff! |
Yes (and no). Yes: It works with any combination of sigma and z-layers. It performs a vertical aggregation of the layers in the Dataset/DataArray regardless of whether they are sigma, z, or a mix of both. No: Currently it only supports aggregation across full layers — you cannot slice at a specific depth (e.g., "give me everything above -5m"). If you need depth-based extraction, you'd first need to select layers manually before applying this operation. |
|
@otzi5300 I tested this out on a different dataset, 3D sigmaZ dfsu with Phytoplankton results. It worked nicely. I think this a useful addition, is a huge time saver from the traditional ways of achieving this. |
|
Terrific and compared to doing calculation with PFS I think it is fast ("mean" operation included). We often calculate eg. "vertical mean, 0-8 meter depth". If there was a functionality that allowed to "ds.sel(z-range=(z_min, z_max))" in addition to the existing "ds.sel(layers=(layer_min, layer_max))", that would be amazing. This is related to the response from @otzi5300 : |
I created a method to extract data from z-val, which extracts a slice at a depth of water either looking from the bottom or down from the surface. We can discuss it more, maybe find a way to implement it in. |
Added vertical mean, max, min
user can now do:
test
Two test files added:
On simple with 3 sigma and 2 z. No sigma stretching
On more complex with 3 sigma and 2 z. On second timesteps, surface level is increase and sigma streched. Bottom layer in element 2 is nan
Computation speed
the loop: for col_idx in range(n_cols) for sigma-z models is slow.
On a laptop, a mesh with:
number of nodes: 131263
number of elements: 192568
number of timesteps: 9504 records (396 days with hourly data)
takes around 1 minutes to do vertical mean on all-layers for one item
min and max takes only a few seconds