Create a Hillshade from a Terrain Raster in Python

This tutorial covers how to create a hillshade from a terrain raster in Python, and demonstrates a few options for visualizing lidar-derived Digital Elevation Models.

First, let's import the required packages:

import dotenv
import os
import numpy as np
import neonutilities as nu
import rasterio as rio
from rasterio.plot import show
import matplotlib.pyplot as plt

As of June 2026, NEON requires an API token for data downloads, to reduce bot scraping and improve user support. Tokens can be generated in NEON data portal user accounts - log in to your account or create one, and go to the API Tokens section. For best practices in storing and using tokens, follow the instructions here. Once you've set up your token as an environment variable, you can load it using the python-dotenv package as follows, optionally specifying the path to the .env file in load_dotenv().

dotenv.load_dotenv()
token = os.environ.get("NEON_TOKEN")

Read in the datasets

Download Lidar Elevation Models from TEAK

To start, we will download the NEON Elevation Models (DTM and DSM) which are provided in geotiff (.tif) format. Use the download_url function below to download the data directly from the cloud storage location.

For more information on these data products, refer to the NEON Data Portal page, linked here: Elevation - LiDAR.

# download the DSM and DTM data to the C:/data directory - change this if desired
nu.by_tile_aop(dpid='DP3.30024.001',
               site='TEAK',
               year=2021,
               easting=320000,
               northing=4092000,
               token=token,
               savepath=r'C:\data')

Provisional NEON data are not included. To download provisional data, use input parameter include_provisional=True.


Continuing will download 3 NEON data files totaling approximately 6.7 MB. Do you want to proceed? (y/n)  y


Downloading 3 NEON data files totaling approximately 6.7 MB

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 3/3 [00:01<00:00,  2.97it/s]

Display the DSM and DTM files that have been downloaded:

dem_dir = os.path.expanduser(r"C:\data\DP3.30024.001")

for root, dirs, files in os.walk(dem_dir):
    for file in files:
        if file.endswith('DTM.tif'):
            dtm_file = os.path.join(root, file)
            print(os.path.join(root, file).replace(os.path.expanduser('~/Downloads/'),'..'))
        if file.endswith('DSM.tif'):
            dsm_file = os.path.join(root, file)
            print(os.path.join(root, file).replace(os.path.expanduser('~/Downloads/'),'..'))

C:\data\DP3.30024.001\neon-aop-products\2021\FullSite\D17\2021_TEAK_5\L3\DiscreteLidar\DSMGtif\NEON_D17_TEAK_DP3_320000_4092000_DSM.tif
C:\data\DP3.30024.001\neon-aop-products\2021\FullSite\D17\2021_TEAK_5\L3\DiscreteLidar\DTMGtif\NEON_D17_TEAK_DP3_320000_4092000_DTM.tif

Calculate Hillshade

Hillshade is used to visualize the hypothetical illumination value (from 0-255) of each pixel on a surface given a specified light source. To calculate hillshade, we need the zenith (altitude) and azimuth of the illumination source, as well as the slope and aspect of the terrain. The formula for hillshade is:

$$Hillshade = 255.0 * (( cos(zenith_I)*cos(slope_T))+(sin(zenith_I)*sin(slope_T)*cos(azimuth_I-aspect_T))$$

Where all angles are in radians.

For more information about how hillshades work, refer to the ESRI ArcGIS page How Hillshade Works.

# function to caluclate hillshade
def hillshade(array,azimuth,angle_altitude):
    azimuth = 360.0 - azimuth 
    
    x, y = np.gradient(array)
    slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
    aspect = np.arctan2(-x, y)
    azm_rad = azimuth*np.pi/180. #azimuth in radians
    alt_rad = angle_altitude*np.pi/180. #altitude in radians
 
    shaded = np.sin(alt_rad)*np.sin(slope) + np.cos(alt_rad)*np.cos(slope)*np.cos((azm_rad - np.pi/2.) - aspect)
    
    return 255*(shaded + 1)/2

dtm_dataset = rio.open(dtm_file)
dtm_data = dtm_dataset.read(1)

fig, ax = plt.subplots(1, 1, figsize=(6,6))
dtm_map = show(dtm_dataset, title='Digital Terrain Model', adjust=False, cmap='terrain', ax=ax); # adjust=False prevents auto-normalization
show(dtm_dataset, adjust=False, contour=True, ax=ax); # overlay the contours
im = dtm_map.get_images()[0]
fig.colorbar(im, label = 'Elevation (m)', ax=ax) # add a colorbar
ax.ticklabel_format(useOffset=False, style='plain') # turn off scientific notation

Now that we have a function to generate hillshade, we need to read in the DTM raster using rasterio and then calculate hillshade using the hillshade function. We can then plot both.

# Use hillshade function on the DTM data array
hs_data = hillshade(dtm_data,225,45)

fig, ax = plt.subplots(1, 1, figsize=(6,6))
ext = [dtm_dataset.bounds.left, dtm_dataset.bounds.right, dtm_dataset.bounds.bottom, dtm_dataset.bounds.top]
plt.imshow(hs_data, extent=ext)
plt.title('TEAK Hillshade')
ax=plt.gca(); ax.ticklabel_format(useOffset=False, style='plain') #do not use scientific notation 
rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90) #rotate x tick labels 90 degrees

#Overlay transparent hillshade on DTM:
fig, ax = plt.subplots(1, 1, figsize=(6,6))
im1 = plt.imshow(dtm_data,cmap='terrain_r',extent=ext); 
cbar = plt.colorbar(); cbar.set_label('Elevation, m',rotation=270,labelpad=20)
im2 = plt.imshow(hs_data,cmap='Greys',alpha=0.8,extent=ext); #plt.colorbar()
ax=plt.gca(); ax.ticklabel_format(useOffset=False, style='plain') #do not use scientific notation 
rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90) #rotate x tick labels 90 degrees
plt.grid('on'); # plt.colorbar(); 
plt.title('TEAK Hillshade + DTM');

Calculate CHM & Overlay on Top of Hillshade

Canopy Height can be simply calculated by subtracting the Digital Terrain Model from the Digital Surface Model. While NEON's CHM is calculated using a slightly more sophisticated "pit-free" algorithm (see the ATBD linked at the top of this tutorial), in this example, we'll calculate the CHM with the simple difference formula. First, read in the DSM data set, which we previously downloaded into the data folder.

dsm_dataset = rio.open(dsm_file)
dsm_data = dsm_dataset.read(1)

# calculate CHM by differencing the terrain model (DTM) from the surface model (DSM):
chm_data = dsm_data - dtm_data;

Plot the Canopy Height Model for reference:

fig, ax = plt.subplots(1, 1, figsize=(6,6))
im1 = plt.imshow(chm_data,cmap='Greens',extent=ext); 
ax=plt.gca(); ax.ticklabel_format(useOffset=False, style='plain') #do not use scientific notation 
fig.colorbar(im1, label = 'CHM (m)', ax=ax) # add a colorbar
ax.set_title('Canopy Height Model (DSM-DTM)');

Finally, we can make a plot to bring together all of these visualizations from earlier in the tutorial.

#Overlay transparent hillshade on DTM:
fig, ax = plt.subplots(1, 1, figsize=(10,10))

#Terrain
im1 = plt.imshow(dtm_data,cmap='terrain',extent=ext); 
cbar1 = plt.colorbar(im1,fraction=0.04,pad=0.08,ax=ax); 
cbar1.set_label('Elevation, m',rotation=270,labelpad=20)

#Hillshade
im2 = plt.imshow(hs_data,cmap='Greys',alpha=.5,extent=ext); 

#Canopy
im3 = plt.imshow(chm_data,cmap='Greens',alpha=0.6,extent=ext); 
cbar2 = plt.colorbar(im3,fraction=0.045,pad=0.04,ax=ax); cbar2.set_label('Canopy Height, m',rotation=270,labelpad=15)

ax=plt.gca(); ax.ticklabel_format(useOffset=False, style='plain') #do not use scientific notation 
rotatexlabels = plt.setp(ax.get_xticklabels(),rotation=90) #rotate x tick labels 90 degrees
plt.grid('on'); # plt.colorbar(); 
plt.title('Terrain, Hillshade, & Canopy Height');