How to Plot Multiple Layers on One Map in GeoPandas

Problem statement

Three layers, three calls, three separate figures:

import matplotlib
matplotlib.use("Agg")
import geopandas as gpd
import matplotlib.pyplot as plt

boundary, roads, sites = (gpd.read_file(f"{name}.gpkg") for name in ("boundary", "roads", "sites"))
boundary.plot()
roads.plot()
sites.plot()

Adding ax=ax gets them onto one figure and introduces the next set of problems:

fig, ax = plt.subplots()
boundary.plot(ax=ax)
roads.plot(ax=ax)
sites.plot(ax=ax)

The boundary is a solid blue mass that hides everything drawn after it. The roads are invisible against it. The sites are there, somewhere. And here, with the ward in British National Grid metres and the roads and bus stops in degrees, the ward sits at x ≈ 430,000 and everything else at x ≈ −1.5 — matplotlib helpfully zooms out to fit both, leaving specks at opposite corners.

A multi-layer map needs four things settled: a shared CRS, a deliberate drawing order, fill and stroke that let layers coexist, and an extent chosen rather than inherited.

Quick answer

import geopandas as gpd
import matplotlib.pyplot as plt

CRS = 27700
boundary = gpd.read_file("boundary.gpkg").to_crs(CRS)
roads    = gpd.read_file("roads.gpkg").to_crs(CRS)
sites    = gpd.read_file("sites.gpkg").to_crs(CRS)

fig, ax = plt.subplots(figsize=(10, 10))

boundary.plot(ax=ax, facecolor="#f1f5f9", edgecolor="#94a3b8",
              linewidth=0.8, zorder=1)
roads.plot(ax=ax, color="#64748b", linewidth=0.6, zorder=2)
sites.plot(ax=ax, color="#ef4444", markersize=45,
           edgecolor="white", linewidth=0.8, zorder=3)

ax.set_xlim(*boundary.total_bounds[[0, 2]])
ax.set_ylim(*boundary.total_bounds[[1, 3]])
ax.set_axis_off()
fig.savefig("quick.png", dpi=80)

One rule prevents most of it: every layer that is context rather than content gets facecolor="none" or a low alpha.

boundary   EPSG:27700     [428641.298, 428899.398, 434245.944, 433503.8]
roads      EPSG:4326      [-1.589, 53.749, -1.472, 53.803]
sites      EPSG:4326      [-1.558, 53.762, -1.496, 53.797]

One CRS, bottom-up order, extent, legend

1. Put every layer in the same CRS

Matplotlib draws whatever coordinates it is given. Two layers in different CRS are two sets of numbers on one pair of axes, and nothing warns you:

layers = {name: gpd.read_file(f"{name}.gpkg") for name in ("boundary", "roads", "sites")}   # as delivered
for name, layer in layers.items():

    print(f"{name:<10} {str(layer.crs):<14} {layer.total_bounds.round(3).tolist()}")

boundary   EPSG:27700     [428641.298, 428899.398, 434245.944, 433503.8]
roads      EPSG:4326      [-1.589, 53.749, -1.472, 53.803]
sites      EPSG:4326      [-1.558, 53.762, -1.496, 53.797]

The mismatch is unmistakable once printed. Reproject everything to one CRS up front:

CRS = 27700
layers = {name: layer.to_crs(CRS) for name, layer in layers.items()}

Choose that CRS deliberately — see choosing a map projection for display. If one layer has crs = None, to_crs raises, and guessing is not a fix; see raster and vector do not line up for the diagnosis.

2. Draw from the bottom up

Matplotlib draws in call order, so a later layer covers an earlier one. That works until you reorder the calls and the map changes. zorder makes it explicit and order-independent:

boundary.plot(ax=ax, zorder=1, ...)
roads.plot(ax=ax, zorder=3, ...)
sites.plot(ax=ax, zorder=4, ...)

Leave gaps between the numbers. Inserting a layer between 1 and 3 needs no renumbering; inserting between 1 and 2 does.

The general principle is large and contextual at the bottom, small and important at the top, because a small feature drawn under a large one is simply gone.

3. Make polygons stop hiding each other

The default plot() gives a solid fill, which is right for exactly one layer per map:

# outline only — context that hides nothing
boundary.plot(ax=ax, facecolor="none", edgecolor="#334155", linewidth=1.2, zorder=5)

# semi-transparent fill — a thematic layer over a basemap
zones.plot(ax=ax, column="risk", cmap="YlOrRd", alpha=0.55, zorder=2)

# hatched — for qualitative overlays such as exclusion areas
excluded.plot(ax=ax, facecolor="none", edgecolor="#ef4444",
              hatch="///", linewidth=0.8, zorder=4)

# just the boundary of a polygon layer, as lines
boundary.boundary.plot(ax=ax, color="#334155", linewidth=1.2, zorder=5)

.boundary returns a GeoSeries of the polygon outlines, which is subtly different from facecolor="none": it is a line layer, so it takes line styling and never fills.

Note that alpha applies to the whole layer including its edges, so a semi-transparent layer has semi-transparent outlines. When you want a solid outline over a translucent fill, plot the layer twice:

zones.plot(ax=ax, column="risk", cmap="YlOrRd", alpha=0.55, zorder=2)
zones.boundary.plot(ax=ax, color="#334155", linewidth=0.4, zorder=3)

4. Set the extent from the layer that defines it

By default matplotlib fits everything drawn, which means one stray feature — a site in the wrong CRS, a national boundary behind a city map — sets the scale:

minx, miny, maxx, maxy = boundary.total_bounds
pad = 0.05 * max(maxx - minx, maxy - miny)
ax.set_xlim(minx - pad, maxx + pad)
ax.set_ylim(miny - pad, maxy + pad)

Set the limits after all layers are drawn, or a later plot() call will autoscale over them.

Clipping the data rather than the view is the other option, and it is better when the extra features would still influence a legend or a classification:

roads_in = gpd.clip(roads, boundary)          # only what is inside, cut at the edge
touching = roads[roads.intersects(boundary.union_all())]
print(f"{len(roads):,} roads in the file; {len(touching):,} touch the ward; clipped length "
      f"{roads_in.length.sum() / 1000:,.1f} km against {touching.length.sum() / 1000:,.1f} km uncut")

9,802 roads in the file; 4,522 touch the ward; clipped length 306.8 km against 324.7 km uncut

gpd.clip cuts geometries at the boundary; roads[roads.intersects(area)] keeps whole features that touch it. The first gives a clean edge, the second keeps features intact — see how to clip spatial data.

5. Build a legend that describes the layers

GeoPandas' automatic legend describes one column of one layer. A multi-layer map needs one built by hand:

from matplotlib.lines import Line2D
from matplotlib.patches import Patch

handles = [
    Patch(facecolor="#f1f5f9", edgecolor="#94a3b8", label="Ward boundary"),
    Line2D([0], [0], color="#64748b", linewidth=1.5, label="Roads"),
    Line2D([0], [0], marker="o", color="none", markerfacecolor="#ef4444",
           markeredgecolor="white", markersize=9, label="Bus stops"),
]
ax.legend(handles=handles, loc="lower right", frameon=True,
          framealpha=0.9, fontsize=9)

Patch for areas, Line2D for lines, and Line2D with color="none" plus a marker for points. Matching each handle's styling to the layer it describes is the whole job — a legend whose swatches do not match the map is worse than none.

Layered maps of Hunslet & Riverside

The data used on this page

One Leeds electoral ward, Hunslet & Riverside — the city centre south of the station and the riverside to the east — with the layers a map of it would carry. Each file is written in the CRS its source delivers, which is what step 1 is about:

• boundary.gpkg — the ward, from the ONS December 2024 boundaries, in British National Grid (EPSG:27700).

• roads.gpkg and parks.gpkg — roads in the ward's bounding box, and park polygons clipped to the ward, from Geofabrik's 1 January 2026 West Yorkshire OpenStreetMap extract, in WGS84.

• sites.gpkg — the NaPTAN bus stops in the ward, in WGS84.

• incidents.gpkg — police.uk street-level crime for June 2026 in the ward, in WGS84.

import io

import geopandas as gpd
import numpy as np
import pandas as pd
import requests

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
boundary = gpd.read_file(f"{ONS}/Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
                         "?where=WD24NM%3D%27Hunslet%20%26%20Riverside%27&outFields=WD24NM&outSR=27700&f=geojson")
boundary = boundary.set_crs(27700, allow_override=True)
boundary.to_file("boundary.gpkg")
ward = boundary.to_crs(4326)
west, south, east, north = ward.total_bounds

OSM = "https://download.geofabrik.de/europe/united-kingdom/england/west-yorkshire-260101-free.shp.zip"
with open("west-yorkshire.zip", "wb") as f:
    f.write(requests.get(OSM, timeout=600).content)
roads = gpd.read_file("zip://west-yorkshire.zip!gis_osm_roads_free_1.shp", bbox=(west, south, east, north),
                      columns=["osm_id", "fclass", "name"])

roads.to_file("roads.gpkg")                        # everything in the ward's bounding box, uncut
landuse = gpd.read_file("zip://west-yorkshire.zip!gis_osm_landuse_a_free_1.shp", bbox=(west, south, east, north))
gpd.clip(landuse[landuse["fclass"] == "park"], ward).to_file("parks.gpkg")

# NaPTAN area 450 as saved on 9 October 2026 (the register changes weekly); live:
# naptan.api.dft.gov.uk/v1/access-nodes?dataFormat=csv&atcoAreaCodes=450
naptan = pd.read_csv(io.StringIO(requests.get("https://data.spatialworkflow.io/naptan-area-450-2026-10-09.csv", timeout=300).text),
                     low_memory=False)
naptan = naptan[(naptan["StopType"] == "BCT") & (naptan["Status"] == "active")]
stops = gpd.GeoDataFrame(naptan[["ATCOCode", "CommonName"]], crs=4326,
                         geometry=gpd.points_from_xy(naptan["Longitude"], naptan["Latitude"]))
stops[stops.within(ward.union_all())].to_file("sites.gpkg")

crimes = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120, data={
    "poly": f"{south},{west}:{south},{east}:{north},{east}:{north},{west}", "date": "2026-06"}).json()
crimes = gpd.GeoDataFrame({"category": [c["category"] for c in crimes]}, crs=4326,
                          geometry=gpd.points_from_xy([float(c["location"]["longitude"]) for c in crimes],
                                                      [float(c["location"]["latitude"]) for c in crimes]))
crimes[crimes.within(ward.union_all())].to_file("incidents.gpkg")

for name in ("boundary", "roads", "parks", "sites", "incidents"):
    layer = gpd.read_file(f"{name}.gpkg")
    print(f"{name:<10} {len(layer):>6,} features  {layer.crs}")

boundary        1 features  EPSG:27700
roads       9,802 features  EPSG:4326
parks          31 features  EPSG:4326
sites         126 features  EPSG:4326
incidents     589 features  EPSG:4326

Example 1: a layered map function

import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib.patches import Patch

DEFAULTS = {
    "Polygon":    {"facecolor": "none", "edgecolor": "#475569", "linewidth": 0.8},
    "LineString": {"color": "#64748b", "linewidth": 0.7},
    "Point":      {"color": "#ef4444", "markersize": 35,
                   "edgecolor": "white", "linewidth": 0.7},
}

def layered_map(spec, *, crs=None, extent_from=None, figsize=(11, 11),
                title=None, legend=True):
    """spec: ordered list of (label, gdf, style_overrides). Bottom layer first."""
    crs = crs or next(g.crs for _, g, _ in spec if g.crs is not None)
    prepared = []
    for i, (label, gdf, style) in enumerate(spec):
        if gdf.crs is None:
            raise ValueError(f"layer '{label}' has no CRS")
        prepared.append((label, gdf.to_crs(crs), style, (i + 1) * 2))

    fig, ax = plt.subplots(figsize=figsize)
    handles = []
    for label, gdf, style, z in prepared:
        kind = gdf.geom_type.iloc[0].replace("Multi", "")
        kwargs = {**DEFAULTS.get(kind, {}), **style, "zorder": z}
        gdf.plot(ax=ax, **kwargs)

        if label:
            if kind == "Point":
                handles.append(Line2D([0], [0], marker="o", color="none",
                                      markerfacecolor=kwargs.get("color", "#ef4444"),
                                      markeredgecolor=kwargs.get("edgecolor", "white"),
                                      markersize=8, label=label))
            elif kind == "LineString":
                handles.append(Line2D([0], [0], color=kwargs.get("color", "#64748b"),
                                      linewidth=max(kwargs.get("linewidth", 1), 1.4),
                                      label=label))
            else:
                handles.append(Patch(facecolor=kwargs.get("facecolor", "none"),
                                     edgecolor=kwargs.get("edgecolor", "#475569"),
                                     hatch=kwargs.get("hatch"), label=label))

    by_label = {label: gdf for label, gdf, _, _ in prepared}
    base = by_label.get(extent_from, prepared[0][1]) if extent_from else prepared[0][1]
    minx, miny, maxx, maxy = base.total_bounds
    pad = 0.05 * max(maxx - minx, maxy - miny)
    ax.set_xlim(minx - pad, maxx + pad)
    ax.set_ylim(miny - pad, maxy + pad)

    if legend and handles:
        ax.legend(handles=handles, loc="lower right", frameon=True,
                  framealpha=0.9, fontsize=9)
    if title:
        ax.set_title(title, fontsize=14, loc="left")
    ax.set_axis_off()
    return fig, ax

fig, ax = layered_map([
    ("Ward",        gpd.read_file("boundary.gpkg"), {"facecolor": "#f8fafc"}),
    ("Green space", gpd.read_file("parks.gpkg"),    {"facecolor": "#dcfce7", "edgecolor": "#86efac"}),
    ("Roads",       gpd.read_file("roads.gpkg"),    {"linewidth": 0.5}),
    ("Bus stops",   gpd.read_file("sites.gpkg"),    {"markersize": 12}),
], crs=27700, extent_from="Ward", title="Hunslet & Riverside: bus stops")
fig.savefig("layered.png", dpi=80)
print([h.get_label() for h in ax.get_legend().legend_handles])

['Ward', 'Green space', 'Roads', 'Bus stops']

Three properties make this usable. The CRS is taken from the first layer that has one and applied to all, so a mismatch cannot survive. zorder comes from list position, so the order is visible in the call rather than in scattered arguments. And the legend handle is built from the same kwargs dict that styled the layer, so the swatch cannot drift out of sync with the map — the failure that makes hand-built legends untrustworthy.

Raising on a CRS-less layer is deliberate: silently assuming it matches is how a layer ends up drawn in the wrong hemisphere.

Example 2: layer-specific styling from an attribute

Real maps style within a layer as well as between layers:

import geopandas as gpd
import matplotlib.pyplot as plt

roads = gpd.read_file("roads.gpkg").to_crs(27700)          # Geofabrik's classes are in `fclass`

ROAD_STYLE = {
    "motorway":   {"color": "#1e3a8a", "linewidth": 2.4, "zorder": 6},
    "primary":    {"color": "#2563eb", "linewidth": 1.6, "zorder": 5},
    "secondary":  {"color": "#60a5fa", "linewidth": 1.0, "zorder": 4},
    "residential": {"color": "#cbd5e1", "linewidth": 0.5, "zorder": 3},
}

fig, ax = plt.subplots(figsize=(12, 12))
gpd.read_file("boundary.gpkg").to_crs(27700).plot(
    ax=ax, facecolor="#f8fafc", edgecolor="#e2e8f0", zorder=1)

drawn = 0
for kind, style in ROAD_STYLE.items():
    subset = roads[roads["fclass"] == kind]
    if subset.empty:
        print(f"  · no {kind} roads in this extent")
        continue
    subset.plot(ax=ax, **style)
    drawn += len(subset)
    print(f"  ✓ {kind:<12} {len(subset):>6,} features")

unstyled = roads[~roads["fclass"].isin(ROAD_STYLE)]
if len(unstyled):
    print(f"  ! {len(unstyled):,} features with an unhandled class: "
          f"{unstyled['fclass'].value_counts().head(5).to_dict()}")
ax.set_axis_off()

  ✓ motorway         51 features
  ✓ primary         307 features
  ✓ secondary        48 features
  ✓ residential   1,228 features
  ! 8,168 features with an unhandled class: {'service': 2863, 'footway': 2278, 'cycleway': 792, 'path': 435, 'tertiary': 420}

The unstyled report is the part worth copying. Styling by dictionary lookup silently drops anything not in the dictionary — here most of the features in the file, footpaths and service roads among them — which is the sort of thing that goes unnoticed until someone asks why a road is absent. Explicit ordering by zorder also means motorways draw over residential streets regardless of loop order.

Example 3: small multiples sharing one style

import geopandas as gpd
import matplotlib.pyplot as plt
import math

def small_multiples(gdf, group_col, context=None, *, crs=27700, ncols=4,
                    figsize_per=(3.2, 3.2), column=None, **plot_kwargs):
    gdf = gdf.to_crs(crs)
    context = context.to_crs(crs) if context is not None else None
    groups = sorted(gdf[group_col].dropna().unique())
    nrows = math.ceil(len(groups) / ncols)

    fig, axes = plt.subplots(nrows, ncols,
                             figsize=(figsize_per[0] * ncols, figsize_per[1] * nrows))
    axes = axes.ravel() if len(groups) > 1 else [axes]

    # one shared colour range so panels are comparable
    if column:
        plot_kwargs.setdefault("vmin", gdf[column].min())
        plot_kwargs.setdefault("vmax", gdf[column].max())

    for ax, key in zip(axes, groups):
        part = gdf[gdf[group_col] == key]
        if context is not None:
            context.plot(ax=ax, facecolor="#f1f5f9", edgecolor="#e2e8f0", zorder=1)
        part.plot(ax=ax, column=column, zorder=2, **plot_kwargs)
        ax.set_title(f"{key}  (n={len(part):,})", fontsize=9)
        ax.set_axis_off()
        if context is not None:
            minx, miny, maxx, maxy = context.total_bounds
            ax.set_xlim(minx, maxx); ax.set_ylim(miny, maxy)

    for ax in axes[len(groups):]:
        ax.set_axis_off()
    plt.tight_layout()
    return fig

incidents = gpd.read_file("incidents.gpkg")
top = incidents["category"].value_counts().head(8).index
fig = small_multiples(incidents[incidents["category"].isin(top)], "category",
                      context=gpd.read_file("boundary.gpkg"),
                      color="#ef4444", markersize=6, alpha=0.6)
fig.savefig("multiples.png", dpi=70)
print(incidents["category"].value_counts().head(8).to_string())

category
violent-crime            200
shoplifting               77
anti-social-behaviour     61
public-order              47
other-theft               37
drugs                     36
criminal-damage-arson     33
vehicle-crime             24

Two shared settings make small multiples comparable. The same vmin/vmax across panels means a colour has the same meaning everywhere — without it each panel scales to its own range and the grid becomes eight unrelated maps. And the same axis limits, taken from the context layer, mean each panel covers the same ground, so a sparse category reads as sparse rather than as zoomed in.

A layered map is a composition

A map with several layers is a composition, and the things that go wrong are composition problems rather than GIS problems.

The CRS requirement comes from matplotlib knowing nothing about geography. gdf.plot(ax=ax) extracts coordinate arrays and draws them; it does not check whether two layers' numbers refer to the same space. Two layers in different CRS are simply two sets of numbers, and matplotlib will faithfully draw one at (430000, 430000) and the other at (-1.5, 53.8), as step 1's printout showed. The autoscaling then produces the characteristic symptom: a map that is 99% empty with two tiny clusters. Nothing raises, because from matplotlib's point of view nothing is wrong.

Occlusion is the second structural issue, and it is a consequence of the default fill. A filled polygon is opaque, so it hides everything beneath it. On a single-layer map that is fine. On a layered map it means the drawing order silently decides what is visible, and since the default order is call order, the map's content depends on the sequence of statements. zorder converts that from an accident into a decision, and facecolor="none" or alpha removes the conflict entirely for layers that are context rather than content.

The design principle underneath is worth stating: visual weight should follow importance, not size. A ward boundary is large and usually the least important thing on the map; a bus stop is a few pixels and often the point. Left to defaults, the large thing dominates. Giving context layers no fill, thin pale strokes and a low zorder, and content layers bright colours, white edges and a high zorder, inverts that in the right direction.

The extent problem is the same problem in a different dimension. Matplotlib's autoscale fits everything drawn, so the layer with the widest bounds sets the scale — usually a context layer, sometimes an error. Setting limits explicitly from the layer that defines the map's subject, after all drawing, gives a reproducible frame. Where extra features would also distort a legend or a classification, clipping the data rather than the view is the stronger fix, since it changes what is in the map rather than only what is visible.

Finally, the legend. GeoPandas' built-in legend describes a single column of a single layer, because that is all it can know about. A multi-layer legend is a manual object, and manual objects drift: someone changes a colour in the plot call and not in the handle. Building both from one dict, as Example 1 does, is a small amount of structure that removes an entire class of quietly wrong output — the same argument as for drawing the classification onto a choropleth.

ax=ax, axis limits, alpha and markersize

• ax=ax is required on every call after the first, or each layer opens its own figure.
• Set xlim/ylim after all plotting, since each plot() autoscales.
• alpha applies to fill and edge together. Plot the layer twice — fill then .boundary — for a solid outline over a translucent fill.
• .boundary gives a line GeoSeries, which takes color and linewidth rather than facecolor and edgecolor.
• markersize in GeoPandas is area in points², so doubling it does not double the visual diameter.
• Mixed geometry types in one layer get one style. Split by geom_type when they need different treatment.
• gpd.clip cuts geometries; a predicate filter keeps whole features. Different edges, different feature counts.
• Leave gaps in zorder so a new layer can be inserted without renumbering.
• plt.close(fig) in a loop, or a batch of maps exhausts memory — see how to batch-generate map images.
• facecolor="none" is a string, not None. None means "use the default", which is a solid fill.

Internal links

• How to plot maps in Python with GeoPandas and Matplotlib — the single-layer basics
• How to add a basemap to a GeoPandas map with contextily — the bottom layer of most maps
• How to add legends, labels and a scale bar to a GeoPandas map — finishing the composition
• Choosing a map projection for display — picking the shared CRS
• Choropleth classification explained — styling within a thematic layer
• How to clip spatial data in Python with GeoPandas — trimming layers to the map extent
• How to batch-generate a map image for every region — the same composition, many times
• Why is my GeoPandas plot blank or empty? — when a layer does not appear at all

FAQ

Why do my layers appear on separate figures?

Every plot() call creates a new figure unless you pass ax=. Create the axes once with plt.subplots() and pass them to every layer.

One layer is drawn miles away from the others.

It is in a different CRS. Print crs and total_bounds for every layer, then reproject them all to one CRS before plotting.

How do I stop a polygon layer hiding everything?

Give it facecolor="none" if it is context, or an alpha around 0.5 if it is thematic. Solid fills are for one layer per map.

How does drawing order work?

Later calls draw over earlier ones. Pass zorder on each layer to make the order explicit and independent of call sequence — higher numbers on top.

How do I build a legend for several layers?

By hand, with Patch for areas and Line2D for lines and points, then ax.legend(handles=[…]). Build the handles from the same styling dict you plotted with so they cannot drift apart.

Why does my map zoom out to include an empty area?

Matplotlib autoscales to fit everything drawn. Set xlim and ylim explicitly from the layer that defines your subject, after all layers are plotted.

Should I clip layers to the map extent?

Clip when the extra features would also affect a legend, a classification or a count. If it is only about what is visible, setting the axis limits is cheaper and keeps the geometries intact.