.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples\prepare\point_interpolation.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_prepare_point_interpolation.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_prepare_point_interpolation.py:


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.

.. GENERATED FROM PYTHON SOURCE LINES 12-15

.. code-block:: Python


    import matplotlib.pyplot as plt








.. GENERATED FROM PYTHON SOURCE LINES 16-17

We'll start with the usual imports

.. GENERATED FROM PYTHON SOURCE LINES 17-23

.. code-block:: Python

    import numpy as np

    import imod










.. GENERATED FROM PYTHON SOURCE LINES 25-34

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
<https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_csv.html#pandas.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.

.. GENERATED FROM PYTHON SOURCE LINES 34-42

.. code-block:: Python


    heads = imod.data.head_observations()

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

    mean_heads.head(5)







.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>Unnamed: 0</th>
          <th>Filternummer</th>
          <th>time</th>
          <th>head</th>
          <th>filt_top</th>
          <th>filt_bot</th>
          <th>Meetpunt tov m NAP</th>
          <th>x</th>
          <th>y</th>
        </tr>
        <tr>
          <th>id</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>B12A1745001</th>
          <td>390.717647</td>
          <td>1.0</td>
          <td>2012-02-12 01:58:35.294117632</td>
          <td>8.806310</td>
          <td>9.29</td>
          <td>8.29</td>
          <td>9.86</td>
          <td>226410.0</td>
          <td>563800.0</td>
        </tr>
        <tr>
          <th>B12A1745002</th>
          <td>473.572368</td>
          <td>2.0</td>
          <td>2011-11-19 12:18:56.842105344</td>
          <td>6.731060</td>
          <td>6.77</td>
          <td>5.77</td>
          <td>9.84</td>
          <td>226410.0</td>
          <td>563800.0</td>
        </tr>
        <tr>
          <th>B12A1747001</th>
          <td>511.247525</td>
          <td>1.0</td>
          <td>2008-07-30 00:00:00.000000000</td>
          <td>6.211782</td>
          <td>5.64</td>
          <td>4.64</td>
          <td>7.07</td>
          <td>225211.0</td>
          <td>562787.0</td>
        </tr>
        <tr>
          <th>B12A1747002</th>
          <td>269.728261</td>
          <td>2.0</td>
          <td>2007-11-15 09:07:49.565217280</td>
          <td>6.187935</td>
          <td>2.28</td>
          <td>1.28</td>
          <td>7.00</td>
          <td>225211.0</td>
          <td>562787.0</td>
        </tr>
        <tr>
          <th>B12A1748001</th>
          <td>342.303704</td>
          <td>1.0</td>
          <td>2007-12-25 16:10:40.000000000</td>
          <td>5.590000</td>
          <td>2.09</td>
          <td>1.09</td>
          <td>6.20</td>
          <td>224736.0</td>
          <td>564392.0</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 43-45

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

.. GENERATED FROM PYTHON SOURCE LINES 45-55

.. code-block:: Python


    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)








.. GENERATED FROM PYTHON SOURCE LINES 56-58

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

.. GENERATED FROM PYTHON SOURCE LINES 58-70

.. code-block:: Python


    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)







.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>Unnamed: 0</th>
          <th>Filternummer</th>
          <th>time</th>
          <th>head</th>
          <th>filt_top</th>
          <th>filt_bot</th>
          <th>Meetpunt tov m NAP</th>
          <th>x</th>
          <th>y</th>
        </tr>
        <tr>
          <th>id</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>B12A1745001</th>
          <td>390.717647</td>
          <td>1.0</td>
          <td>2012-02-12 01:58:35.294117632</td>
          <td>8.806310</td>
          <td>9.29</td>
          <td>8.29</td>
          <td>9.86</td>
          <td>226410.0</td>
          <td>563800.0</td>
        </tr>
        <tr>
          <th>B12A1745002</th>
          <td>473.572368</td>
          <td>2.0</td>
          <td>2011-11-19 12:18:56.842105344</td>
          <td>6.731060</td>
          <td>6.77</td>
          <td>5.77</td>
          <td>9.84</td>
          <td>226410.0</td>
          <td>563800.0</td>
        </tr>
        <tr>
          <th>B12A1805001</th>
          <td>0.000000</td>
          <td>1.0</td>
          <td>2012-07-20 00:00:00.000000000</td>
          <td>NaN</td>
          <td>9.06</td>
          <td>8.56</td>
          <td>10.36</td>
          <td>226134.0</td>
          <td>563906.0</td>
        </tr>
        <tr>
          <th>B12A1805002</th>
          <td>670.838150</td>
          <td>2.0</td>
          <td>2014-07-16 13:19:04.508670464</td>
          <td>6.472890</td>
          <td>5.77</td>
          <td>4.77</td>
          <td>10.33</td>
          <td>226134.0</td>
          <td>563906.0</td>
        </tr>
        <tr>
          <th>B12A1806001</th>
          <td>1359.800000</td>
          <td>1.0</td>
          <td>2016-05-24 22:24:00.000000000</td>
          <td>9.542414</td>
          <td>9.76</td>
          <td>9.26</td>
          <td>10.96</td>
          <td>226509.0</td>
          <td>563571.0</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 71-77

The previous head information needs to be assigned to the model grid.
imod-python has a tool called :doc:`/api/generated/select/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.

.. GENERATED FROM PYTHON SOURCE LINES 77-92

.. code-block:: Python


    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)




.. image-sg:: /examples/prepare/images/sphx_glr_point_interpolation_001.png
   :alt: dx = 100.0, dy = -100.0
   :srcset: /examples/prepare/images/sphx_glr_point_interpolation_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.image.AxesImage object at 0x00000204D78D3250>



.. GENERATED FROM PYTHON SOURCE LINES 93-99

The previous information is still only available at certain points, so it
needs to be interpolated. The iMOD Python tool
:doc:`/api/generated/prepare/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.

.. GENERATED FROM PYTHON SOURCE LINES 99-108

.. code-block:: Python


    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)




.. image-sg:: /examples/prepare/images/sphx_glr_point_interpolation_002.png
   :alt: dx = 100.0, dy = -100.0
   :srcset: /examples/prepare/images/sphx_glr_point_interpolation_002.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.image.AxesImage object at 0x00000204D97524D0>



.. GENERATED FROM PYTHON SOURCE LINES 109-117

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.

.. GENERATED FROM PYTHON SOURCE LINES 117-131

.. code-block:: Python


    # 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
    )




.. image-sg:: /examples/prepare/images/sphx_glr_point_interpolation_003.png
   :alt: dx = 100.0, dy = -100.0
   :srcset: /examples/prepare/images/sphx_glr_point_interpolation_003.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <matplotlib.image.AxesImage object at 0x00000204CDDF3690>




.. rst-class:: sphx-glr-timing

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


.. _sphx_glr_download_examples_prepare_point_interpolation.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: point_interpolation.ipynb <point_interpolation.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: point_interpolation.py <point_interpolation.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: point_interpolation.zip <point_interpolation.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_