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()
gdal.org/progr...
2) Use the GDAL/OGR Python API and DEMProcessing()
gdal.org/python/
3) Use richdem.TerrainAttribute()
richdem.readth...
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()

[@AfagRizayeva]

Thank you for your videos! It's the third video that has been really helpful for me - appreciate it

[@9lsz2]

Very nice channel !! A few years ago it would have saved me a lot of time and troubles learning on my own. I will definitely recommended it to colleagues and other students. Keep up the good work :)

[@makingsenseremotely6207]

Thank you! I also had many difficulties learning GDAL a couple of years back. So now with this channel, I am trying to reduce the frustration for future students :)

[@seyedpeyman6092]

really really useful , thank you ! I learned a lot things .

[@ridethree]

this is the best channel ever!

[@nibeditasinha290]

Hi, your videos are extremely helpful for processing DEM. Could you please make a video on how to generate flow direction, flow accumulation, river network and watershed from DEM using gdal or richdem

[@payelmallick2580]

can we run this cmd command in jupyter notebook/lab as well?

[@amillatali2589]

thank you so much

[@LuckyLAK17]

Hi many thanks for the video. Very interesting allthough I am not yet quite there with my project. Is there a way how to combine hillshading with colormap to show the terrain relief in way that the shading conveys a 3D effect while looking at a 2D map?

[@shopcrazy3623]

How can I look for the optional arguments, like you used "compute edges" in one of the tool. Is there any documentation about this?

[@aleksandardimovski8589]

Thank you for the extremely useful explanation. I have a minor problem and I really hope someone in this chat could give me a hand. Basically, I have a large DEM file, which I would first like to crop using a smaller polygon and then calculate the slop. For that, I am using the rasterio library and commands rasterio.open and rasterio.mask.
The issue I have is that when opening the original DEM file with gdal.Open, the values are normal, like in your example. However, if I use the file that i cropped or manipulated in any way ( I also tried reprojecting with gdal.Warp) the values don't make sense - It is all -9999 and 0.
Do you have an idea what might be the problem?

[@ivanapetrovska8135]

Can you make a video about DSM or DEM land cover types classification using machine learning algorithms?

[@asdfgh7486]

how can we do aspect analysis and slope analysis?

[@KONDAGENIDEVA]

From where will get dem.tif file?

[@charlesiroegbu8949]

Is it possible to get the dem.tif file you used?

[@XinhLe]

It seems that it does not work on windows 10?

[@makingsenseremotely6207]

It works on Windows 10 - I just tested it. Maybe there is a problem with your GDAL installation?

[@jingli7286]

What's the projection of your dem data?

[@makingsenseremotely6207]

WGS 84/ UTM zone 19S