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

.. only:: html

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

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

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

.. _sphx_glr_basic_examples_plot_infer_spatial_correlation.py:


Spatial correlation of errors
=============================

Digital elevation models have errors that are spatially correlated due to instrument or processing effects. Here, we
rely on a non-stationary spatial statistics framework to estimate and model spatial correlations in elevation error.
We use a sum of variogram forms to model this correlation, with stable terrain as an error proxy for moving terrain.

**Reference**: `Hugonnet et al. (2022) <https://doi.org/10.1109/jstars.2022.3188922>`_, Figure 5 and Equations 13–16.

.. GENERATED FROM PYTHON SOURCE LINES 11-15

.. code-block:: default

    import geoutils as gu

    import xdem








.. GENERATED FROM PYTHON SOURCE LINES 17-20

We load a difference of DEMs at Longyearbyen, already coregistered using :ref:`coregistration-nuthkaab` as shown in
the :ref:`sphx_glr_basic_examples_plot_nuth_kaab.py` example. We also load the glacier outlines here corresponding to
moving terrain.

.. GENERATED FROM PYTHON SOURCE LINES 20-23

.. code-block:: default

    dh = xdem.DEM(xdem.examples.get_path("longyearbyen_ddem"))
    glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))








.. GENERATED FROM PYTHON SOURCE LINES 24-27

Then, we run the pipeline for inference of elevation heteroscedasticity from stable terrain (*Note: we pass a*
``random_state`` *argument to ensure a fixed, reproducible random subsampling in this example*). We ask for a fit with
a Gaussian model for short range (as it is passed first), and Spherical for long range (as it is passed second):

.. GENERATED FROM PYTHON SOURCE LINES 27-35

.. code-block:: default

    (
        df_empirical_variogram,
        df_model_params,
        spatial_corr_function,
    ) = xdem.spatialstats.infer_spatial_correlation_from_stable(
        dvalues=dh, list_models=["Gaussian", "Spherical"], unstable_mask=glacier_outlines, random_state=42
    )








.. GENERATED FROM PYTHON SOURCE LINES 36-42

The first output corresponds to the dataframe of the empirical variogram, by default estimated using Dowd's estimator
and the circular sampling scheme of ``skgstat.RasterEquidistantMetricSpace`` (Fig. S13 of Hugonnet et al. (2022)). The
``lags`` columns is the upper bound of spatial lag bins (lower bound of first bin being 0), the ``exp`` column is the
"experimental" variance value of the variogram in that bin, the ``count`` the number of pairwise samples, and
``err_exp`` the 1-sigma error of the "experimental" variance, if more than one variogram is estimated with the
``n_variograms`` parameter.

.. GENERATED FROM PYTHON SOURCE LINES 42-44

.. code-block:: default

    df_empirical_variogram






.. 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>exp</th>
          <th>lags</th>
          <th>count</th>
          <th>err_exp</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>4.664756</td>
          <td>28.284271</td>
          <td>10</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>1</th>
          <td>10.093535</td>
          <td>40.000000</td>
          <td>35</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>2</th>
          <td>12.749274</td>
          <td>56.568542</td>
          <td>45</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>3</th>
          <td>10.508531</td>
          <td>80.000000</td>
          <td>80</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>4</th>
          <td>11.057757</td>
          <td>113.137085</td>
          <td>161</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>5</th>
          <td>11.257549</td>
          <td>160.000000</td>
          <td>323</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>6</th>
          <td>11.856904</td>
          <td>226.274170</td>
          <td>699</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>7</th>
          <td>11.290362</td>
          <td>320.000000</td>
          <td>1233</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>8</th>
          <td>12.570223</td>
          <td>452.548340</td>
          <td>2648</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>9</th>
          <td>12.341127</td>
          <td>640.000000</td>
          <td>5019</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>10</th>
          <td>12.702629</td>
          <td>905.096680</td>
          <td>9659</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>11</th>
          <td>12.733932</td>
          <td>1280.000000</td>
          <td>16745</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>12</th>
          <td>13.092382</td>
          <td>1810.193360</td>
          <td>25450</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>13</th>
          <td>13.194880</td>
          <td>2560.000000</td>
          <td>29348</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>14</th>
          <td>13.508054</td>
          <td>3620.386720</td>
          <td>27663</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>15</th>
          <td>13.287467</td>
          <td>5120.000000</td>
          <td>25641</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>16</th>
          <td>12.907434</td>
          <td>7240.773439</td>
          <td>26203</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>17</th>
          <td>13.402171</td>
          <td>10240.000000</td>
          <td>25447</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>18</th>
          <td>13.957255</td>
          <td>14481.546879</td>
          <td>26786</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>19</th>
          <td>14.954976</td>
          <td>20480.000000</td>
          <td>31023</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>20</th>
          <td>17.591130</td>
          <td>28963.093757</td>
          <td>26464</td>
          <td>NaN</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 45-47

The second output is the dataframe of optimized model parameters (``range``, ``sill``, and possibly ``smoothness``)
for a sum of gaussian and spherical models:

.. GENERATED FROM PYTHON SOURCE LINES 47-49

.. code-block:: default

    df_model_params






.. 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>model</th>
          <th>range</th>
          <th>psill</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>gaussian</td>
          <td>66.532318</td>
          <td>11.962027</td>
        </tr>
        <tr>
          <th>0</th>
          <td>spherical</td>
          <td>28963.093757</td>
          <td>4.206344</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 50-51

The third output is the spatial correlation function with spatial lags, derived from the variogram:

.. GENERATED FROM PYTHON SOURCE LINES 51-58

.. code-block:: default

    for spatial_lag in [0, 100, 1000, 10000, 30000]:
        print(
            "Errors are correlated at {:.1f}% for a {:,.0f} m spatial lag".format(
                spatial_corr_function(spatial_lag) * 100, spatial_lag
            )
        )





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

 .. code-block:: none

    Errors are correlated at 100.0% for a 0 m spatial lag
    Errors are correlated at 25.9% for a 100 m spatial lag
    Errors are correlated at 24.7% for a 1,000 m spatial lag
    Errors are correlated at 13.1% for a 10,000 m spatial lag
    Errors are correlated at 0.0% for a 30,000 m spatial lag




.. GENERATED FROM PYTHON SOURCE LINES 59-60

We can plot the empirical variogram and its model on a non-linear X-axis to identify the multi-scale correlations.

.. GENERATED FROM PYTHON SOURCE LINES 60-68

.. code-block:: default

    xdem.spatialstats.plot_variogram(
        df=df_empirical_variogram,
        list_fit_fun=[xdem.spatialstats.get_variogram_model_func(df_model_params)],
        xlabel="Spatial lag (m)",
        ylabel="Variance of\nelevation differences (m)",
        xscale_range_split=[100, 1000],
    )




.. image-sg:: /basic_examples/images/sphx_glr_plot_infer_spatial_correlation_001.png
   :alt: plot infer spatial correlation
   :srcset: /basic_examples/images/sphx_glr_plot_infer_spatial_correlation_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 69-73

This pipeline will not always work optimally with default parameters: variogram sampling is more robust with a lot of
samples but takes long computing times, and the fitting might require multiple tries for forms and possibly bounds
and first guesses to help the least-squares optimization. **To learn how to tune more parameters and use the
subfunctions, see the gallery example:** :ref:`sphx_glr_advanced_examples_plot_variogram_estimation_modelling.py`!


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

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


.. _sphx_glr_download_basic_examples_plot_infer_spatial_correlation.py:

.. only:: html

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




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

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

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

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


.. only:: html

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

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