Use this Skill for historical GIS: georeferencing historical maps with GDAL, digitizing boundaries, temporal territory change overlays, and historical population interpolation.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
A direct command skips the review prompt. Inspect the source before running it.
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
