Re-creating a terrain map with Python
In a previous post, I documented the layers of a terrain map that I made using ArcGIS Pro while following along with a great tutorial by John Nelson. I wanted to try to automate the creation of the map using Python libraries.
The starting point is the DEM that I created for the original map.
That DEM by itself is a thing of beauty, but let’s get to layering.
Applying a custom colormap to the DEM
For the initial elevation-style coloring that I did in the ArcGIS Pro version, I used ‘Elevation #7.’ I didn’t see anything close to this scheme in the standard color maps available in matplotlib, so I decided to create a custom colormap based on it.
First I took a screen capture of the ArcGIS Pro colormap and saved it as a PNG.
Next I wrote the following code to extract the colors from it and return a colormap.
from PIL import Image
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.pyplot as plt
def extract_colors_from_image(image_path: str) -> np.ndarray:
"""Extract colors from the center line of an image."""
image = Image.open(image_path).convert('RGB')
width, height = image.size
colors = [image.getpixel((x, int(height / 2))) for x in range(width)]
return np.array(colors) / 255.0
colors = extract_colors_from_image('input/ramp.png')
color_ramp = LinearSegmentedColormap.from_list("custom_ramp", colors)To check it out, we can show the colormap with matplotlib.
plt.imshow([colors], aspect='auto')
plt.title('Custom Color Ramp')
plt.axis('off')
plt.show()
To see how this looks on the area, let’s load the DEM and preview the results.
from osgeo import gdal
import numpy as np
def normalize_array(array: np.ndarray, lower_percentile: float = 0.0, upper_percentile: float = 100.0) -> np.ndarray:
"""Normalize array values to a specified percentile range."""
lower_val = np.percentile(array, lower_percentile)
upper_val = np.percentile(array, upper_percentile)
clipped = np.clip(array, lower_val, upper_val)
return (clipped - lower_val) / (upper_val - lower_val)
dem_path = 'input/oblique-clip.tif'
dem_data = gdal.Open(dem_path).ReadAsArray()
normalized_dem = normalize_array(dem_data)
plt.imshow(normalized_dem, cmap=color_ramp)
plt.colorbar()
plt.title('DEM with Custom Color Ramp')
plt.show()
This looks great. It’s pretty much identical to how things looked in ArcGIS Pro. Let’s go with it.
def apply_colormap(array: np.ndarray, cmap: LinearSegmentedColormap) -> np.ndarray:
"""Apply a matplotlib colormap to an array."""
colormap = plt.get_cmap(cmap) if isinstance(cmap, str) else cmap
normed_data = (array - array.min()) / (array.max() - array.min())
colored = colormap(normed_data)
return (colored[:, :, :3] * 255).astype('uint8')
terrain_dem = apply_colormap(normalized_dem, color_ramp)
terrain_image = Image.fromarray(terrain_dem)
terrain_image.save('output/1-terrain.png')
Warming the terrain
In John Nelson’s tutorial, he used the ‘inferno’ colormap in ArcGIS Pro for his warming layer. The equivalent in matplotlib is the ‘plasma’ color map. We can just use its name to apply it to the DEM.
warm_dem = apply_colormap(normalized_dem, 'plasma')
warm_image = Image.fromarray(warm_dem)
warm_image.save('output/2-warming.png')
Now we’ve got two versions of the DEM: the custom colormap and the plasma version. In ArcGIS Pro, the Soft Light blend mode was used to combine them. Soft light combines the brightness values from the layers and enhances contrast. In Python, we can use the PIL library to do this.
from PIL import ImageChops
warming_blended = ImageChops.soft_light(terrain_image, warm_image)
warming_blended.save('output/3-warming_blended.png')
Now it’s time for some hillshade.
Applying a traditional hillshade
I wrote two functions to handle hillshade generation. The first generates a single hillshade and the second uses that function to generate and combine multiple hillshades (used in the next step for the multidirectional hillshade layer). I’m not sure if this is close to what’s happening under the hood in ArcGIS Pro, but the results were pretty close.
def generate_hillshade(dem_path: str, azimuth: float, altitude: float, zFactor: float = 1.0) -> np.ndarray:
"""Generate a hillshade for a specific azimuth and altitude."""
options = gdal.DEMProcessingOptions(format='GTiff', computeEdges=True, zFactor=zFactor, azimuth=azimuth, altitude=altitude)
hillshade_path = f'output/temp_hillshade_{azimuth}_{altitude}.tif'
gdal.DEMProcessing(hillshade_path, dem_path, 'hillshade', options=options)
return gdal.Open(hillshade_path).ReadAsArray()
def combine_hillshades(dem_path: str, azimuths: List[float], altitudes: List[float], weights: Optional[List[float]] = None) -> np.ndarray:
"""Combine hillshades from multiple directions."""
hillshades = [generate_hillshade(dem_path, az, alt) for az, alt in zip(azimuths, altitudes)]
if weights is None:
weights = np.ones(len(azimuths))
weights = np.array(weights) / np.sum(weights)
combined_hillshade = np.average(hillshades, axis=0, weights=weights)
combined_hillshade = np.clip(combined_hillshade, 0, 255)
return combined_hillshadeI didn’t bother cleaning up the intermediate hillshades because I wanted to inspect them when playing with outputs.
The following generates a traditional (single azimuth and altitude) hillshade and then blends it into the warmed terrain using the overlay blend mode. Overlay will make dark areas darker and light areas lighter.
def rgb_image_from_array(array: np.ndarray) -> Image.Image:
"""Create an RGB image from an array."""
return Image.fromarray((array * 255 / array.max()).astype('uint8')).convert("RGB")
traditional_hillshade = combine_hillshades(dem_path, [315], [45])
traditional_hillshade_image = rgb_image_from_array(traditional_hillshade)
traditional_hillshade_image.save('output/4-traditional_hillshade.png')
traditional_hillshade_blended = ImageChops.overlay(warming_blended, traditional_hillshade_image)
traditional_hillshade_blended.save('output/5-hillshade_blended.png')
Looking good and still tracking to the ArcGIS Pro version.
Applying a multidirectional hillshade
This step uses the same functions as the previous step except this time we feed in four different azimuths and altitudes to be averaged together. Multidirectional hillshade simulates the real world a bit better, making for a more realistic effect.
azimuths = [45, 135, 225, 315]
altitudes = [45, 45, 45, 45]
multidirectional_hillshade = combine_hillshades(dem_path, azimuths, altitudes)
multidirectional_hillshade_image = rgb_image_from_array(multidirectional_hillshade)
multidirectional_hillshade_image.save('output/6-multidirectional_hillshade.png')
multidirectional_hillshade_blended = ImageChops.multiply(traditional_hillshade_blended, multidirectional_hillshade_image)
multidirectional_hillshade_blended.save('output/7-blended-multidirectional_hillshade.png')
Multiply multiplies the color values of the blend layer with the base layer, resulting in a darker overall image where the colors are combined. This will amplify the shadows a bit.
Low-light hillshade
The last hillshade is another traditional hillshade, but at a lower light angle. This one is blended using soft light.
low_light_hillshade = combine_hillshades(dem_path, [315], [25])
low_light_hillshade_image = rgb_image_from_array(low_light_hillshade)
low_light_hillshade_image.save('output/8-low_light_hillshade.png')
low_light_hillshade_blended = ImageChops.soft_light(multidirectional_hillshade_blended, low_light_hillshade_image)
low_light_hillshade_blended.save('output/9-low_light_hillshade_blended.png')
Lighting
The next layer is a lighting layer meant to brighten the peaks and darken the valleys. For this we can use an existing white-to-black (binary reversed) colormap from matplotlib.
lighting_dem = apply_colormap(normalized_dem, 'binary_r')
lighting_image = Image.fromarray(lighting_dem)
lighting_image.save('output/10-lighting.png')
lighting_blended = ImageChops.soft_light(low_light_hillshade_blended, lighting_image)
lighting_blended.save('output/11-lighting_blended.png')
Mist time!
The next two steps are my favorite. First, mist.
The mist layer requires another custom colormap. It needs to be semi-transparent white at low elevations, shifting quickly to fully transparent white before the elevation values get too high. Getting the transition just right took a bit of experimentation.
colors = [
(1.0, 1.0, 1.0, 0.85),
(1.0, 1.0, 1.0, 0.45),
(1.0, 1.0, 1.0, 0.15),
(1.0, 1.0, 1.0, 0.0),
(1.0, 1.0, 1.0, 0.0)
]
positions = [0.0, 0.15, 0.25, .35, 1.0]
mist_custom_cmap = LinearSegmentedColormap.from_list("custom_cmap", list(zip(positions, colors)))Each element in the colors list is a tuple containing values for RGBA (Red, Green, Blue, Alpha) color components. All colors have three ones, which is white. The only variation is in the alpha component, which determines transparency. The first record of 0.85 means 85% opaque (slight transparency).
The positions list specifies the exact points, ranging from 0 to 1, where each color from the colors list will be placed in the colormap. Each position value corresponds to a specific point in the gradient, determining where the transition between colors occurs. In this custom colormap, we gradually change the opacity of white from 85% at the start to full transparency by 35% of the way through, effectively making it fully transparent by the time we're one-third through the range from low to high elevations.
Let’s apply it.
mist_dem = mist_custom_cmap(normalized_dem)
mist_image = Image.fromarray((mist_dem * 255).astype(np.uint8), 'RGBA')
mist_image.save("output/12-dem_with_transparent_mist.png")
terrain_image = Image.open('output/11-lighting_blended.png').convert('RGBA')
combined_terrain_and_mist = Image.alpha_composite(terrain_image, mist_image).convert('RGB')
combined_terrain_and_mist.save('output/13-combined_terrain_and_mist.png')
And finally… slope as hillshade
For the final step, we generate a slope raster and apply the ‘cividis’ colormap to highlight steep areas.
def generate_slope_array(dem_path: str) -> Tuple[np.ndarray, Optional[float]]:
"""Generate a slope array from a DEM."""
dem_data = gdal.Open(dem_path, gdal.GA_ReadOnly)
slope_ds = gdal.DEMProcessing('', dem_data, 'slope', format='MEM', scale=1)
slope_band = slope_ds.GetRasterBand(1)
slope_array = slope_band.ReadAsArray()
no_data_value = slope_band.GetNoDataValue()
if no_data_value is not None:
mask = slope_array == no_data_value
slope_array = np.ma.masked_where(mask, slope_array)
min_val = np.min(slope_array)
max_val = np.max(slope_array)
norm_slope_array = (slope_array - min_val) / (max_val - min_val)
norm_slope_array = np.ma.filled(norm_slope_array, 0)
return norm_slope_array
cividis_cmap = plt.get_cmap('cividis')
slope_dem = cividis_cmap(generate_slope_array(dem_path))
slope_image = rgb_image_from_array(slope_dem)
slope_image.save('output/14-slope.png', 'png')
The final result
The final step is to blend the slope into the misty terrain with soft light.
slope_blended = ImageChops.soft_light(combined_terrain_and_mist, slope_image)
slope_blended.save('output/15-slope_blended.png')
I’m pretty happy with the results. I was able to get pretty close to the ArcGIS Pro version using just a few Python libraries.
Here are some close-ups.
Opportunities for improvement
I’d like to experiment with one hillshade step. Generating them incrementally allows you to see the effects and experiment, but using weights could get close to the end result in one step or fewer steps at least. I also feel like blending the standalone low-light hillshade undid some of the beauty of the multidirectional hillshade. Averaging with weights might help achieve a better balance of effects.
I’m new to image processing at this level, so I wasn’t always sure what channels/modes/normalization to use for the different effects. I’d like to review, refactor and understand that code better and see if there are opportunities to simplify.
Thanks to John Nelson for the original tutorial and the inspiration!
Complete code can be found on GitHub.