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

.. only:: html

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

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

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

.. _sphx_glr_advanced_examples_plot_heterosc_estimation_modelling.py:


Estimation and modelling of heteroscedasticity
==============================================

Digital elevation models have a precision that can vary with terrain and instrument-related variables. This variability
in variance is called `heteroscedasticy <https://en.wikipedia.org/wiki/Homoscedasticity_and_heteroscedasticity>`_,
and rarely accounted for in DEM studies (see :ref:`intro`). Quantifying elevation heteroscedasticity is essential to
use stable terrain as an error proxy for moving terrain, and standardize data towards a stationary variance, necessary
to apply spatial statistics (see :ref:`spatialstats`).

Here, we show an advanced example in which we look for terrain-dependent explanatory variables to explain the
heteroscedasticity for a DEM difference at Longyearbyen. We use `data binning <https://en.wikipedia.org/wiki/Data_binning>`_
and robust statistics in N-dimension with :func:`xdem.spatialstats.nd_binning`, apply a N-dimensional interpolation with
:func:`xdem.spatialstats.interp_nd_binning`, and scale our interpolant function with a two-step standardization
:func:`xdem.spatialstats.two_step_standardization` to produce the final elevation error function.

**References**: `Hugonnet et al. (2021) <https://doi.org/10.1038/s41586-021-03436-z>`_, Equation 1, Extended Data Fig.
3a and `Hugonnet et al. (2022) <https://doi.org/10.1109/jstars.2022.3188922>`_, Figs. 4 and S6–S9. Equations 7 or 8 can
be used to convert elevation change errors into elevation errors.

.. GENERATED FROM PYTHON SOURCE LINES 21-28

.. code-block:: Python

    import geoutils as gu

    import matplotlib.pyplot as plt
    import numpy as np

    import xdem








.. GENERATED FROM PYTHON SOURCE LINES 30-32

Here, we detail the steps used by ``xdem.spatialstats.infer_heteroscedasticity_from_stable`` exemplified in
:ref:`sphx_glr_basic_examples_plot_infer_heterosc.py`. First, we load example files and create a glacier mask.

.. GENERATED FROM PYTHON SOURCE LINES 32-38

.. code-block:: Python


    ref_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
    dh = xdem.DEM(xdem.examples.get_path("longyearbyen_ddem"))
    glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))
    mask_glacier = glacier_outlines.create_mask(dh)








.. GENERATED FROM PYTHON SOURCE LINES 39-41

We derive terrain attributes from the reference DEM (see :ref:`sphx_glr_basic_examples_plot_terrain_attributes.py`),
which we will use to explore the variability in elevation error.

.. GENERATED FROM PYTHON SOURCE LINES 41-45

.. code-block:: Python

    slope, aspect, planc, profc = xdem.terrain.get_terrain_attribute(
        dem=ref_dem, attribute=["slope", "aspect", "planform_curvature", "profile_curvature"]
    )








.. GENERATED FROM PYTHON SOURCE LINES 46-47

We convert to arrays and keep only stable terrain for the analysis of variability

.. GENERATED FROM PYTHON SOURCE LINES 47-53

.. code-block:: Python

    dh_arr = dh[~mask_glacier].filled(np.nan)
    slope_arr = slope[~mask_glacier].filled(np.nan)
    aspect_arr = aspect[~mask_glacier].filled(np.nan)
    planc_arr = planc[~mask_glacier].filled(np.nan)
    profc_arr = profc[~mask_glacier].filled(np.nan)








.. GENERATED FROM PYTHON SOURCE LINES 54-57

We use :func:`xdem.spatialstats.nd_binning` to perform N-dimensional binning on all those terrain variables, with uniform
bin length divided by 30. We use the NMAD as a robust measure of `statistical dispersion <https://en.wikipedia.org/wiki/Statistical_dispersion>`_
(see :ref:`robuststats-meanstd`).

.. GENERATED FROM PYTHON SOURCE LINES 57-66

.. code-block:: Python


    df = xdem.spatialstats.nd_binning(
        values=dh_arr,
        list_var=[slope_arr, aspect_arr, planc_arr, profc_arr],
        list_var_names=["slope", "aspect", "planc", "profc"],
        statistics=["count", xdem.spatialstats.nmad],
        list_var_bins=30,
    )








.. GENERATED FROM PYTHON SOURCE LINES 67-70

We obtain a dataframe with the 1D binning results for each variable, the 2D binning results for all combinations of
variables and the N-D (here 4D) binning with all variables.
Overview of the dataframe structure for the 1D binning:

.. GENERATED FROM PYTHON SOURCE LINES 70-72

.. code-block:: Python

    df[df.nd == 1]






.. 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>nd</th>
          <th>count</th>
          <th>nmad</th>
          <th>slope</th>
          <th>aspect</th>
          <th>planc</th>
          <th>profc</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>1</td>
          <td>146414.0</td>
          <td>1.960866</td>
          <td>[0.0, 1.9261366800016269)</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>1</th>
          <td>1</td>
          <td>65489.0</td>
          <td>2.287088</td>
          <td>[1.9261366800016269, 3.8522733600032537)</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>2</th>
          <td>1</td>
          <td>65545.0</td>
          <td>2.174368</td>
          <td>[3.8522733600032537, 5.77841004000488)</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>3</th>
          <td>1</td>
          <td>65616.0</td>
          <td>2.158894</td>
          <td>[5.77841004000488, 7.704546720006507)</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>4</th>
          <td>1</td>
          <td>67128.0</td>
          <td>2.169402</td>
          <td>[7.704546720006507, 9.630683400008134)</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>25</th>
          <td>1</td>
          <td>124.0</td>
          <td>3.435265</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>[4.952611862982318, 5.6439102379436346)</td>
        </tr>
        <tr>
          <th>26</th>
          <td>1</td>
          <td>48.0</td>
          <td>3.653727</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>[5.6439102379436346, 6.335208612904951)</td>
        </tr>
        <tr>
          <th>27</th>
          <td>1</td>
          <td>22.0</td>
          <td>2.731598</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>[6.335208612904951, 7.026506987866272)</td>
        </tr>
        <tr>
          <th>28</th>
          <td>1</td>
          <td>6.0</td>
          <td>3.132506</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>[7.026506987866272, 7.717805362827589)</td>
        </tr>
        <tr>
          <th>29</th>
          <td>1</td>
          <td>2.0</td>
          <td>0.175111</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>[7.717805362827589, 8.409103737788904)</td>
        </tr>
      </tbody>
    </table>
    <p>120 rows × 7 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 73-74

And for the 4D binning:

.. GENERATED FROM PYTHON SOURCE LINES 74-76

.. code-block:: Python

    df[df.nd == 4]






.. 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>nd</th>
          <th>count</th>
          <th>nmad</th>
          <th>slope</th>
          <th>aspect</th>
          <th>planc</th>
          <th>profc</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>4</td>
          <td>0.0</td>
          <td>NaN</td>
          <td>[0.0, 1.9261366800016269)</td>
          <td>[0.0, 11.99999532134599)</td>
          <td>[-11.621762873481751, -10.770174675527628)</td>
          <td>[-12.329847511050636, -11.638549136089317)</td>
        </tr>
        <tr>
          <th>1</th>
          <td>4</td>
          <td>0.0</td>
          <td>NaN</td>
          <td>[0.0, 1.9261366800016269)</td>
          <td>[0.0, 11.99999532134599)</td>
          <td>[-11.621762873481751, -10.770174675527628)</td>
          <td>[-11.638549136089317, -10.947250761128)</td>
        </tr>
        <tr>
          <th>2</th>
          <td>4</td>
          <td>0.0</td>
          <td>NaN</td>
          <td>[0.0, 1.9261366800016269)</td>
          <td>[0.0, 11.99999532134599)</td>
          <td>[-11.621762873481751, -10.770174675527628)</td>
          <td>[-10.947250761128, -10.255952386166681)</td>
        </tr>
        <tr>
          <th>3</th>
          <td>4</td>
          <td>0.0</td>
          <td>NaN</td>
          <td>[0.0, 1.9261366800016269)</td>
          <td>[0.0, 11.99999532134599)</td>
          <td>[-11.621762873481751, -10.770174675527628)</td>
          <td>[-10.255952386166681, -9.564654011205363)</td>
        </tr>
        <tr>
          <th>4</th>
          <td>4</td>
          <td>0.0</td>
          <td>NaN</td>
          <td>[0.0, 1.9261366800016269)</td>
          <td>[0.0, 11.99999532134599)</td>
          <td>[-11.621762873481751, -10.770174675527628)</td>
          <td>[-9.564654011205363, -8.873355636244046)</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>809995</th>
          <td>4</td>
          <td>0.0</td>
          <td>NaN</td>
          <td>[55.85796372004718, 57.78410040004881)</td>
          <td>[347.9998643190337, 359.9998596403797)</td>
          <td>[13.074294867187836, 13.925883065141957)</td>
          <td>[4.952611862982318, 5.6439102379436346)</td>
        </tr>
        <tr>
          <th>809996</th>
          <td>4</td>
          <td>0.0</td>
          <td>NaN</td>
          <td>[55.85796372004718, 57.78410040004881)</td>
          <td>[347.9998643190337, 359.9998596403797)</td>
          <td>[13.074294867187836, 13.925883065141957)</td>
          <td>[5.6439102379436346, 6.335208612904951)</td>
        </tr>
        <tr>
          <th>809997</th>
          <td>4</td>
          <td>0.0</td>
          <td>NaN</td>
          <td>[55.85796372004718, 57.78410040004881)</td>
          <td>[347.9998643190337, 359.9998596403797)</td>
          <td>[13.074294867187836, 13.925883065141957)</td>
          <td>[6.335208612904951, 7.026506987866272)</td>
        </tr>
        <tr>
          <th>809998</th>
          <td>4</td>
          <td>0.0</td>
          <td>NaN</td>
          <td>[55.85796372004718, 57.78410040004881)</td>
          <td>[347.9998643190337, 359.9998596403797)</td>
          <td>[13.074294867187836, 13.925883065141957)</td>
          <td>[7.026506987866272, 7.717805362827589)</td>
        </tr>
        <tr>
          <th>809999</th>
          <td>4</td>
          <td>0.0</td>
          <td>NaN</td>
          <td>[55.85796372004718, 57.78410040004881)</td>
          <td>[347.9998643190337, 359.9998596403797)</td>
          <td>[13.074294867187836, 13.925883065141957)</td>
          <td>[7.717805362827589, 8.409103737788904)</td>
        </tr>
      </tbody>
    </table>
    <p>810000 rows × 7 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 77-81

We can now visualize the results of the 1D binning of the computed NMAD of elevation differences with each variable
using :func:`xdem.spatialstats.plot_1d_binning`.
We can start with the slope that has been long known to be related to the elevation measurement error (e.g.,
`Toutin (2002) <https://doi.org/10.1109/TGRS.2002.802878>`_).

.. GENERATED FROM PYTHON SOURCE LINES 81-85

.. code-block:: Python

    xdem.spatialstats.plot_1d_binning(
        df, var_name="slope", statistic_name="nmad", label_var="Slope (degrees)", label_statistic="NMAD of dh (m)"
    )




.. image-sg:: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_001.png
   :alt: plot heterosc estimation modelling
   :srcset: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 86-90

We identify a clear variability, with the dispersion estimated from the NMAD increasing from ~2 meters for nearly flat
slopes to above 12 meters for slopes steeper than 50°.

What about the aspect?

.. GENERATED FROM PYTHON SOURCE LINES 90-93

.. code-block:: Python


    xdem.spatialstats.plot_1d_binning(df, "aspect", "nmad", "Aspect (degrees)", "NMAD of dh (m)")




.. image-sg:: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_002.png
   :alt: plot heterosc estimation modelling
   :srcset: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 94-97

There is no variability with the aspect that shows a dispersion averaging 2-3 meters, i.e. that of the complete sample.

What about the plan curvature?

.. GENERATED FROM PYTHON SOURCE LINES 97-100

.. code-block:: Python


    xdem.spatialstats.plot_1d_binning(df, "planc", "nmad", "Planform curvature (100 m$^{-1}$)", "NMAD of dh (m)")




.. image-sg:: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_003.png
   :alt: plot heterosc estimation modelling
   :srcset: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 101-108

The relation with the plan curvature remains ambiguous.
We should better define our bins to avoid sampling bins with too many or too few samples. For this, we can partition
the data in quantiles in :func:`xdem.spatialstats.nd_binning`.
*Note: we need a higher number of bins to work with quantiles and still resolve the edges of the distribution. As
with many dimensions the ND bin size increases exponentially, we avoid binning all variables at the same
time and instead bin one at a time.*
We define 1000 quantile bins of size 0.001 (equivalent to 0.1% percentile bins) for the profile curvature:

.. GENERATED FROM PYTHON SOURCE LINES 108-118

.. code-block:: Python


    df = xdem.spatialstats.nd_binning(
        values=dh_arr,
        list_var=[profc_arr],
        list_var_names=["profc"],
        statistics=["count", np.nanmedian, xdem.spatialstats.nmad],
        list_var_bins=[np.nanquantile(profc_arr, np.linspace(0, 1, 1000))],
    )
    xdem.spatialstats.plot_1d_binning(df, "profc", "nmad", "Profile curvature (100 m$^{-1}$)", "NMAD of dh (m)")




.. image-sg:: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_004.png
   :alt: plot heterosc estimation modelling
   :srcset: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 119-123

We clearly identify a variability with the profile curvature, from 2 meters for low curvatures to above 4 meters
for higher positive or negative curvature.

What about the role of the plan curvature?

.. GENERATED FROM PYTHON SOURCE LINES 123-133

.. code-block:: Python


    df = xdem.spatialstats.nd_binning(
        values=dh_arr,
        list_var=[planc_arr],
        list_var_names=["planc"],
        statistics=["count", np.nanmedian, xdem.spatialstats.nmad],
        list_var_bins=[np.nanquantile(planc_arr, np.linspace(0, 1, 1000))],
    )
    xdem.spatialstats.plot_1d_binning(df, "planc", "nmad", "Planform curvature (100 m$^{-1}$)", "NMAD of dh (m)")




.. image-sg:: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_005.png
   :alt: plot heterosc estimation modelling
   :srcset: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 134-136

The plan curvature shows a similar relation. Those are symmetrical with 0, and almost equal for both types of curvature.
To simplify the analysis, we here combine those curvatures into the maximum absolute curvature:

.. GENERATED FROM PYTHON SOURCE LINES 136-147

.. code-block:: Python


    maxc_arr = np.maximum(np.abs(planc_arr), np.abs(profc_arr))
    df = xdem.spatialstats.nd_binning(
        values=dh_arr,
        list_var=[maxc_arr],
        list_var_names=["maxc"],
        statistics=["count", np.nanmedian, xdem.spatialstats.nmad],
        list_var_bins=[np.nanquantile(maxc_arr, np.linspace(0, 1, 1000))],
    )
    xdem.spatialstats.plot_1d_binning(df, "maxc", "nmad", "Maximum absolute curvature (100 m$^{-1}$)", "NMAD of dh (m)")




.. image-sg:: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_006.png
   :alt: plot heterosc estimation modelling
   :srcset: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_006.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 148-155

Here's our simplified relation! We now have both slope and maximum absolute curvature with clear variability of
the elevation error.

**But, one might wonder: high curvatures might occur more often around steep slopes than flat slope,
so what if those two dependencies are actually one and the same?**

We need to explore the variability with both slope and curvature at the same time:

.. GENERATED FROM PYTHON SOURCE LINES 155-174

.. code-block:: Python


    df = xdem.spatialstats.nd_binning(
        values=dh_arr,
        list_var=[slope_arr, maxc_arr],
        list_var_names=["slope", "maxc"],
        statistics=["count", np.nanmedian, xdem.spatialstats.nmad],
        list_var_bins=30,
    )

    xdem.spatialstats.plot_2d_binning(
        df,
        var_name_1="slope",
        var_name_2="maxc",
        statistic_name="nmad",
        label_var_name_1="Slope (degrees)",
        label_var_name_2="Maximum absolute curvature (100 m$^{-1}$)",
        label_statistic="NMAD of dh (m)",
    )




.. image-sg:: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_007.png
   :alt: plot heterosc estimation modelling
   :srcset: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_007.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 175-179

We can see that part of the variability seems to be independent, but with the uniform bins it is hard to tell much
more.

If we use custom quantiles for both binning variables, and adjust the plot scale:

.. GENERATED FROM PYTHON SOURCE LINES 179-221

.. code-block:: Python


    custom_bin_slope = np.unique(
        np.concatenate(
            [
                np.nanquantile(slope_arr, np.linspace(0, 0.95, 20)),
                np.nanquantile(slope_arr, np.linspace(0.96, 0.99, 5)),
                np.nanquantile(slope_arr, np.linspace(0.991, 1, 10)),
            ]
        )
    )

    custom_bin_curvature = np.unique(
        np.concatenate(
            [
                np.nanquantile(maxc_arr, np.linspace(0, 0.95, 20)),
                np.nanquantile(maxc_arr, np.linspace(0.96, 0.99, 5)),
                np.nanquantile(maxc_arr, np.linspace(0.991, 1, 10)),
            ]
        )
    )

    df = xdem.spatialstats.nd_binning(
        values=dh_arr,
        list_var=[slope_arr, maxc_arr],
        list_var_names=["slope", "maxc"],
        statistics=["count", np.nanmedian, xdem.spatialstats.nmad],
        list_var_bins=[custom_bin_slope, custom_bin_curvature],
    )
    xdem.spatialstats.plot_2d_binning(
        df,
        "slope",
        "maxc",
        "nmad",
        "Slope (degrees)",
        "Maximum absolute curvature (100 m$^{-1}$)",
        "NMAD of dh (m)",
        scale_var_2="log",
        vmin=2,
        vmax=10,
    )





.. image-sg:: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_008.png
   :alt: plot heterosc estimation modelling
   :srcset: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_008.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 222-234

We identify clearly that the two variables have an independent effect on the precision, with

- *high curvatures and flat slopes* that have larger errors than *low curvatures and flat slopes*
- *steep slopes and low curvatures* that have larger errors than *low curvatures and flat slopes* as well

We also identify that, steep slopes (> 40°) only correspond to high curvature, while the opposite is not true, hence
the importance of mapping the variability in two dimensions.

Now we need to account for the heteroscedasticity identified. For this, the simplest approach is a numerical
approximation i.e. a piecewise linear interpolation/extrapolation based on the binning results available through
the function :func:`xdem.spatialstats.interp_nd_binning`. To ensure that only robust statistic values are used
in the interpolation, we set a ``min_count`` value at 30 samples.

.. GENERATED FROM PYTHON SOURCE LINES 234-239

.. code-block:: Python


    unscaled_dh_err_fun = xdem.spatialstats.interp_nd_binning(
        df, list_var_names=["slope", "maxc"], statistic="nmad", min_count=30
    )








.. GENERATED FROM PYTHON SOURCE LINES 240-244

The output is an interpolant function of slope and curvature that predicts the elevation error at any point. However,
this predicted error might have a spread slightly off from that of the data:

We compare the spread of the elevation difference on stable terrain and the average predicted error:

.. GENERATED FROM PYTHON SOURCE LINES 244-253

.. code-block:: Python

    dh_err_stable = unscaled_dh_err_fun((slope_arr, maxc_arr))

    print(
        "The spread of elevation difference is {:.2f} "
        "compared to a mean predicted elevation error of {:.2f}.".format(
            xdem.spatialstats.nmad(dh_arr), np.nanmean(dh_err_stable)
        )
    )





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

 .. code-block:: none

    The spread of elevation difference is 2.51 compared to a mean predicted elevation error of 2.56.




.. GENERATED FROM PYTHON SOURCE LINES 254-256

Thus, we rescale the function to exactly match the spread on stable terrain using the
:func:`xdem.spatialstats.two_step_standardization` function, and get our final error function.

.. GENERATED FROM PYTHON SOURCE LINES 256-267

.. code-block:: Python


    zscores, dh_err_fun = xdem.spatialstats.two_step_standardization(
        dh_arr, list_var=[slope_arr, maxc_arr], unscaled_error_fun=unscaled_dh_err_fun
    )

    for s, c in [(0.0, 0.1), (50.0, 0.1), (0.0, 20.0), (50.0, 20.0)]:
        print(
            "Elevation measurement error for slope of {:.0f} degrees, "
            "curvature of {:.2f} m-1: {:.1f}".format(s, c / 100, dh_err_fun((s, c))) + " meters."
        )





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

 .. code-block:: none

    Elevation measurement error for slope of 0 degrees, curvature of 0.00 m-1: 2.1 meters.
    Elevation measurement error for slope of 50 degrees, curvature of 0.00 m-1: 5.4 meters.
    Elevation measurement error for slope of 0 degrees, curvature of 0.20 m-1: 8.3 meters.
    Elevation measurement error for slope of 50 degrees, curvature of 0.20 m-1: 14.4 meters.




.. GENERATED FROM PYTHON SOURCE LINES 268-269

This function can be used to estimate the spatial distribution of the elevation error on the extent of our DEMs:

.. GENERATED FROM PYTHON SOURCE LINES 269-273

.. code-block:: Python

    maxc = np.maximum(np.abs(profc), np.abs(planc))
    errors = dh.copy(new_array=dh_err_fun((slope.data, maxc.data)))

    errors.plot(cmap="Reds", vmin=2, vmax=8, cbar_title=r"Elevation error ($1\sigma$, m)")



.. image-sg:: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_009.png
   :alt: plot heterosc estimation modelling
   :srcset: /advanced_examples/images/sphx_glr_plot_heterosc_estimation_modelling_009.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_advanced_examples_plot_heterosc_estimation_modelling.py:

.. only:: html

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

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

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

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

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


.. only:: html

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

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