How to Create Detailed Terrain Maps in Python with Rasterio, GeoPandas, and Cartopy

This article introduces how to draw a terrain map using Python. Required third‑party modules include rasterio, geopandas, and cartopy, which can be installed via pip.

1 Example Code

1.1 Import Required Modules

import rasterio
import geopandas as gpd
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap

1.2 Set Map Font and Size

plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams["font.size"] = 14

1.3 Draw the Map

# Set projection: Mercator
projection = ccrs.Mercator()

# Create figure and axis
fig, ax = plt.subplots(figsize=(20, 10), subplot_kw={"projection": projection})

# Define map extent (coordinates in the chosen projection)
ax.set_extent([13530000, 14630000, 4960000, 5850000], crs=projection)

# Read vector file
shp = gpd.read_file("Data/Jilin_Mercator.shp")
shp.plot(ax=ax, transform=projection, edgecolor="black", linewidth=1, facecolor="none")

# Create custom colormap
colors = ["#369121", "#95C769", "#FFFFBF", "#E6865A", "#D14E30", "#BA1414"]
n_bins = 100
cmap_name = "green_brown"
cm = LinearSegmentedColormap.from_list(cmap_name, colors, N=n_bins)

# Read terrain raster data
dataset = rasterio.open("Data/DEM_Jilin_Mercator.tif")
data = dataset.read(1)  # first band
nodata_value = dataset.nodata
# Mask NoData values
data = np.ma.masked_where(data == nodata_value, data)

# Add raster to map
extent = [dataset.bounds.left, dataset.bounds.right, dataset.bounds.bottom, dataset.bounds.top]
im = ax.imshow(data, origin="upper", extent=extent, transform=projection, cmap=cm)

# Draw gridlines and labels
gl = ax.gridlines(draw_labels=True, linestyle="--", color="#4F4F4F")
gl.xlocator = plt.FixedLocator(range(120, 135, 3))
gl.ylocator = plt.FixedLocator(range(40, 50, 2))

# Add colorbar
cbar = plt.colorbar(im, ax=ax, orientation="horizontal", shrink=0.4, pad=0.1)
cbar.set_label("Elevation (m)", labelpad=10)
cbar.ax.xaxis.set_label_position("top")

# Save and show
plt.savefig("Pic.jpg", dpi=600)
plt.show()

2 Resulting Map

Link: https://www.cnblogs.com/qsgeo/p/18030978 (© Original author, removed if infringing)