.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "basic_examples/plot_spatial_error_propagation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_basic_examples_plot_spatial_error_propagation.py: Spatial propagation of elevation errors ======================================= Propagating elevation errors spatially accounting for heteroscedasticity and spatial correlation is complex. It requires computing the pairwise correlations between all points of an area of interest (be it for a sum, mean, or other operation), which is computationally intensive. Here, we rely on published formulations to perform computationally-efficient spatial propagation for the mean of elevation (or elevation differences) in an area. **References**: `Hugonnet et al. (2022) `_, Figure S16, Equations 17–19 and `Rolstad et al. (2009) `_, Equation 8. .. GENERATED FROM PYTHON SOURCE LINES 13-20 .. code-block:: default import geoutils as gu import matplotlib.pyplot as plt import numpy as np import xdem .. GENERATED FROM PYTHON SOURCE LINES 22-25 We load the same data, and perform the same calculations on heteroscedasticity and spatial correlations of errors as in the :ref:`sphx_glr_basic_examples_plot_infer_heterosc.py` and :ref:`sphx_glr_basic_examples_plot_infer_spatial_correlation.py` examples. .. GENERATED FROM PYTHON SOURCE LINES 25-34 .. code-block:: default dh = xdem.DEM(xdem.examples.get_path("longyearbyen_ddem")) ref_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem")) glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines")) slope, maximum_curvature = xdem.terrain.get_terrain_attribute(ref_dem, attribute=["slope", "maximum_curvature"]) errors, df_binning, error_function = xdem.spatialstats.infer_heteroscedasticity_from_stable( dvalues=dh, list_var=[slope, maximum_curvature], list_var_names=["slope", "maxc"], unstable_mask=glacier_outlines ) .. GENERATED FROM PYTHON SOURCE LINES 35-37 We use the error map to standardize the elevation differences before variogram estimation, following Equation 12 of Hugonnet et al. (2022), which is more robust as it removes the variance variability due to heteroscedasticity. .. GENERATED FROM PYTHON SOURCE LINES 37-42 .. code-block:: default zscores = dh / errors emp_variogram, params_variogram_model, spatial_corr_function = xdem.spatialstats.infer_spatial_correlation_from_stable( dvalues=zscores, list_models=["Gaussian", "Spherical"], unstable_mask=glacier_outlines, random_state=42 ) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/xdem-erikm/conda/latest/lib/python3.11/site-packages/scipy/optimize/_lsq/trf.py:234: RuntimeWarning: divide by zero encountered in divide Delta = norm(x0 * scale_inv / v**0.5) /home/docs/checkouts/readthedocs.org/user_builds/xdem-erikm/conda/latest/lib/python3.11/site-packages/scipy/optimize/_lsq/trf.py:363: RuntimeWarning: invalid value encountered in scalar divide alpha *= Delta / Delta_new /home/docs/checkouts/readthedocs.org/user_builds/xdem-erikm/conda/latest/lib/python3.11/site-packages/scipy/optimize/_lsq/common.py:49: RuntimeWarning: invalid value encountered in scalar divide t2 = c / q /home/docs/checkouts/readthedocs.org/user_builds/xdem-erikm/conda/latest/lib/python3.11/site-packages/scipy/optimize/_lsq/trf.py:363: RuntimeWarning: invalid value encountered in scalar multiply alpha *= Delta / Delta_new /home/docs/checkouts/readthedocs.org/user_builds/xdem-erikm/conda/latest/lib/python3.11/site-packages/scipy/optimize/_lsq/common.py:398: RuntimeWarning: invalid value encountered in cast return min_step, np.equal(steps, min_step) * np.sign(s).astype(int) .. GENERATED FROM PYTHON SOURCE LINES 43-46 With our estimated heteroscedasticity and spatial correlation, we can now perform the spatial propagation of errors. We select two glaciers intersecting this elevation change map in Svalbard. The best estimation of their standard error is done by directly providing the shapefile, which relies on Equation 18 of Hugonnet et al. (2022). .. GENERATED FROM PYTHON SOURCE LINES 46-57 .. code-block:: default areas = [ glacier_outlines.ds[glacier_outlines.ds["NAME"] == "Brombreen"], glacier_outlines.ds[glacier_outlines.ds["NAME"] == "Medalsbreen"], ] stderr_glaciers = xdem.spatialstats.spatial_error_propagation( areas=areas, errors=errors, params_variogram_model=params_variogram_model ) for glacier_name, stderr_gla in [("Brombreen", stderr_glaciers[0]), ("Medalsbreen", stderr_glaciers[1])]: print(f"The error (1-sigma) in mean elevation change for {glacier_name} is {stderr_gla:.2f} meters.") .. rst-class:: sphx-glr-script-out .. code-block:: none The error (1-sigma) in mean elevation change for Brombreen is 0.80 meters. The error (1-sigma) in mean elevation change for Medalsbreen is 0.87 meters. .. GENERATED FROM PYTHON SOURCE LINES 58-61 When passing a numerical area value, we compute an approximation with disk shape from Equation 8 of Rolstad et al. (2009). This approximation is practical to visualize changes in elevation error when averaging over different area sizes, but is less accurate to estimate the standard error of a certain area shape. .. GENERATED FROM PYTHON SOURCE LINES 61-91 .. code-block:: default areas = 10 ** np.linspace(1, 12) stderrs = xdem.spatialstats.spatial_error_propagation( areas=areas, errors=errors, params_variogram_model=params_variogram_model ) plt.plot(areas / 10**6, stderrs) plt.xlabel("Averaging area (km²)") plt.ylabel("Standard error (m)") plt.vlines( x=np.pi * params_variogram_model["range"].values[0] ** 2 / 10**6, ymin=np.min(stderrs), ymax=np.max(stderrs), colors="red", linestyles="dashed", label="Disk area with radius the\n1st correlation range of {:,.0f} meters".format( params_variogram_model["range"].values[0] ), ) plt.vlines( x=np.pi * params_variogram_model["range"].values[1] ** 2 / 10**6, ymin=np.min(stderrs), ymax=np.max(stderrs), colors="blue", linestyles="dashed", label="Disk area with radius the\n2nd correlation range of {:,.0f} meters".format( params_variogram_model["range"].values[1] ), ) plt.xscale("log") plt.legend() plt.show() .. image-sg:: /basic_examples/images/sphx_glr_plot_spatial_error_propagation_001.png :alt: plot spatial error propagation :srcset: /basic_examples/images/sphx_glr_plot_spatial_error_propagation_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 35.484 seconds) .. _sphx_glr_download_basic_examples_plot_spatial_error_propagation.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_spatial_error_propagation.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_spatial_error_propagation.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_