Head point interpolation#

The starting heads to be used for a model could be based on the interpolation of x-y head measurements.

TIP: In order to have better interpolation results, an area larger than the model domain should be considered.

import matplotlib.pyplot as plt

We’ll start with the usual imports

import numpy as np

import imod

Head measurements information has been obtained from the Dinoloket website for the case study area. This data consists on a .csv file (read using Pandas pd.read_csv ) that contains the following columns: id, time, head, filt_top, filt_bot, elevation, x and y. For this example, all the measurements at different depths and along time for a same id were averaged, to have one reference head value in each point.

heads = imod.data.head_observations()

mean_heads = heads.groupby("id").mean()  # Calculating mean values by id

mean_heads.head(5)
/builds/deltares/imod/imod-python/examples/prepare/point_interpolation.py:37: FutureWarning: The default value of numeric_only in DataFrameGroupBy.mean is deprecated. In a future version, numeric_only will default to False. Either specify numeric_only or select only columns which should be valid for the function.
  mean_heads = heads.groupby("id").mean()  # Calculating mean values by id
Unnamed: 0 Filternummer head filt_top filt_bot Meetpunt tov m NAP x y
id
B12A1745001 390.717647 1.0 8.806310 9.29 8.29 9.86 226410.0 563800.0
B12A1745002 473.572368 2.0 6.731060 6.77 5.77 9.84 226410.0 563800.0
B12A1747001 511.247525 1.0 6.211782 5.64 4.64 7.07 225211.0 562787.0
B12A1747002 269.728261 2.0 6.187935 2.28 1.28 7.00 225211.0 562787.0
B12A1748001 342.303704 1.0 5.590000 2.09 1.09 6.20 224736.0 564392.0


Next we’ll require a grid to interpolate on. iMOD Python has some useful utility functions to generate an empty grid.

xmin = 225_500.0
xmax = 240_000.0
ymin = 559_000.0
ymax = 564_000.0
dx = 100
dy = -100

grid = imod.util.empty_2d(dx, xmin, xmax, dy, ymin, ymax)

Before we can select points in a grid, we first have to remove the points outside the domain.

points_outside_grid = (
    (mean_heads["x"] < xmin)
    | (mean_heads["x"] > xmax)
    | (mean_heads["y"] < ymin)
    | (mean_heads["y"] > ymax)
)
mean_heads_in_grid = mean_heads.loc[~points_outside_grid]

mean_heads_in_grid.head(5)
Unnamed: 0 Filternummer head filt_top filt_bot Meetpunt tov m NAP x y
id
B12A1745001 390.717647 1.0 8.806310 9.29 8.29 9.86 226410.0 563800.0
B12A1745002 473.572368 2.0 6.731060 6.77 5.77 9.84 226410.0 563800.0
B12A1805001 0.000000 1.0 NaN 9.06 8.56 10.36 226134.0 563906.0
B12A1805002 670.838150 2.0 6.472890 5.77 4.77 10.33 226134.0 563906.0
B12A1806001 1359.800000 1.0 9.542414 9.76 9.26 10.96 226509.0 563571.0


The previous head information needs to be assigned to the model grid. imod-python has a tool called imod.select.points_set_values, which assigns values based on x-y coordinates to a previously defined array. In this case, the array is the starting_heads_larger, the values are the mean calculated heads and the x and y are the coordinates corresponding to the heads.

x = mean_heads_in_grid["x"]
y = mean_heads_in_grid["y"]

heads_grid = imod.select.points_set_values(
    grid,
    values=mean_heads_in_grid["head"].to_list(),
    x=x.to_list(),
    y=y.to_list(),
)

# Plotting the points
fig, ax = plt.subplots()
heads_grid.plot.imshow(ax=ax)
dx = 100.0, dy = -100.0
<matplotlib.image.AxesImage object at 0x7ff23fae7250>

The previous information is still only available at certain points, so it needs to be interpolated. The iMOD Python tool imod.prepare.laplace_interpolate will be used to do an interpolation of the previously indicated head values. It is possible to assign interpolation parameters such as the number of iterations and the closing criteria.

interpolated_heads = imod.prepare.laplace_interpolate(
    heads_grid, close=0.001, mxiter=150, iter1=100
)

# Plotting the interpolation results
fig, ax = plt.subplots()
interpolated_heads.plot.imshow(ax=ax)
dx = 100.0, dy = -100.0
<matplotlib.image.AxesImage object at 0x7ff240a22590>

It might be nice to have the locations of the boreholes plotted together with the interpolation, so that we can better judge the quality of the interpolation.

One way to show this is to derive a grid with whether a cell contains an observation. Consequently we can plot this as an overlay, setting cells without an observation as transparent. imshow accepts a matrix with alpha values as well, which allows us to set transparency per cell.

# Derive grid with which cell has an observation
is_observation = (~np.isnan(heads_grid)).astype(np.float64)

fig, ax = plt.subplots()
# Plot the interpolation results
interpolated_heads.plot.imshow(ax=ax)

# We'll plot on the same axis with transparency to use ``is_observation``
# as overlay.
is_observation.plot.imshow(
    ax=ax, add_colorbar=False, cmap="gray_r", alpha=is_observation.values
)
dx = 100.0, dy = -100.0
<matplotlib.image.AxesImage object at 0x7ff23fcfee30>

Total running time of the script: ( 0 minutes 12.660 seconds)

Gallery generated by Sphinx-Gallery