DEM differencing

Subtracting a DEM with another one should be easy.

xDEM allows to use any operator on xdem.DEM objects, such as + or - 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 geoutils.Raster.reproject() and xdem.DEM.to_vcrs() methods are used for this.

import geoutils as gu

import xdem

We load two DEMs near Longyearbyen, Svalbard.

dem_2009 = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))
dem_1990 = xdem.DEM(xdem.examples.get_path("longyearbyen_tba_dem_coreg"))

We can print the information about the DEMs for a “sanity check”.

dem_2009.info()
dem_1990.info()

Driver:               GTiff
Opened from file:     /home/docs/checkouts/readthedocs.org/user_builds/xdem/checkouts/latest/examples/data/Longyearbyen/data/DEM_2009_ref.tif
Filename:             /home/docs/checkouts/readthedocs.org/user_builds/xdem/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/checkouts/latest/examples/data/Longyearbyen/processed/DEM_1990_coreg.tif
Filename:             /home/docs/checkouts/readthedocs.org/user_builds/xdem/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

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 geoutils.Raster.reproject():

dem_1990 = dem_1990.reproject(dem_2009)

/home/docs/checkouts/readthedocs.org/user_builds/xdem/conda/latest/lib/python3.11/site-packages/geoutils/raster/raster.py:2731: UserWarning: Output projection, bounds and grid size are identical -> returning self (not a copy!)
  warnings.warn(

Oops! GeoUtils just warned us that dem_1990 did not need reprojection. We can hide this output with .reproject(..., silent=True). By default, 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 geoutils.Raster.info() to print statistics.

ddem = dem_2009 - dem_1990

ddem.info(stats=True)

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.19
[STD DEV]:          5.56

It is a new xdem.DEM instance, loaded in memory. Let’s visualize it, with some glacier outlines.

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

And we save the output to file.

ddem.save("temp.tif")