Use this Skill for historical GIS: georeferencing historical maps with GDAL, digitizing boundaries, temporal territory change overlays, and historical population interpolation.
Use this Skill for historical GIS: georeferencing historical maps with GDAL, digitizing boundaries, temporal territory change overlays, and historical population interpolation.
Historical GIS: Georeferencing and Spatial Humanities
TL;DR — Georeference historical maps using GDAL, digitize territorial boundaries
into GeoDataFrames, compare multi-period territories with Hausdorff distance, interpolate
historical population from census points, and animate territorial change over centuries.
When to Use
Use this Skill when you need to:
Align a historical map scan to a modern coordinate reference system (CRS)
Digitize boundaries of historical polities, dioceses, or trade zones
Overlay territorial extents from multiple centuries in a single plot
Interpolate population density from sparse historical census data to polygon areas
Animate territorial change over time for scholarly visualization or publication
Do not use this Skill for:
Modern satellite or aerial image registration (use QGIS auto-registration tools)
Vector map editing at scale (use QGIS or ArcGIS Pro workflows)
Background
Georeferencing assigns real-world coordinates to a raster image (map scan) using
Ground Control Points (GCPs) — identifiable features on the historical map matched
to known modern coordinates (church spires, river confluences, coast outlines).
1st order (affine, 6 params), 2nd order (quadratic, 12), 3rd order (20)
RMS error
Root mean square residual of GCP reprojection; <2 px target for 1st order
EPSG:4326
WGS84 geographic CRS; EPSG:3857 = Web Mercator for display
Hausdorff distance
Max of directed distances between two boundary curves
Voronoi / Thiessen
Partition space so each polygon contains all points nearer to its site
GDAL gdal_translate embeds GCPs into the raster metadata; gdalwarp then resamples
to produce a properly georeferenced GeoTIFF that can be loaded in any GIS tool.
import subprocess
import os
from pathlib import Path
defwrite_gcp_file(gcps: list[dict], gcp_path: str) -> None:
"""
Write a GCP text file compatible with GDAL gdal_translate -gcp flag.
Each GCP dict must contain:
pixel_x (float): column in the source image
pixel_y (float): row in the source image
geo_x (float): target longitude (EPSG:4326)
geo_y (float): target latitude (EPSG:4326)
label (str): human-readable identifier
Args:
gcps: List of GCP dicts.
gcp_path: Output path for the generated shell-script gcp file.
"""
lines = ["#!/usr/bin/env bash", "# Auto-generated GCPs for gdal_translate", ""]
for g in gcps:
lines.append(
f"-gcp {g['pixel_x']}{g['pixel_y']}{g['geo_x']}{g['geo_y']}"f" # {g.get('label', '')}"
)
Path(gcp_path).write_text("\n".join(lines), encoding="utf-8")
print(f"GCP file written to {gcp_path}")
defgeoreference_map(
input_raster: str,
output_raster: str,
gcps: list[dict],
order: int = 1,
target_crs: str = "EPSG:4326",
resampling: str = "bilinear",
) -> dict:
"""
Georeference a historical map scan using GDAL GCPs + polynomial warp.
Steps:
1. gdal_translate: embed GCPs into intermediate GeoTIFF.
2. gdalwarp: polynomial transform and resample to target CRS.
Args:
input_raster: Path to the raw scan (JPG, PNG, TIFF).
output_raster: Path to write the georeferenced GeoTIFF.
gcps: List of GCP dicts (pixel_x, pixel_y, geo_x, geo_y, label).
order: Polynomial order (1=affine, 2=quadratic, 3=cubic).
target_crs: EPSG code string for the output CRS.
resampling: Resampling algorithm: "bilinear" | "cubic" | "near".
Returns:
Dict with intermediate_path, output_path, rms_estimate.
"""
intermediate = input_raster.replace(".", "_gcps.")
ifnot intermediate.endswith(".tif"):
intermediate += ".tif"# Build -gcp arguments
gcp_args = []
for g in gcps:
gcp_args.extend([
"-gcp",
str(g["pixel_x"]), str(g["pixel_y"]),
str(g["geo_x"]), str(g["geo_y"]),
])
# Step 1: embed GCPs
cmd_translate = (
["gdal_translate", "-of", "GTiff"]
+ gcp_args
+ [input_raster, intermediate]
)
result = subprocess.run(cmd_translate, capture_output=True, text=True)
if result.returncode != 0:
raise RuntimeError(f"gdal_translate failed:\n{result.stderr}")
# Step 2: warp to target CRS
cmd_warp = [
"gdalwarp",
"-order", str(order),
"-r", resampling,
"-t_srs", target_crs,
"-overwrite",
intermediate, output_raster,
]
result = subprocess.run(cmd_warp, capture_output=True, text=True)
if result.returncode != 0:
raise RuntimeError(f"gdalwarp failed:\n{result.stderr}")
# Rough RMS estimate from residuals (manual if needed)print(f"Georeferenced raster written to {output_raster}")
return {
"intermediate_path": intermediate,
"output_path": output_raster,
"gcp_count": len(gcps),
"order": order,
"target_crs": target_crs,
}
Step 2 — Multi-Period Territory Overlay Plot
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from shapely.geometry import Polygon
import numpy as np
defcreate_historical_territory(
name: str,
coords: list[tuple[float, float]],
period: str,
color: str,
) -> dict:
"""
Construct a historical territory record from a polygon coordinate list.
Args:
name: Territory name (e.g. "Holy Roman Empire").
coords: List of (lon, lat) tuples defining the boundary.
period: Century or date label (e.g. "1250 CE").
color: Matplotlib color string for this territory.
Returns:
Dict with keys: name, period, geometry, color.
"""
poly = Polygon(coords)
return {"name": name, "period": period, "geometry": poly, "color": color}
defplot_temporal_territories(
territories: list[dict],
title: str = "Historical Territory Change",
output_path: str = None,
alpha: float = 0.35,
) -> None:
"""
Overlay historical territories from multiple periods on a single map.
Args:
territories: List of territory dicts from create_historical_territory().
title: Plot title string.
output_path: If given, save figure here.
alpha: Polygon fill transparency.
"""
gdf = gpd.GeoDataFrame(territories, crs="EPSG:4326")
fig, ax = plt.subplots(figsize=(12, 8))
periods = gdf["period"].unique()
legend_patches = []
for territory in territories:
geom = territory["geometry"]
x, y = geom.exterior.xy
ax.fill(x, y, alpha=alpha, color=territory["color"], linewidth=1.5,
edgecolor=territory["color"])
# Build legend: one entry per territory name+periodfor _, row in gdf.iterrows():
patch = mpatches.Patch(
color=row["color"], alpha=0.7,
label=f"{row['name']} ({row['period']})"
)
legend_patches.append(patch)
ax.legend(handles=legend_patches, loc="lower left", fontsize=8)
ax.set_title(title, fontsize=14, fontweight="bold")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.grid(True, alpha=0.3)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=150, bbox_inches="tight")
print(f"Territory map saved to {output_path}")
plt.show()
defhausdorff_distance_territories(
geom_a: Polygon,
geom_b: Polygon,
) -> float:
"""
Compute the Hausdorff distance between two historical boundary polygons.
The Hausdorff distance is the maximum distance from any point on boundary A
to the nearest point on boundary B (symmetric max). Values are in CRS units
(degrees for EPSG:4326; project to a metric CRS for meaningful km values).
Args:
geom_a: Shapely Polygon for period A boundary.
geom_b: Shapely Polygon for period B boundary.
Returns:
Hausdorff distance as float in CRS units.
"""return geom_a.hausdorff_distance(geom_b)
Step 3 — Historical Population Interpolation Between Census Years
from scipy.interpolate import griddata
import rasterio
from rasterio.transform import from_bounds
definterpolate_historical_population(
census_points: list[dict],
bbox: tuple[float, float, float, float],
output_raster: str,
resolution: float = 0.1,
method: str = "linear",
) -> np.ndarray:
"""
Interpolate historical census point data to a population density raster.
Args:
census_points: List of dicts with keys: lon, lat, population, year.
bbox: Bounding box (min_lon, min_lat, max_lon, max_lat).
output_raster: Path to write the output GeoTIFF.
resolution: Grid cell size in CRS units (degrees for WGS84).
method: Scipy interpolation method: "linear" | "cubic" | "nearest".
Returns:
Interpolated grid as 2D numpy array.
"""
min_lon, min_lat, max_lon, max_lat = bbox
lon_vals = np.array([p["lon"] for p in census_points])
lat_vals = np.array([p["lat"] for p in census_points])
pop_vals = np.array([p["population"] for p in census_points], dtype=float)
# Normalise population to density per sq degree
grid_lon = np.arange(min_lon, max_lon, resolution)
grid_lat = np.arange(min_lat, max_lat, resolution)
grid_x, grid_y = np.meshgrid(grid_lon, grid_lat)
grid_pop = griddata(
points=np.column_stack([lon_vals, lat_vals]),
values=pop_vals,
xi=(grid_x, grid_y),
method=method,
fill_value=0.0,
)
grid_pop = np.nan_to_num(grid_pop, nan=0.0)
# Write as single-band GeoTIFF
transform = from_bounds(min_lon, min_lat, max_lon, max_lat,
grid_pop.shape[1], grid_pop.shape[0])
with rasterio.open(
output_raster, "w",
driver="GTiff",
height=grid_pop.shape[0],
width=grid_pop.shape[1],
count=1,
dtype=grid_pop.dtype,
crs="EPSG:4326",
transform=transform,
) as dst:
dst.write(grid_pop, 1)
print(f"Population raster saved to {output_raster}")
return grid_pop
Advanced Usage
Territory Change Animation
import matplotlib.animation as animation
defanimate_territory_change(
territory_frames: list[list[dict]],
period_labels: list[str],
output_gif: str,
fps: int = 2,
) -> None:
"""
Create an animated GIF showing territory change across periods.
Args:
territory_frames: List of territory lists, one per time step.
period_labels: Labels for each time step (same length as territory_frames).
output_gif: Path to write the output GIF.
fps: Frames per second.
"""
fig, ax = plt.subplots(figsize=(10, 7))
defupdate(frame_idx: int):
ax.clear()
territories = territory_frames[frame_idx]
for t in territories:
geom = t["geometry"]
x, y = geom.exterior.xy
ax.fill(x, y, alpha=0.4, color=t["color"], edgecolor=t["color"], linewidth=1.5)
ax.set_title(f"Territory — {period_labels[frame_idx]}", fontsize=14)
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.grid(True, alpha=0.3)
ani = animation.FuncAnimation(
fig, update, frames=len(territory_frames),
interval=1000 // fps, repeat=True,
)
ani.save(output_gif, writer="pillow", fps=fps)
print(f"Animation saved to {output_gif}")
plt.close(fig)
Troubleshooting
Problem
Cause
Fix
gdalwarp: command not found
GDAL not on PATH
export PATH=/usr/bin:$PATH; verify with which gdalwarp
High RMS error (>5 px) after 1st order
Non-linear map projection or distortion
Use order=2 or add more GCPs near edges
CRSError in rasterio
CRS string not recognised
Use EPSG integer: crs=4326; verify PROJ data dir
Shapely TopologicalError on Polygon
Self-intersecting coordinates
Apply .buffer(0) to fix topology
Population grid shows NaN islands
Insufficient census points for interpolation area
Switch to method="nearest" or add synthetic boundary points
