https://twitter.com/allixender/status/1325843393050710016
The DEM is a pre-clipped and reprojected (in QGIS) GeoTiff from the awesome MERIT-Hydro global hydrologically pre-processed DEM:
http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_Hydro/
It would be quite straight-forward to do the clipping and reprojecting in Rasterio, too (e.g. here).
But now we do an actual fully scripted stream delineation in Python with PySheds:
https://github.com/mdbartos/pysheds
https://www.hatarilabs.com/ih-en/watershed-and-stream-network-delimitation-with-python-and-pysheds-tutorial
Hillshade with Earthpy: https://earthpy.readthedocs.io/en/latest
Plotting it all together wit hsome Seaborn and Matplotlib color tinkering.
The Jupyter notebook (view here) is provided in my GitHub repo.
It would be interesting to redo the exercise with e.g. pygeoprocessing