Automating Multi-Layer Legend Creation with GeoPandas

Call ax.get_legend_handles_labels() once after all layers are drawn, deduplicate by label string while preserving insertion order, remove Matplotlib’s auto-generated legend fragments, then rebuild a single ax.legend() — this gives you a deterministic, unified legend across polygon, line, and point layers in under 20 lines of Python.

Core Algorithm and Workflow

GeoPandas delegates rendering to Matplotlib’s object-oriented API. Each gdf.plot(legend=True) call appends new Patch, Line2D, or PathCollection artist handles to the active Axes object — but it does not merge them with handles from previous calls. The result without intervention: stacked legend boxes, duplicated categorical labels, and missing symbology for continuous ramps.

The aggregation algorithm has five deterministic steps:

1. Shared Axes initialization. Create one fig, ax = plt.subplots() pair. Pass ax=ax to every .plot() call so all layers share the same coordinate frame and handle registry.
2. Handle registration. Pass legend=True to each .plot() call. This forces GeoPandas to generate proxy artists for each unique category value, registering them in ax’s internal _get_lines or legend_handles state.
3. Single-pass extraction. After the final .plot(), call raw_handles, raw_labels = ax.get_legend_handles_labels(). This retrieves every registered handle — across all layers — as flat lists.
4. Order-preserving deduplication. Iterate zip(raw_handles, raw_labels) with a seen: set[str] tracker. Append to unique_pairs on the first occurrence of each label; skip repeats. This runs in O(n) time and preserves the draw order, which matches layer stacking.
5. Legend reconstruction. Remove Matplotlib’s auto-generated legend artifact with ax.get_legend().remove(), then call ax.legend(unique_handles, unique_labels, **kwargs) with explicit positioning, title, and font size.

The diagram below illustrates the data flow from individual layer renders to the unified legend object:

Production-Ready Python Implementation

The script below creates synthetic polygon, line, and point layers, plots them on a shared Axes, and assembles a single unified legend. Run it headlessly under the Agg backend for CI/CD compatibility.

import matplotlib
matplotlib.use("Agg")  # must precede pyplot import on headless runners

import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import LineString, Point, Polygon

# ---------------------------------------------------------------------------
# 1. Synthetic multi-layer data (replace with real GeoDataFrames in production)
# ---------------------------------------------------------------------------
polygons = gpd.GeoDataFrame(
    {"zone": ["Zone A", "Zone B", "Zone C"]},
    geometry=[
        Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]),
        Polygon([(2, 1), (4, 1), (4, 3), (2, 3)]),
        Polygon([(1, 2), (3, 2), (3, 4), (1, 4)]),
    ],
    crs="EPSG:4326",
)

lines = gpd.GeoDataFrame(
    {"route": ["Highway 1", "Highway 2"]},
    geometry=[
        LineString([(0.5, 0.5), (3.5, 2.5)]),
        LineString([(1.5, 3.5), (3.5, 1.5)]),
    ],
    crs="EPSG:4326",
)

points = gpd.GeoDataFrame(
    {"facility": ["Station X", "Station Y", "Station Z"]},
    geometry=[Point(1, 1), Point(3, 2), Point(2, 3)],
    crs="EPSG:4326",
)

# ---------------------------------------------------------------------------
# 2. Shared Axes — every layer shares this coordinate frame
# ---------------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(10, 8))
plt.rcParams["font.family"] = "DejaVu Sans"  # stable across Linux CI runners

# ---------------------------------------------------------------------------
# 3. Plot layers; legend=True registers handles with the shared Axes
# ---------------------------------------------------------------------------
polygons.plot(
    ax=ax, column="zone", cmap="Pastel1", legend=True,
    edgecolor="black", linewidth=1.2,
)
lines.plot(ax=ax, color="darkred", linewidth=2.5, legend=True, label="Highway")
points.plot(ax=ax, color="navy", marker="s", markersize=80, legend=True, label="Station")

# ---------------------------------------------------------------------------
# 4. Extract all registered handles and labels in a single pass
# ---------------------------------------------------------------------------
raw_handles, raw_labels = ax.get_legend_handles_labels()

# Preserve draw order; drop the first repeated occurrence of each label string
seen: set[str] = set()
unique_handles, unique_labels = [], []
for handle, label in zip(raw_handles, raw_labels):
    if label not in seen:
        seen.add(label)
        unique_handles.append(handle)
        unique_labels.append(label)

# ---------------------------------------------------------------------------
# 5. Remove Matplotlib's auto-generated legend artifact and rebuild
# ---------------------------------------------------------------------------
existing_legend = ax.get_legend()
if existing_legend is not None:
    existing_legend.remove()

ax.legend(
    unique_handles,
    unique_labels,
    loc="upper right",
    frameon=True,
    framealpha=0.95,
    fontsize=10,
    title="Map Features",
    title_fontsize=11,
    fancybox=False,   # flat corners aid WCAG contrast compliance
    shadow=False,
)

# ---------------------------------------------------------------------------
# 6. Finalize and export
# ---------------------------------------------------------------------------
ax.set_axis_off()
ax.set_title("Multi-Layer Geospatial Map", fontsize=14, pad=15)
fig.savefig("multi_layer_legend.png", dpi=300, bbox_inches="tight", facecolor="white")
plt.close(fig)

Handle type reference

ax.legend() accepts all three discrete artist types in a single call. Continuous colormaps produce a Colorbar rather than a legend handle — retrieve the ScalarMappable from ax.collections[-1] and call fig.colorbar() with location='right' so it sits outside the categorical legend area.

Performance Tuning and Cartographic Best Practices

• O(n) deduplication scales cleanly. The seen set lookup is O(1); iterating raw_handles once is O(n) where n is the total number of registered handles. For 10+ layers with 20+ categories each, this remains sub-millisecond. Avoid sorting handles inside the deduplication loop — sort once after the unique list is built if you need alphabetical label order.

• Prefix labels when symbology differs between layers. If two layers share a label string (e.g., both produce "Urban"), the deduplication step drops one handle and loses its symbol. Assign layer-scoped labels before plotting: gdf["_legend_label"] = "Urban (fill)" for polygons and "Urban (outline)" for lines, using a temporary column rather than mutating the original data.

• Cache figure state for batch map runs. When generating 50+ maps in a loop, call plt.rcParams.update({...}) once before the loop rather than inside it, and reuse the same fig, ax pair with ax.cla() between iterations. This avoids the overhead of plt.subplots() for each output file, which becomes significant at scale.

• Set bbox_inches="tight" consistently. Without it, the exported PNG may clip the legend when it extends outside the default figure boundary. Pair this with an explicit facecolor="white" to prevent transparent backgrounds in CMYK workflows that expect opaque whites — relevant if the maps feed into print-ready export pipelines.

• Validate legend-data parity before export. After rebuilding the legend, assert len(unique_labels) == len(expected_categories) where expected_categories is derived from pd.concat([gdf[col] for gdf, col in layer_column_pairs]).unique(). A mismatch indicates a silent deduplication collision and prevents exporting corrupt cartographic output.

Integration and Next Steps

This legend aggregation pattern slots directly into automated map generation pipelines. The reconstructed ax.legend() object is export-ready for both raster (fig.savefig(..., dpi=300)) and vector (format="svg") outputs, which matters when the downstream step involves high-resolution vector export for print or atlas production.

In GitHub Actions or Jenkins pipelines, pair this script with matplotlib.use("Agg") at the top of the module and a MPLCONFIGDIR environment variable pointed at a writable temp path to prevent font cache errors on read-only runners. Pass plt.rcParams["font.family"] = "DejaVu Sans" explicitly — this font ships with Matplotlib and eliminates fallback mismatches across Ubuntu 22.04 and 24.04 LTS runner images.

For maps that also require label collision avoidance, run the collision detection step after legend reconstruction. The legend bounding box changes the available canvas area for labels, and STRtree-based placement engines need the final legend geometry to compute valid label positions.

When applying rule-based styling engines — for instance, JSON-driven style rules that assign colors and markers per feature class — integrate the label-prefix strategy described above so the rule engine controls both the symbology and the human-readable legend label in a single configuration object. This eliminates the manual mapping step between render-time color values and legend text.

Related

• Dynamic Legend Generation — parent cluster covering handle extraction, continuous colormap proxies, and legend-data parity validation across the full automated export workflow.
• Solving Label Overlap in Dense Urban Maps with Python — STRtree-based spatial indexing to place map labels without collision, relevant after legend geometry is finalized.
• Building a JSON-Based Rule Styling Engine for QGIS — declarative styling rules that can drive both symbology and legend labels from a single configuration source.

Back to Programmatic Map Styling and Label Automation