.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "basic_examples/plot_dem_subtraction.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_dem_subtraction.py: DEM differencing ================ Subtracting a DEM with another one should be easy. xDEM allows to use any operator on :class:`xdem.DEM` objects, such as :func:`+` or :func:`-` as well as most NumPy functions while respecting nodata values and checking that georeferencing is consistent. This functionality is inherited from `GeoUtils' Raster class `_. Before DEMs can be compared, they need to be reprojected to the same grid and have the same 3D CRSs. The :func:`geoutils.Raster.reproject` and :func:`xdem.DEM.to_vcrs` methods are used for this. .. GENERATED FROM PYTHON SOURCE LINES 13-17 .. code-block:: Python import geoutils as gu import xdem .. GENERATED FROM PYTHON SOURCE LINES 18-19 We load two DEMs near Longyearbyen, Svalbard. .. GENERATED FROM PYTHON SOURCE LINES 19-23 .. code-block:: Python dem_2009 = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem")) dem_1990 = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem_coreg")) .. GENERATED FROM PYTHON SOURCE LINES 24-25 We can print the information about the DEMs for a "sanity check". .. GENERATED FROM PYTHON SOURCE LINES 25-29 .. code-block:: Python dem_2009.info() dem_1990.info() .. rst-class:: sphx-glr-script-out .. code-block:: none Driver: GTiff Opened from file: /home/docs/checkouts/readthedocs.org/user_builds/xdem-adehecq/checkouts/latest/examples/data/Longyearbyen/data/DEM_2009_ref.tif Filename: /home/docs/checkouts/readthedocs.org/user_builds/xdem-adehecq/checkouts/latest/examples/data/Longyearbyen/data/DEM_2009_ref.tif Loaded? True Modified since load? False Grid size: 1332, 985 Number of bands: 1 Data types: float32 Coordinate system: ['EPSG:25833'] Nodata value: -9999.0 Pixel interpretation: Area Pixel size: 20.0, 20.0 Upper left corner: 502810.0, 8654330.0 Lower right corner: 529450.0, 8674030.0 Driver: GTiff Opened from file: /home/docs/checkouts/readthedocs.org/user_builds/xdem-adehecq/checkouts/latest/examples/data/Longyearbyen/processed/DEM_1990_coreg.tif Filename: /home/docs/checkouts/readthedocs.org/user_builds/xdem-adehecq/checkouts/latest/examples/data/Longyearbyen/processed/DEM_1990_coreg.tif Loaded? True Modified since load? False Grid size: 1332, 985 Number of bands: 1 Data types: float32 Coordinate system: ['EPSG:25833'] Nodata value: -9999.0 Pixel interpretation: Area Pixel size: 20.0, 20.0 Upper left corner: 502810.0, 8654330.0 Lower right corner: 529450.0, 8674030.0 .. GENERATED FROM PYTHON SOURCE LINES 30-32 In this particular case, the two DEMs are already on the same grid (they have the same bounds, resolution and coordinate system). If they don't, we need to reproject one DEM to fit the other using :func:`geoutils.Raster.reproject`: .. GENERATED FROM PYTHON SOURCE LINES 32-35 .. code-block:: Python dem_1990 = dem_1990.reproject(dem_2009) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/xdem-adehecq/conda/latest/lib/python3.11/site-packages/geoutils/raster/raster.py:2768: UserWarning: Output projection, bounds and grid size are identical -> returning self (not a copy!) warnings.warn( .. GENERATED FROM PYTHON SOURCE LINES 36-42 Oops! GeoUtils just warned us that ``dem_1990`` did not need reprojection. We can hide this output with ``.reproject(..., silent=True)``. By default, :func:`xdem.DEM.reproject` uses "bilinear" resampling (assuming resampling is needed). Other options are detailed at `geoutils.Raster.reproject() `_ and `rasterio.enums.Resampling `_. We now compute the difference by simply substracting, passing `stats=True` to :func:`geoutils.Raster.info` to print statistics. .. GENERATED FROM PYTHON SOURCE LINES 42-47 .. code-block:: Python ddem = dem_2009 - dem_1990 ddem.info(stats=True) .. rst-class:: sphx-glr-script-out .. code-block:: none Driver: None Opened from file: None Filename: None Loaded? True Modified since load? True Grid size: 1332, 985 Number of bands: 1 Data types: float32 Coordinate system: ['EPSG:25833'] Nodata value: -9999.0 Pixel interpretation: Area Pixel size: 20.0, 20.0 Upper left corner: 502810.0, 8654330.0 Lower right corner: 529450.0, 8674030.0 [MAXIMUM]: 49.87 [MINIMUM]: -50.31 [MEDIAN]: -0.36 [MEAN]: -1.18 [STD DEV]: 5.56 .. GENERATED FROM PYTHON SOURCE LINES 48-50 It is a new :class:`xdem.DEM` instance, loaded in memory. Let's visualize it, with some glacier outlines. .. GENERATED FROM PYTHON SOURCE LINES 50-57 .. code-block:: Python # Load the outlines glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines")) glacier_outlines = glacier_outlines.crop(ddem, clip=True) ddem.plot(cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)") glacier_outlines.plot(ref_crs=ddem, fc="none", ec="k") .. image-sg:: /basic_examples/images/sphx_glr_plot_dem_subtraction_001.png :alt: plot dem subtraction :srcset: /basic_examples/images/sphx_glr_plot_dem_subtraction_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 58-59 And we save the output to file. .. GENERATED FROM PYTHON SOURCE LINES 59-61 .. code-block:: Python ddem.save("temp.tif") .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.984 seconds) .. _sphx_glr_download_basic_examples_plot_dem_subtraction.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_dem_subtraction.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_dem_subtraction.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_