Comparing modelled heads to observations
This notebook showcases some methods to compare modelled heads to head observations.
from pathlib import Path
import hydropandas as hpd
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr
import nlmod
First read in model results from the IJmuiden model (03_local_grid_refinement.ipynb).
model_ws = Path("../examples/03_local_grid_refinement")
model_name = "IJm_planeten"
# read model dataset
ds = xr.open_dataset(model_ws / f"{model_name}.nc")
ds.attrs["model_ws"] = model_ws
# read heads
head = nlmod.gwf.output.get_heads_da(ds)
Compute the groundwater level and plot the result for the first timestep.
# compute the groundwater level in each time step
gwl = nlmod.gwf.output.get_gwl_from_wet_cells(head)
# plot the heads in the first aquifer
ax = nlmod.plot.map_array(
gwl.isel(time=0), ds=ds, cmap="RdBu_r", colorbar_label="head [m+NAP]"
)
Load the measurements and plot the locations of the observation wells.
# df = pd.read_pickle("./data/20250428_bro_ijmuiden_np1_26_4.pklz", compression="zip")
df = hpd.read_excel("../examples/data/20250428_bro_ijmuiden.xlsx")
df.head()
| x | y | location | filename | source | unit | tube_nr | screen_top | screen_bottom | ground_level | tube_top | metadata_available | obs | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| name | |||||||||||||
| GMW000000004503_001 | 101959.9 | 497438.0 | GMW000000004503 | NaN | BRO | mNAP | 1 | 4.770 | 3.770 | 7.640 | 7.490 | True | GroundwaterObs GMW000000004503_001 -----metada... |
| GMW000000004504_001 | 104691.0 | 497251.5 | GMW000000004504 | NaN | BRO | mNAP | 1 | 0.000 | -1.000 | 2.130 | 2.000 | True | GroundwaterObs GMW000000004504_001 -----metada... |
| GMW000000004505_001 | 100251.4 | 496653.8 | GMW000000004505 | NaN | BRO | mNAP | 1 | 2.936 | 1.936 | 4.795 | 4.696 | True | GroundwaterObs GMW000000004505_001 -----metada... |
| GMW000000007301_001 | 104348.2 | 495624.2 | GMW000000007301 | NaN | BRO | mNAP | 1 | 1.810 | 0.810 | 3.720 | 3.670 | True | GroundwaterObs GMW000000007301_001 -----metada... |
| GMW000000013100_001 | 104554.0 | 499841.0 | GMW000000013100 | NaN | BRO | mNAP | 1 | 1.880 | 0.880 | 3.970 | 3.880 | True | GroundwaterObs GMW000000013100_001 -----metada... |
f, ax = nlmod.plot.get_map(nlmod.grid.get_extent(ds), background=True)
ax.plot(df.x, df.y, "ko");
Get the modeled heads
Get the heads from the cells in which the observation wells are located.
For this we use the nlmod.layers.get_modellayers_indexer() method which takes
a model dataset (defining the model grid) and a dataframe (with the observation
well metadata) as input.
The dataframe must define the x,y locations of the
observation wells, and the top and bottom of the screens. By default it is
assumed these column names follow the hydropandas standard: x, y
screen_top and screen_bottom.
idx = nlmod.layers.get_modellayers_indexer(ds, df)
idx
<xarray.Dataset> Size: 5kB
Dimensions: (name: 114)
Coordinates:
* name (name) object 912B 'GMW000000004503_001' ... 'GMW000000083586_001'
Data variables:
icell2d (name) int64 912B 2569 2796 3409 4624 195 ... 5773 5226 5227 5227
layer (name) <U8 4kB 'NASC' 'NAZA' 'NASC' 'NAZA' ... 'NAZA' 'NAZA' 'NAZA'This indexer can be used directly (if no warnings were raised or if drop_nan_layers=True) to obtain the heads in the cells with observation wells.
Note
If warnings were raised, this means there are observation wells for which the corresponding model layer could not be determined (these probably lie above or below the model). In this case the modellayer is returned as a float array and contains NaNs.
Some post-processing will be necessary to be able to use the indexer e.g. dropping the NaN values:
idx.dropna("name", subset=["layer"])
Additionally, the layer might also have to be renamed to get the layer names corresponding to the layer indices:
idx["layer"].values = ds["layer"].values[idx["layer"].astype(int)]
Try using the indexer to get the modelled heads for each observation well
hsim = head.sel(**idx)
hsim
<xarray.DataArray 'head' (time: 6, name: 114)> Size: 5kB
array([[ 1.09010564, 1.39928643, 1.43258501, 2.85551797, 1.3533678 ,
1.70784255, 1. , 1.67610717, 1.17642116, 2.01915568,
1.28301974, 0.71477995, 2.2319453 , 1.23639789, 1.12483195,
2.73074945, 2.77703465, 1.97130616, 1.92459264, 3.3052092 ,
2.05190367, 2.65102168, 1.7624203 , 3.44951217, 2.10696258,
2.43449446, 2.39039195, 1.48672321, 2.24327515, 1.93501609,
2.84198741, 3.06245864, 1.76152455, 3.12525763, 2.88956555,
2.9754541 , 1.63717556, 1.4790152 , 1.02034932, 3.05452561,
2.02688873, 3.66685363, 1. , 0.92840141, 1.51531802,
1.73170309, 1.80284579, 1.99039628, 2.00050052, 1.80284579,
2.05197211, 1.68635636, 2.87407379, 2.05190367, 3.3052092 ,
2.73074945, 3.85172902, 1.92459264, 2.06842626, 2.02688873,
1.33703625, 1.87005095, 3.4194908 , 1.53744633, 2.06842626,
2.44120241, 2.06833218, 1.71826081, 1.22035302, 1.63920745,
1. , 3.31632905, 1.28921923, 2.03424576, 2.84198741,
1. , 2.46191195, 2.34438063, 3.62343098, 1.0099972 ,
1.8980955 , 1.38976512, 1.27599336, 1. , 1.93582334,
2.03424576, 2.62847225, 1.06025254, 1.99741901, 2.84842844,
2.55693482, 2.05190367, 1. , 2.65100458, 3.44022101,
1.3533678 , 1.18944211, 1.07827447, 1.26945472, 1.44430704,
...
2.2694793 , 2.29630247, 1.58817229, 1.60078215, 2.65306581,
1.79795084, 2.22360158, 1.57941582, 2.78650805, 1.85949189,
2.07410103, 2.00880122, 1.26710019, 1.80979743, 1.61562886,
2.33185057, 2.47698567, 1.46521123, 2.51892283, 2.35328355,
2.4070142 , 1.34196625, 1.19554568, 0.99834095, 2.42925621,
1.79141855, 2.93248157, 1. , 0.91416313, 1.3265559 ,
1.5817896 , 1.54048038, 1.62184451, 1.63258747, 1.54048038,
1.73357427, 1.41636526, 2.36663525, 1.79795084, 2.65306581,
2.2694793 , 3.06237739, 1.60078215, 1.81145578, 1.79141855,
1.13855734, 1.56683194, 2.75231449, 1.31038901, 1.81145578,
1.97083442, 1.78831813, 1.41455254, 0.99280699, 1.48828668,
1. , 2.62425383, 1.1770168 , 1.69446926, 2.33185057,
1. , 2.0602053 , 1.96525277, 2.89144956, 0.99525614,
1.57853711, 1.18693709, 1.09129788, 1. , 1.58483885,
1.69446926, 2.11394479, 0.87294285, 1.64055847, 2.27585644,
2.04268253, 1.79795084, 1. , 2.22062489, 2.76070181,
1.2588423 , 1.13174039, 1.04789346, 1.14875473, 1.27271541,
1.30556235, 1.30556235, 2.13716386, 1.43564494, 1.04419843,
1.13174039, 1.12064264, 1.29200496, 1.31886252, 1.34351739,
1.10531659, 2.07410103, 2.03859008, 2.03859008]])
Coordinates:
* time (time) datetime64[ns] 48B 2015-01-02 2015-01-03 ... 2015-01-07
* name (name) object 912B 'GMW000000004503_001' ... 'GMW00000008358...
y (name) float64 912B 4.974e+05 4.972e+05 ... 4.95e+05 4.95e+05
x (name) float64 912B 1.02e+05 1.046e+05 ... 1.046e+05 1.046e+05
layer (name) <U8 4kB 'NASC' 'NAZA' 'NASC' ... 'NAZA' 'NAZA' 'NAZA'
spatial_ref int64 8B 0
Attributes:
units: m NAPGet and plot the result for the a random observation well.
i = 20
hsim.isel(name=i)
<xarray.DataArray 'head' (time: 6)> Size: 48B
array([2.05190367, 1.10360664, 1.51119841, 0.75373465, 0.92884497,
1.79795084])
Coordinates:
* time (time) datetime64[ns] 48B 2015-01-02 2015-01-03 ... 2015-01-07
y float64 8B 4.948e+05
x float64 8B 1.048e+05
spatial_ref int64 8B 0
layer <U8 32B 'NAZA'
name <U19 76B 'GMW000000022508_001'
Attributes:
units: m NAPhsim.isel(name=i).plot(marker="o", figsize=(10, 3))
# plot observations
df.obs.loc[hsim["name"].values[i]].loc[
pd.Timestamp(hsim.time[0].item()) : pd.Timestamp(hsim.time[-1].item())
].plot(y="value", ax=plt.gca(), marker="x");
/home/docs/checkouts/readthedocs.org/user_builds/nlmod/envs/latest/lib/python3.11/site-packages/pandas/plotting/_matplotlib/core.py:981: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
return ax.plot(*args, **kwds)
hmean = hsim.mean("time")
f, ax = nlmod.plot.get_map(nlmod.grid.get_extent(ds), background=True)
cm = ax.scatter(
hmean.x,
hmean.y,
s=75,
c=hmean.values,
cmap="RdBu_r",
edgecolors="k",
linewidths=0.75,
)
cbar = f.colorbar(cm, ax=ax, label="head [m+NAP]", shrink=0.85)
Interpolating heads
It is also possible to interpolate the heads at the locations of the observation wells. For this we need the original x, y coordinates of the observation wells as well as the layer each well is measuring in.
To get all this information it can be useful to use full_output=True in
nlmod.layers.get_modellayers_indexer(). This returns every variable that is
necessary to compute the layer for each observation well.
Note
This also returns the model layer for the screen_top and screen_bottom separately, allowing you to identify observation wells spanning multiple layers.
idx_full = nlmod.layers.get_modellayers_indexer(ds, df, full_output=True)
idx_full
<xarray.Dataset> Size: 48kB
Dimensions: (name: 114, ilayer: 40)
Coordinates:
* name (name) object 912B 'GMW000000004503_001' ... 'GMW00000008...
Dimensions without coordinates: ilayer
Data variables:
x (name) float64 912B 1.02e+05 1.047e+05 ... 1.047e+05
y (name) float64 912B 4.974e+05 4.973e+05 ... 4.95e+05
screen_top (name) float64 912B 4.77 0.0 2.936 1.81 ... 1.43 1.43 1.42
screen_bottom (name) float64 912B 3.77 -1.0 1.936 0.81 ... 0.43 0.43 0.42
icell2d (name) int64 912B 2569 2796 3409 4624 ... 5226 5227 5227
top (name) float64 912B 7.75 2.25 4.75 3.75 ... 2.75 3.75 3.75
botm (ilayer, name) float64 36kB 7.25 1.75 4.25 ... -221.0 -221.0
modellayer_top (name) float64 912B 1.0 3.0 1.0 3.0 3.0 ... 3.0 3.0 3.0 3.0
modellayer_bot (name) float64 912B 1.0 3.0 1.0 3.0 3.0 ... 3.0 3.0 3.0 3.0
layer (name) <U8 4kB 'NASC' 'NAZA' 'NASC' ... 'NAZA' 'NAZA' 'NAZA'Now we can use nlmod.observations.interpolate_points_ds() to compute the interpolated heads
at each observation well. The first argument is the data array we want to
interpolate. The second argument is a dataset containing information about the
location and layer for each observation well.
We need to pass the correct names for each variable:
x,y: the coordinate names for the locations of the computed heads, the default is"x"and"y"xi,yi: the coordinate names of the observation wells inidx_full, the default is"x"and"y"layer: the layer dimension, the default is “layer”
Our data matches the default so we don’t need to adjust anything.
Note
For structured grids the returned x and y-coordinates in nlmod.layers.get_modellayers_indexer() are the coordinates corresponding to the cell centers. This way the result can be directly used for indexing a data array. The original locations of the observation wells are stored under x_obs y_obs. When using structured grids make sure to pass the correct coordinate names for xi and yi to the interpolate function.
hsim_i = nlmod.observations.interpolate_to_points(head, idx_full)
hsim_i
<xarray.DataArray 'interpolated' (time: 6, name: 114)> Size: 5kB
array([[ 1.08859729e+00, 1.36521920e+00, 1.42751658e+00,
2.85389136e+00, 1.37450418e+00, 1.76154497e+00,
nan, 1.64195433e+00, 1.19219718e+00,
2.07280201e+00, 1.23036662e+00, 9.30510621e-01,
2.26381644e+00, 1.21727456e+00, nan,
2.70654863e+00, 2.77809513e+00, 1.96798861e+00,
1.94790693e+00, 3.33954987e+00, 2.08525086e+00,
2.67232589e+00, nan, 3.38342789e+00,
2.11742862e+00, 2.46980307e+00, 2.50048847e+00,
1.47805184e+00, 2.09834926e+00, 1.89407371e+00,
2.82413937e+00, 3.05591631e+00, 1.77468731e+00,
3.10724708e+00, 2.86047726e+00, 2.94856112e+00,
1.67952780e+00, 1.53767624e+00, 1.04937076e+00,
3.11264800e+00, 2.03097566e+00, 3.66380940e+00,
nan, 9.95164375e-01, 1.41296683e+00,
1.75872794e+00, 1.76225203e+00, 1.99109200e+00,
1.97064926e+00, 1.81611675e+00, 2.02639847e+00,
1.69922668e+00, 2.87118934e+00, 2.08506546e+00,
3.33885099e+00, 2.70632991e+00, 3.93601652e+00,
1.94759966e+00, 2.03463387e+00, 2.02753480e+00,
...
2.67664520e+00, 2.25326040e+00, 3.12115629e+00,
1.61926297e+00, 1.78481012e+00, 1.79129381e+00,
1.14657668e+00, 1.56878268e+00, 2.66875481e+00,
1.30479933e+00, 1.78497901e+00, 1.95173814e+00,
1.78370263e+00, 1.48088561e+00, 9.76096800e-01,
1.55242419e+00, 1.16793263e+00, 2.56136955e+00,
1.11225426e+00, 1.68503166e+00, 2.32187535e+00,
nan, 2.34150870e+00, 1.91622742e+00,
2.86290990e+00, 7.97971884e-01, 1.54601620e+00,
1.18557309e+00, 1.12077942e+00, nan,
1.51191366e+00, 1.69747638e+00, 2.07975944e+00,
8.44669413e-01, 1.66516857e+00, 2.24044108e+00,
2.21606459e+00, 1.80844882e+00, nan,
2.13127501e+00, 2.79910667e+00, 1.26098566e+00,
nan, 1.02394673e+00, 1.15441746e+00,
1.28009684e+00, 1.32470402e+00, 1.32094573e+00,
2.14710639e+00, 1.33778002e+00, 1.05090884e+00,
nan, 1.12347079e+00, 1.29243104e+00,
1.30924184e+00, 1.32922267e+00, 1.11826164e+00,
2.08133773e+00, 2.02348235e+00, 2.03132784e+00]])
Coordinates:
* time (time) datetime64[ns] 48B 2015-01-02 2015-01-03 ... 2015-01-07
* name (name) object 912B 'GMW000000004503_001' ... 'GMW00000008358...
spatial_ref int64 8B 0
Attributes:
units: m NAPCompare the interpolated result to the earlier result.
hsim.isel(name=i).plot(marker="o", figsize=(10, 3))
hsim_i.isel(name=i).plot(marker="o", ax=plt.gca());
Plot the location of the observation well in the grid:
obswell = idx_full.isel(name=i)
extent = [obswell.x - 200, obswell.x + 200, obswell.y - 200, obswell.y + 200]
f, ax = nlmod.plot.get_map(extent, background=True, figsize=6)
nlmod.plot.modelgrid(ds, ax=ax)
ax.plot(head.x, head.y, "k.", label="cell centers")
ax.plot(obswell.x, obswell.y, "ro", markersize=10, label=obswell.name.item())
ax.legend(loc=(0, 1), frameon=False, ncol=2, fontsize="small");
