Lesson 1. Layer a raster dataset over a hillshade in Python to create a beautiful basemap that represents topography.


Customize Plots of Raster Data - Scientists guide to plotting data in python textbook course module

Welcome to the first lesson in the Customize Plots of Raster Data module. This module covers how overlay rasters to create visualizations and how to make interactive plots.

Chapter Four - Customize Raster Plots

In this chapter, you will learn how to create and customize raster plots in Python using earthpy, matplotlib, and folium.

Learning Objectives

After completing this chapter, you will be able to:

  • Overlay two rasters in Python to create a plot.
  • Create interactive map in Jupyter Notebook using the folium package for Python.
  • Overlay a raster on an interactive map created with folium.
  • Customize a raster map in Python using matplotlib.

What You Need

You need Python and Jupyer Notebook to complete this chapter. You should also have an earth-analytics directory setup on your computer with a data subdirectory within it. You should have completed the lesson on Setting Up the Earth Analytics Python Conda Environment..

Download Colorado Flood Teaching Data Subset data

Overlay Rasters in Python

In this lesson, you will learn about overlaying rasters on top of a hillshade for nicer looking plots in Python.

To overlay a raster, you will plot two different raster datasets in the same plot in matplotlib. You will use alpha to adjust the transparency of one of your rasters so the terrain hillshade gives the raster texture!

You will also turn off the legend for the hillshade plot, as the legend we want to see is the DEM elevation values.

What is a Hillshade?

A hillshade is a representation of the earth’s surface as it would look with shade and shadows from the sun. You often render a hillshade using a greyscale colorramp.

Hillshades make nice underlays for other data as they emphasize the topography visually. This adds depth to your map!

To begin, open up both the Digital Terrain Model and the Digital terrain model hillshade files.

import os
import matplotlib.pyplot as plt
import numpy as np
import rasterio as rio
import earthpy as et
import earthpy.plot as ep

# Import data from EarthPy
data = et.data.get_data('colorado-flood')

# Set working directory
os.chdir(os.path.join(et.io.HOME, 'earth-analytics'))
Downloading from https://ndownloader.figshare.com/files/16371473
Extracted output to /root/earth-analytics/data/colorado-flood/.
# Open raster DTM data
lidar_dem_path = os.path.join("data", "colorado-flood", "spatial", 
                              "boulder-leehill-rd", "pre-flood", 
                              "lidar", "pre_DTM.tif")

with rio.open(lidar_dem_path) as lidar_dem:
    lidar_dem_im = lidar_dem.read(1, masked=True)

# Open dem hillshade
lidar_hs_path = os.path.join("data", "colorado-flood", "spatial", 
                              "boulder-leehill-rd", "pre-flood", 
                              "lidar", "pre_DTM_hill.tif")

with rio.open(lidar_hs_path) as lidar_dem_hill:
    lidar_dem_hill = lidar_dem_hill.read(1, masked=True)

To plot both layers together, you add a alpha value to the dem image. This value makes the image more transparent. Below an alpha of .5 (50%) is applied. Play around with the alpha value to see how it impacts your map.

fig, ax = plt.subplots(figsize=(10, 6))

ep.plot_bands(lidar_dem_im, ax=ax, cmap='viridis_r',
              title="Lidar Digital Elevation Model (DEM)\n overlayed on top of a hillshade")

ax.imshow(lidar_dem_hill, cmap='Greys', alpha=.5)
ax.set_axis_off()

plt.show()
Plot of the Digital Elevation Model overlayed on top of a hillshade.
Plot of the Digital Elevation Model overlayed on top of a hillshade.

Leave a Comment