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 geoutils as gu
import matplotlib.pyplot as plt

import xdem

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”

print(dem_2009)
print(dem_1990)

[[585.1568603515625 593.7670288085938 599.27587890625 ...
  276.78131103515625 287.1056823730469 296.2648620605469]
 [585.667236328125 594.3804321289062 603.4343872070312 ...
  276.9458923339844 288.22772216796875 298.573486328125]
 [584.7866821289062 594.0465087890625 602.7151489257812 ...
  275.3355407714844 285.9632263183594 296.79010009765625]
 ...
 [360.0013122558594 359.192138671875 358.2886657714844 ...
  558.6853637695312 561.381103515625 565.4931640625]
 [360.65240478515625 359.96295166015625 359.01483154296875 ...
  543.9490966796875 546.5569458007812 550.02880859375]
 [361.43194580078125 360.63262939453125 359.747802734375 ...
  529.3741455078125 532.504150390625 537.7479248046875]]
[[589.1026000976562 593.223876953125 598.860595703125 ... 277.697265625
  285.8686218261719 295.0356140136719]
 [590.661376953125 595.2307739257812 601.6672973632812 ...
  277.97552490234375 286.8202209472656 297.48150634765625]
 [589.582763671875 594.1261596679688 601.8610229492188 ...
  277.18865966796875 285.37164306640625 296.2841796875]
 ...
 [382.8089599609375 382.2327575683594 381.5018615722656 ...
  550.736083984375 556.1953735351562 564.4788818359375]
 [383.3118896484375 382.7854309082031 382.1239929199219 ...
  536.96337890625 544.261962890625 555.515625]
 [384.0756530761719 383.80059814453125 383.3339538574219 ...
  526.4791870117188 531.260986328125 541.6585693359375]]

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)

/home/docs/checkouts/readthedocs.org/user_builds/xdem-erikm/conda/latest/lib/python3.11/site-packages/geoutils/raster/raster.py:2099: 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 “bilinear” resampling (assuming resampling is needed). Other options are “nearest” (fast but inaccurate), “cubic_spline”, “lanczos” 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)

[[-3.94573974609375 0.54315185546875 0.415283203125 ... -0.91595458984375
  1.237060546875 1.229248046875]
 [-4.994140625 -0.850341796875 1.76708984375 ... -1.029632568359375
  1.407501220703125 1.09197998046875]
 [-4.79608154296875 -0.07965087890625 0.8541259765625 ...
  -1.853118896484375 0.591583251953125 0.50592041015625]
 ...
 [-22.807647705078125 -23.040618896484375 -23.21319580078125 ...
  7.94927978515625 5.18572998046875 1.0142822265625]
 [-22.65948486328125 -22.822479248046875 -23.109161376953125 ...
  6.9857177734375 2.29498291015625 -5.48681640625]
 [-22.643707275390625 -23.16796875 -23.586151123046875 ...
  2.89495849609375 1.2431640625 -3.91064453125]]

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

ddem.show(cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)")

Let’s add some glacier outlines

# Load the outlines
glacier_outlines = gu.Vector(xdem.examples.get_path("longyearbyen_glacier_outlines"))

# Need to create a common matplotlib Axes to plot both on the same figure
# The xlim/ylim commands are necessary only because outlines extend further than the raster extent
ax = plt.subplot(111)
ddem.show(ax=ax, cmap="coolwarm_r", vmin=-20, vmax=20, cbar_title="Elevation change (m)")
glacier_outlines.ds.plot(ax=ax, fc="none", ec="k")
plt.xlim(ddem.bounds.left, ddem.bounds.right)
plt.ylim(ddem.bounds.bottom, ddem.bounds.top)
plt.title("With glacier outlines")
plt.show()

For missing values, xdem provides a number of interpolation methods which are shown in the other examples.

Saving the output to a file is also very simple

ddem.save("temp.tif")

… and that’s it!