DEM subtraction

Subtracting one DEM with another should be easy! This is why xdem (with functionality from geoutils) allows directly using the - or + operators on xdem.DEM objects.

Before DEMs can be compared, they need to be reprojected/resampled/cropped to fit the same grid. The xdem.DEM.reproject() method takes care of this.

import matplotlib.pyplot as plt
import xdem

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

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

print(dem_2009)
print(dem_1990)

Out:

Driver:               GTiff
Opened from file:     /home/docs/checkouts/readthedocs.org/user_builds/xdem-erikm/envs/example-gallery/lib/python3.7/site-packages/examples/Longyearbyen/data/DEM_2009_ref.tif
Filename:             /home/docs/checkouts/readthedocs.org/user_builds/xdem-erikm/envs/example-gallery/lib/python3.7/site-packages/examples/Longyearbyen/data/DEM_2009_ref.tif
Raster modified since disk load?  False
Size:                 1332, 985
Number of bands:      1
Data types:           ('float32',)
Coordinate System:    EPSG:25833
NoData Value:         -3.4028234663852886e+38
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-erikm/envs/example-gallery/lib/python3.7/site-packages/examples/Longyearbyen/processed/DEM_1990_coreg.tif
Filename:             /home/docs/checkouts/readthedocs.org/user_builds/xdem-erikm/envs/example-gallery/lib/python3.7/site-packages/examples/Longyearbyen/processed/DEM_1990_coreg.tif
Raster modified since disk load?  False
Size:                 1332, 985
Number of bands:      1
Data types:           ('float32',)
Coordinate System:    EPSG:25833
NoData Value:         -3.4028234663852886e+38
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. xdem.DEM.reproject() is a multi-purpose method that ensures a fit each time:

_ = dem_1990.reproject(dem_2009)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/xdem-erikm/envs/example-gallery/lib/python3.7/site-packages/geoutils/georaster.py:913: UserWarning: Output projection, bounds and size are identical -> return self (not a copy!)
  warnings.warn("Output projection, bounds and size are identical -> return self (not a copy!)")

Oops! xdem just warned us that dem_1990 did not need reprojection, but we asked it to anyway. To hide this prompt, add .reproject(..., silent=True). By default, xdem.DEM.reproject() uses “nearest neighbour” resampling (assuming resampling is needed). This approach is great for speed, but it may often be erroneous when resampling. Other more accurate methods are “bilinear”, “cubic_spline” and others. See geoutils.Raster.reproject() and rasterio.enums.Resampling for more information about reprojection.

Now, we are ready to generate the dDEM:

ddem = dem_2009 - dem_1990

print(ddem)

Out:

Driver:               GTiff
Opened from file:     None
Filename:             /vsimem/5b478c2e-8d8d-4de5-bc98-456cb19cf369/5b478c2e-8d8d-4de5-bc98-456cb19cf369.tif
Raster modified since disk load?  None
Size:                 1332, 985
Number of bands:      1
Data types:           ('float32',)
Coordinate System:    EPSG:25833
NoData Value:         -3.4028234663852886e+38
Pixel Size:           20.0, 20.0
Upper Left Corner:    502810.0, 8654330.0
Lower Right Corner:   529450.0, 8674030.0

It is a new xdem.DEM instance, loaded in memory. Let’s visualize it:

plt.figure(figsize=(8, 5))
plt.imshow(ddem.data.squeeze(), cmap="coolwarm_r", vmin=-20, vmax=20)

plt.show()

… and that’s it! For missing values, xdem provides a number of interpolation methods which are shown in the other examples.