add new methods for geulen by OnnoEbbens · Pull Request #556 · gwmod/nlmod

OnnoEbbens · Mar 27, 2026

Add two new methods for merging geotop geulen into a geotop stratigraphic model.

OnnoEbbens · Mar 28, 2026

This PR is ready for review. Although I think it will be pretty hard to review.

I tried to take into account all possible scenarios for this tool, but still I added a warning because I am not sure the method works for all edge cases (there are so many!).

There are a few assumptions in the code:

A geul unit with a higher number (e.g. 6420) is always below a geul unit with a lower number (e.g. 6010).
A geul unit is never interrupted by another unit in the vertical. Meaning there is never another unit in between geul units when looking at a column. This does not hold the other way around, a geul unit can be both above and below another unit.

In the PR there are now 2 functions split_layers_on_geul and split_layers_on_geul_optimal. As the name suggest I think the last one is better as it is slightly faster and the final layer model has less layers. However the code is definitely more complex and more steps are involved. I mainly left both methods in so we can easily go back to that method if there are problems with the 'optimal' one.

The 15_geotop.ipynb has code to show the result from both methods. This might be the easiest way to check the results. I will remove one of the methods and clean up the notebook once we choose.

Copilot

Pull request overview

Adds support for converting GeoTOP stratigraphy with “geulen” (paleochannels) into a consistent layered model by optionally splitting layers, and introduces a new plotting helper to visualize GeoTOP stratigraphy in cross-sections.

Changes:

Add two new geul-handling algorithms (split_layers_on_geul and split_layers_on_geul_optimal) and integrate them into to_model_layers(method_geulen="split_layers").
Generalize GeoTOP cross-section plotting to a variable-based helper and add geotop_strat_in_cross_section.
Update GeoTOP documentation notebook to demonstrate the new stratigraphy cross-section plot and the new split_layers method.

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 9 comments.

File	Description
`nlmod/read/geotop.py`	Adds layer-splitting logic for geulen and wires it into `to_model_layers`.
`nlmod/plot/plot.py`	Adds stratigraphy cross-section plotting and generalizes legend/cross-section helpers.
`nlmod/plot/__init__.py`	Exposes `geotop_strat_in_cross_section` via the public plotting API.
`docs/data_sources/15_geotop.ipynb`	Demonstrates new plotting and experimental `split_layers` geul handling.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

OnnoEbbens · Mar 30, 2026

@dbrakenhoff and @rubencalje, can you wait with the review? Copilot has some good points that I want to solve first.

…otop_geulen

OnnoEbbens · Apr 30, 2026

@rubencalje & @dbrakenhoff this PR is finally ready for review. Failing test not related to this PR.

The PR adds the function split_layers_on_geul which makes it possible to add geulen as separate layers instead of adding geulen to the layer above or the layer below. This comes at the cost of adding extra layers to the layer model. The idea is easiest to understand by having a look at the cross section below. In that example 4 extra layers are added using the split_layers_on_geul function.

Note: there was already an option "add_as_layer" in the nlmod.read.geotop.to_model_layers but that raised an error if geulen are both above and below other stratigraphic layers. I think we can keep that option for the cases where a geul can be interpreted as a "normal" stratigraphic layer.

units in original gt: 16
layers in gtl1: 14
layers in gtl2: 14
layers in gtl3: 20

dbrakenhoff · May 1, 2026

Really nice work! Approved, but I did run into a small bug maybe, so maybe needs one more check? I wrote some tiny code to test a few examples. It produces these kinds of figures (perhaps the basis of some testing...).

I haven't exhausted all the possible, examples so we probably still need to add a few. But I saw that adding a geul at the bottom is currently causing an index error. Perhaps you can check what's going on for that edge case?

# %%
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

from nlmod.read.geotop import split_layers_on_geul


def plot_stratigraphy(strat_original, strat_modified):

    unique_values = np.unique(np.concatenate([strat_original, strat_modified]))
    boundaries = np.append(unique_values, unique_values[-1] + 1)
    norm = mpl.colors.BoundaryNorm(boundaries, len(unique_values))
    cmap = plt.get_cmap("Dark2", len(unique_values))

    fig, (axl, axr) = plt.subplots(1, 2, sharex=True, sharey=True)
    axl.matshow(strat.squeeze(), norm=norm, cmap=cmap)
    axr.matshow(strat_modified.squeeze(), norm=norm, cmap=cmap)
    axl.set_title("Original Stratigraphy")
    axr.set_title("Modified Stratigraphy")
    # Loop through the array to add text
    for istrat, iax in zip([strat_original, strat_modified], (axl, axr)):
        rows, _, cols = strat_original.shape
        for i in range(rows):
            for j in range(cols):
                val = istrat[i, 0, j]
                # Use white text for dark colors, black for light colors
                color = "white"
                iax.text(
                    j,
                    i,
                    str(val.item()),
                    va="center",
                    ha="center",
                    color=color,
                    fontweight="bold",
                )
    return


# %% example split single layer
strat = xr.DataArray(
    np.array([[0, 0, 0], [0, 10, 10], [0, 0, 0]])[:, None, :],
    coords={"z": [0, -1, -2], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)

# %%
# %% example geul at top
strat = xr.DataArray(
    np.array([[0, 10, 10], [1, 1, 1], [2, 2, 2]])[:, None, :],
    coords={"z": [0, -1, -2], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1, 2]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)

# %% example geul at bottom
strat = xr.DataArray(
    np.array([[0, 0, 0], [1, 1, 1], [2, 10, 10]])[:, None, :],
    coords={"z": [0, -1, -2], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1, 2]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)

# %% example simple add between layers
strat = xr.DataArray(
    np.array([[0, 0, 0], [0, 10, 10], [1, 1, 1]])[:, None, :],
    coords={"z": [0, -1, -2], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)

# %% example add between layers with geul at top and bottom
strat = xr.DataArray(
    np.array([[0, 0, 10], [10, 10, 10], [10, 1, 1]])[:, None, :],
    coords={"z": [0, -1, -2], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)

# %% example 2 geulen
strat = xr.DataArray(
    np.array([[0, 10, 0], [1, 20, 1], [2, 2, 2]])[:, None, :],
    coords={"z": [0, -1, -2], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1, 2]
geulen = [10, 20]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)

…otop_geulen

OnnoEbbens · May 1, 2026

Thanks a lot! I was struggling with the tests myself, your approach is very clean. I added the tests and made some improvements for the edge cases.

I expanded your script and pasted it below in case we ever need to debug this.

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

from nlmod.read.geotop import split_layers_on_geul


def plot_stratigraphy(strat_original, strat_modified):

    unique_values = np.unique(np.concatenate([strat_original, strat_modified]))
    boundaries = np.append(unique_values, unique_values[-1] + 1)
    norm = mpl.colors.BoundaryNorm(boundaries, len(unique_values))
    cmap = plt.get_cmap("Dark2", len(unique_values))

    fig, (axl, axr) = plt.subplots(1, 2, sharex=True, sharey=True)
    axl.matshow(strat.squeeze(), norm=norm, cmap=cmap)
    axr.matshow(strat_modified.squeeze(), norm=norm, cmap=cmap)
    axl.set_title("Original Stratigraphy")
    axr.set_title("Modified Stratigraphy")
    # Loop through the array to add text
    for istrat, iax in zip([strat_original, strat_modified], (axl, axr)):
        rows, _, cols = strat_original.shape
        for i in range(rows):
            for j in range(cols):
                val = istrat[i, 0, j]
                # Use white text for dark colors, black for light colors
                color = "white"
                iax.text(
                    j,
                    i,
                    str(val.item()),
                    va="center",
                    ha="center",
                    color=color,
                    fontweight="bold",
                )
    return

# %% example 2 geulen
strat = xr.DataArray(
    np.array([[0, 10, 0], [1, 20, 1], [2, 2, 2]])[:, None, :],
    coords={"z": [0, -0.5, -1.0], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1, 2]
geulen = [10, 20]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)
print(new_unit_order)

# %% example geul at top
strat = xr.DataArray(
    np.array([[0, 10, 10], [1, 1, 1], [2, 2, 2]])[:, None, :],
    coords={"z": [0, -0.5, -1.0], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1, 2]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)
print(new_unit_order)

# %% example geul between and below
strat = xr.DataArray(
    np.array([[0, 0, 0], [10, 1, 1], [0,10, 10]])[:, None, :],
    coords={"z": [0, -0.5, -1.0], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)
print(new_unit_order)

# %% example split single layer
strat = xr.DataArray(
    np.array([[0, 0, 0], [0, 10, 10], [0, 0, 0]])[:, None, :],
    coords={"z": [0, -0.5, -1.0], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)
print(new_unit_order)

# %% example geul at bottom
strat = xr.DataArray(
    np.array([[0, 0, 0], [1, 1, 1], [2, 10, 10]])[:, None, :],
    coords={"z": [0, -0.5, -1.0], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1, 2]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)
print(new_unit_order)

# %% example simple add between layers
strat = xr.DataArray(
    np.array([[0, 0, 0], [0, 10, 10], [1, 1, 1]])[:, None, :],
    coords={"z": [0, -0.5, -1.0], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)
print(new_unit_order)

# %% example add between layers with geul at top and bottom
strat = xr.DataArray(
    np.array([[0, 0, 10], [10, 10, 10], [10, 1, 1]])[:, None, :],
    coords={"z": [0, -0.5, -1.0], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)
print(new_unit_order)

# %% example geul between and below
strat = xr.DataArray(
    np.array([[0, 0, 0], [10, 1, 1], [0,10, 10], [0,2, 2]])[:, None, :],
    coords={"z": [0, -0.5, -1.0, -1.5], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1,2]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)
print(new_unit_order)

# %% example geul between and below
strat = xr.DataArray(
    np.array([[0, 0, 0], [10, 1, 1], [0,10, 10]])[:, None, :],
    coords={"z": [0, -0.5, -1.0], "y": [100], "x": [100, 200, 300]},
    dims=("z", "y", "x"),
)
units_no_geul = [0, 1]
geulen = [10]

strat_modified, new_unit_order = split_layers_on_geul(strat, units_no_geul, geulen)
plot_stratigraphy(strat, strat_modified)
print(new_unit_order)

dbrakenhoff · May 5, 2026

Alright, looks good! Just a quick final check:

This new layer order is what we would expect right? And I guess it doesn't really matter since 0 and 1 never overlap here? Because it does look a bit odd that 0 is now placed earlier in the new order than 1.

OnnoEbbens · May 5, 2026

Yes, the order is correct and the naming is definitely odd.

It is not easy to change this in the split_layers_on_geul function. I think the easiest would be to change this during the conversion of the unit numbers to stratigraphic names.

OnnoEbbens · May 8, 2026

Fixed it in the to_model_layer function so the layer order in the end result makes more sense.

OnnoEbbens added 4 commits March 27, 2026 16:59

add new methods for geulen

22ac418

ruff

10dc138

codacy

b14a6e6

minor changes

6fa302f

add warning

efcfd22

OnnoEbbens marked this pull request as draft March 28, 2026 16:36

OnnoEbbens marked this pull request as ready for review March 28, 2026 16:36

OnnoEbbens requested review from Copilot, dbrakenhoff and rubencalje March 28, 2026 16:36

Copilot started reviewing on behalf of OnnoEbbens March 28, 2026 16:37 View session

Copilot AI reviewed Mar 28, 2026

View reviewed changes

Update nlmod/plot/plot.py

043a64a

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

OnnoEbbens added 6 commits April 29, 2026 14:14

Update method

2dbf45d

Merge branch 'geotop_geulen' of github.com:ArtesiaWater/nlmod into ge…

bdba7a3

…otop_geulen

remove output nb

76862c0

codacy

8a3ef6d

tiny fix

c98447d

textual

4939d65

better debug message

f7d54e7

dbrakenhoff reviewed May 1, 2026

View reviewed changes

Comment thread nlmod/read/geotop.py

dbrakenhoff reviewed May 1, 2026

View reviewed changes

Comment thread nlmod/read/geotop.py Outdated

dbrakenhoff reviewed May 1, 2026

View reviewed changes

Comment thread nlmod/read/geotop.py

Update nlmod/read/geotop.py

655703a

dbrakenhoff approved these changes May 1, 2026

View reviewed changes

add tests and fix edge cases

593be39

Merge branch 'geotop_geulen' of github.com:ArtesiaWater/nlmod into ge…

c3c7448

…otop_geulen

ValueError instead of assert

6c183a5

fix split layer order

060a9cf

Search code, repositories, users, issues, pull requests...

Conversation

OnnoEbbens commented Mar 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

OnnoEbbens commented Mar 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

OnnoEbbens commented Mar 30, 2026

Uh oh!

OnnoEbbens commented Apr 30, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

dbrakenhoff commented May 1, 2026

Uh oh!

OnnoEbbens commented May 1, 2026

Uh oh!

dbrakenhoff commented May 5, 2026

Uh oh!

OnnoEbbens commented May 5, 2026

Uh oh!

OnnoEbbens commented May 8, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

OnnoEbbens commented Mar 27, 2026 •

edited

Loading

OnnoEbbens commented Mar 28, 2026 •

edited

Loading