Processing DEMs with GDAL in Python

In this tutorial I explain various approaches to calculating terrain attributes such as the slope, aspect and hillshade of a digital elevation model (DEM) in Python. The different options are:
1) Run gdaldem with os.system() or subprocess.call()
https://gdal.org/programs/gdaldem.html
2) Use the GDAL/OGR Python API and DEMProcessing()
https://gdal.org/python/
3) Use richdem.TerrainAttribute()
https://richdem.readthedocs.io/en/lat...

Code:

from osgeo import gdal
import subprocess
import os
import richdem as rd
import matplotlib.pyplot as plt

using os, subprocess
cmd = "gdaldem slope -compute_edges dem.tif slope1.tif" # replace slope with aspect or hillshade
os.system(cmd)
subprocess.check_call(cmd.split())
slp1 = gdal.Open("slope1.tif")
slp1Array = slp1.GetRasterBand(1).ReadAsArray()

using DEMProcessing
dem = gdal.Open("dem.tif")
slp2 = gdal.DEMProcessing("slope2.tif", dem, "slope", computeEdges=True) # replace "slope" with "aspect" or "hillshade"
slp2Array = slp2.GetRasterBand(1).ReadAsArray()

close your datasets!
slp1 = slp2 = dem = None

using richdem
dem = rd.LoadGDAL("dem.tif")
slp3 = rd.TerrainAttribute(dem, attrib="slope_degrees") # replace "slope_degrees" with "slope_riserun", "aspect" ...
rd.SaveGDAL("slope3.tif", slp3)

visualize (example)
plt.figure()
plt.imshow(slp2Array)
plt.colorbar()
plt.show()