⚙️ How-To Maps & visualisation Geometry & topology Making Maps · step 2 of 13 Updated

How to Plot Multiple Layers on One Map in GeoPandas

Problem statement

Three layers, three calls, three separate figures:

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 if one layer is in a different CRS, it is drawn 5,000 km away and matplotlib helpfully zooms out to fit both, leaving two 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()
Stack showing drawing order from background polygons at the bottom to labels at the top.
Large and contextual at the bottom, small and important at the top.
Layer type zorder Style
basemap or background fill 1 pale, low contrast
area polygons 2 fill with alpha, or facecolor="none"
lines 3 thin, mid-grey
points 4 bright, white edge
labels and annotations 5 always last

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

Step-by-step solution

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 = {"boundary": boundary, "roads": roads, "sites": sites}
for name, layer in layers.items():
    print(f"{name:<10} {str(layer.crs):<14} {[round(v) for v in layer.total_bounds]}")
boundary   EPSG:27700     [351204, 381009, 407881, 445902]
roads      EPSG:27700     [351980, 382104, 407002, 445118]
sites      EPSG:4326      [-2, 53, -1, 53]

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

Grid comparing solid fill, alpha fill, outline-only and hatched fill for overlapping polygons.
Three of these let a layer beneath show through. The default does not.

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

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="Monitoring sites"),
    Patch(facecolor="none", edgecolor="#ef4444", hatch="///", label="Exclusion zone"),
]
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.

Code examples

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("wards.gpkg"),    {"facecolor": "#f8fafc"}),
    ("Green space", gpd.read_file("parks.gpkg"),   {"facecolor": "#dcfce7",
                                                    "edgecolor": "#86efac"}),
    ("Roads",      gpd.read_file("roads.gpkg"),    {"linewidth": 0.5}),
    ("Sites",      gpd.read_file("sites.gpkg"),    {"markersize": 50}),
], crs=27700, extent_from="Ward", title="Monitoring network")

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)

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("wards.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["highway"] == 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["highway"].isin(ROAD_STYLE)]
if len(unstyled):
    print(f"  ! {len(unstyled):,} features with an unhandled highway type: "
          f"{sorted(unstyled['highway'].unique())[:6]}")
ax.set_axis_off()
  ✓ motorway         204 features
  ✓ primary        1,882 features
  ✓ secondary      4,110 features
  ✓ residential   28,904 features
  ! 6,204 features with an unhandled highway type: ['footway', 'service', 'track']

The unstyled report is the part worth copying. Styling by dictionary lookup silently drops anything not in the dictionary, and 6,204 missing features 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

fig = small_multiples(gpd.read_file("incidents.gpkg"), "category",
                      context=gpd.read_file("boundary.gpkg"),
                      color="#ef4444", markersize=6, alpha=0.6)

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.

Explanation

Panels contrasting default styling, where the largest layer dominates, with styling where the most important layer does.
Defaults give visual weight to whichever layer is biggest. Usually that is the least interesting one.

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 (350000, 400000) and the other at (-2.5, 53.5). 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 monitoring site 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.

Edge cases or notes

  • 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.

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.