How to Add a Basemap to a GeoPandas Map with contextily

Problem statement

Your polygons float on white. Nobody can tell whether they are in Manchester or Marseille:

import matplotlib
matplotlib.use("Agg")

gdf.plot()

The obvious fix is one line, and it usually fails on the first attempt:

import warnings

import contextily as cx
from pyproj import Transformer

ax = gdf.plot()
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    cx.add_basemap(ax)
x0, x1, y0, y1 = ax.images[0].get_extent()
lon, lat = Transformer.from_crs(3857, 4326, always_xy=True).transform((x0 + x1) / 2, (y0 + y1) / 2)
print(f"the data: Leeds, in EPSG:27700; the tiles fetched: centred on {lat:.2f}°N, {lon:.2f}°E")

the data: Leeds, in EPSG:27700; the tiles fetched: centred on 3.95°N, 3.87°E

Either nothing appears, or the tiles land somewhere else entirely, or you get:

HTTPError: Tile URL resulted in a 404 error

or the map draws a grey blur because the zoom level was inferred from a bounding box in degrees when the tiles expect metres.

That is what happened above. contextily read the ward coordinates — metres east and north of a point off the Isles of Scilly — as Web Mercator metres from the equator and the prime meridian, and fetched tiles for a patch of the Gulf of Guinea. Almost all of it comes from one requirement: tile basemaps are served in Web Mercator (EPSG:3857), and your data probably is not.

Quick answer

import geopandas as gpd
import contextily as cx
import matplotlib.pyplot as plt

gdf = gpd.read_file("wards.gpkg").to_crs(3857)      # ← the whole fix

fig, ax = plt.subplots(figsize=(10, 10))
gdf.plot(ax=ax, alpha=0.6, edgecolor="white", linewidth=0.5)
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron)
ax.set_axis_off()
fig.savefig("quick.png", dpi=80)

If reprojecting the data is not acceptable, reproject the tiles instead:

ax = gdf.plot()                                    # data stays in EPSG:27700
cx.add_basemap(ax, crs=gdf.crs, source=cx.providers.CartoDB.Positron)

crs= warps the fetched tiles into your CRS. It is slower and slightly blurry, and it is the right answer when the map's coordinates matter.

EPSG:27700
[-200422. 7113349. -143641. 7159912.]

Web Mercator or crs=, data first, provider, zoom, cache

1. Reproject to Web Mercator — or tell contextily not to

Every XYZ tile service uses EPSG:3857. add_basemap assumes the axes are already in it unless told otherwise:

gdf = gpd.read_file("wards.gpkg")
print(gdf.crs)
gdf3857 = gdf.to_crs(3857)
print(gdf3857.total_bounds.round(0))

EPSG:27700
[-200422. 7113349. -143641. 7159912.]

Web Mercator coordinates are metres from the intersection of the equator and the prime meridian, running to about ±20,037,508. Numbers in the millions with that magnitude are the sign you are in 3857.

The alternative — crs=gdf.crs — makes contextily fetch tiles and reproject them:

ax = gdf.plot(alpha=0.6)
cx.add_basemap(ax, crs=gdf.crs, source=cx.providers.CartoDB.Positron)

Reproject the data unless the axis coordinates themselves matter — for a map with a scale bar in metres, or one whose coordinates a reader will use, keep your CRS and warp the tiles.

2. Plot the data first, then the basemap

add_basemap reads the axes' current extent to decide which tiles to fetch. On an empty axis there is no extent:

# ❌ nothing to work from
fig, ax = plt.subplots()
cx.add_basemap(ax)                         # fetches tiles for the whole world, or fails

# ✅ the data sets the extent
fig, ax = plt.subplots()
gdf.plot(ax=ax)
cx.add_basemap(ax)

The basemap is drawn on top of what is already there, so make the data layers semi-transparent, or push the basemap behind:

gdf.plot(ax=ax, alpha=0.6, zorder=2)
cx.add_basemap(ax, zorder=1)

zorder is the explicit form and worth using once a map has more than two layers — see how to plot multiple layers on one map.

3. Choose a provider that suits the map

import contextily as cx

for name in ["CartoDB.Positron", "CartoDB.DarkMatter", "OpenStreetMap.Mapnik",
             "Esri.WorldImagery", "Esri.WorldTopoMap", "CartoDB.VoyagerNoLabels"]:
    provider = cx.providers.query_name(name)
    key = " · needs an API key" if "apikey" in provider else ""
    print(f"{name:<24} {provider.get('attribution', '')[:40]}{key}")

CartoDB.Positron         (C) OpenStreetMap contributors (C) CARTO · needs an API key
CartoDB.DarkMatter       (C) OpenStreetMap contributors (C) CARTO · needs an API key
OpenStreetMap.Mapnik     (C) OpenStreetMap contributors
Esri.WorldImagery        Tiles (C) Esri -- Source: Esri, i-cubed,
Esri.WorldTopoMap        Tiles (C) Esri -- Esri, DeLorme, NAVTEQ,
CartoDB.VoyagerNoLabels  (C) OpenStreetMap contributors (C) CARTO · needs an API key

The CARTO entries are the ones the xyzservices catalogue now marks as needing a key — its 2026 releases warn "CARTO tiles now require an API key" every time one is used. When this page was run the tiles still downloaded without one, but a pipeline that depends on that is one change away from failing. Register for a key, or pick a provider that does not need one, before anything scheduled relies on Positron.

For a choropleth, a busy basemap competes with the data. Positron exists precisely for this, and PositronNoLabels when you add your own labels.

Attribution is a licence condition, not a courtesy. OpenStreetMap tiles require it:

provider = cx.providers.CartoDB.Positron
cx.add_basemap(ax, source=provider)
ax.annotate(provider["attribution"], xy=(0.005, 0.005), xycoords="axes fraction",
            fontsize=6, color="#334155", backgroundcolor="#ffffffaa")

add_basemap adds an attribution line by default; attribution=False turns it off, and attribution="…" overrides the text. Do not turn it off without a reason.

4. Control the zoom level

cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, zoom=12)

Zoom runs 0 (the whole world in one tile) to about 19 (building level). contextily infers it from the extent, and the inference is usually right — but worth overriding when:

• Text is too small or too large. Tiles are raster images with baked-in labels, so their text does not scale with your figure. A higher zoom gives smaller, denser labels.
• The output is for print. A screen-resolution basemap looks soft at 300 dpi. Raise the zoom by one or two and increase the figure size to match.
• You are fetching too many tiles. Each level quadruples the count. Zoom 16 over a city can be thousands of requests.

import math

def tile_count(bounds_3857, zoom):
    """How many 256-px tiles this extent needs at this zoom."""
    world = 20037508.342789244 * 2
    span = world / (2 ** zoom)
    minx, miny, maxx, maxy = bounds_3857
    return (math.ceil((maxx - minx) / span) + 1) * (math.ceil((maxy - miny) / span) + 1)

for z in range(8, 17, 2):
    print(f"zoom {z:>2}  {tile_count(gdf3857.total_bounds, z):>7,} tiles")

zoom  8        4 tiles
zoom 10        9 tiles
zoom 12       42 tiles
zoom 14      525 tiles
zoom 16    7,332 tiles

Thousands of tile requests for one map of one city is abusive to a free service and slow for you. contextily caches to disk, so the second run is cheap — but the first is not, and public tile servers rate-limit.

5. Cache, and be a good citizen

cx.set_cache_dir("./.tilecache")           # persists between runs

Without a cache directory, every run re-fetches. With one, a repeated map is instant and no requests leave your machine.

The OpenStreetMap Foundation's tile usage policy prohibits bulk downloading and heavy automated use. For a batch job producing hundreds of maps, either use a commercial provider with an API key, or fetch the extent once and reuse the image:

# fetch once, reuse for every map in the batch
import os


img, ext = cx.bounds2img(*gdf3857.total_bounds, zoom=12,
                         source=cx.providers.CartoDB.Positron)
os.makedirs("out", exist_ok=True)
for name, part in gdf3857.groupby("ward_name"):
    fig, ax = plt.subplots(figsize=(8, 8))
    ax.imshow(img, extent=ext)             # the cached array, not a new request
    part.plot(ax=ax, facecolor="none", edgecolor="red", linewidth=2)
    ax.set_xlim(ext[0], ext[1]); ax.set_ylim(ext[2], ext[3])
    ax.set_axis_off()
    fig.savefig(f"out/{name}.png", dpi=80, bbox_inches="tight")
    plt.close(fig)
print(f"{len(os.listdir('out'))} ward maps from one {img.shape[1]}×{img.shape[0]} px basemap fetch")

33 ward maps from one 1792×1280 px basemap fetch

Leeds wards and bus stops on tiles

The data used on this page

Leeds's 33 electoral wards (ONS, December 2024) in British National Grid, with 2021 Census population density from Nomis; and, for example 2, the active NaPTAN bus stops in the city-centre ward of Little London & Woodhouse. Both are written in EPSG:27700, the CRS they would arrive in from a British source, which is what makes the basemap go wrong.

import io

import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import requests

ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
gdf = gpd.read_file(f"{ONS}/Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
                    "?where=LAD24NM%3D%27Leeds%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson")
gdf = gdf.set_crs(27700, allow_override=True).rename(columns={"WD24CD": "ward_code", "WD24NM": "ward_name"})
pop = pd.read_csv(io.StringIO(requests.get("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv", timeout=120,
                  params={"date": "latest", "geography": ",".join(gdf["ward_code"]), "c2021_restype_3": 0,
                          "measures": 20100, "select": "geography_code,obs_value"}).text))
gdf["density"] = gdf["ward_code"].map(pop.set_index("GEOGRAPHY_CODE")["OBS_VALUE"]) / (gdf.area / 1e6)
gdf.to_file("wards.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=27700,
                         geometry=gpd.points_from_xy(naptan["Easting"], naptan["Northing"]))
centre = gdf[gdf["ward_name"] == "Little London & Woodhouse"].union_all()
stops[stops.within(centre)].to_file("sites.gpkg")
print(f"{len(gdf)} wards, {len(gpd.read_file('sites.gpkg'))} bus stops in the city-centre ward")

33 wards, 178 bus stops in the city-centre ward

Example 1: a reusable basemap function

import geopandas as gpd
import contextily as cx
import matplotlib.pyplot as plt

cx.set_cache_dir("./.tilecache")

PROVIDERS = {
    "light": cx.providers.CartoDB.Positron,
    "light_plain": cx.providers.CartoDB.PositronNoLabels,
    "dark": cx.providers.CartoDB.DarkMatter,
    "osm": cx.providers.OpenStreetMap.Mapnik,
    "satellite": cx.providers.Esri.WorldImagery,
    "terrain": cx.providers.Esri.WorldTopoMap,
}

def map_with_basemap(gdf, *, style="light", zoom=None, pad=0.08, figsize=(10, 10),
                     reproject=True, **plot_kwargs):
    """Plot a layer over a tile basemap, handling the CRS and the extent."""
    source = PROVIDERS.get(style, style)

    if reproject:
        layer, tile_crs = gdf.to_crs(3857), None
    else:
        layer, tile_crs = gdf, gdf.crs

    fig, ax = plt.subplots(figsize=figsize)
    plot_kwargs.setdefault("alpha", 0.65)
    plot_kwargs.setdefault("edgecolor", "white")
    plot_kwargs.setdefault("linewidth", 0.5)
    layer.plot(ax=ax, zorder=2, **plot_kwargs)

    # pad the extent so features are not flush against the frame
    minx, miny, maxx, maxy = layer.total_bounds
    dx, dy = (maxx - minx) * pad, (maxy - miny) * pad
    ax.set_xlim(minx - dx, maxx + dx)
    ax.set_ylim(miny - dy, maxy + dy)

    cx.add_basemap(ax, source=source, zoom=zoom or "auto", crs=tile_crs, zorder=1)
    ax.set_axis_off()
    return fig, ax

fig, ax = map_with_basemap(
    gpd.read_file("wards.gpkg"),
    style="light", column="density", scheme="quantiles", k=5,
    cmap="YlOrRd", legend=True,
)
fig.savefig("density.png", dpi=100, bbox_inches="tight")

Two details repay themselves. Setting the limits after plotting and before add_basemap means the basemap covers the padded extent rather than the tight data bounds, so there is no white border. And zorder is explicit on both layers, so adding a third later cannot silently reorder them.

Example 2: keeping your own CRS

When the axes must be in a projected national grid — for a scale bar, or coordinates a reader will use:

import geopandas as gpd, contextily as cx, matplotlib.pyplot as plt

gdf = gpd.read_file("sites.gpkg").to_crs(27700)      # British National Grid

fig, ax = plt.subplots(figsize=(10, 10))
gdf.plot(ax=ax, color="#ef4444", markersize=40, zorder=3)

cx.add_basemap(ax, crs=gdf.crs,                       # ← warp the tiles, not the data
               source=cx.providers.CartoDB.Positron, zoom=13, zorder=1)

# now a scale bar is meaningful, because the axes really are in metres
from matplotlib.patches import Rectangle
x0, x1 = ax.get_xlim(); y0, y1 = ax.get_ylim()
bar_m = 1000
bx, by = x0 + (x1 - x0) * 0.05, y0 + (y1 - y0) * 0.05
ax.add_patch(Rectangle((bx, by), bar_m, (y1 - y0) * 0.006,
                       facecolor="black", zorder=4))
ax.text(bx + bar_m / 2, by + (y1 - y0) * 0.012, f"{bar_m//1000} km",
        ha="center", fontsize=9, zorder=4)
ax.set_axis_off()

This is the case that justifies the slower path. A scale bar drawn on Web Mercator axes is wrong by a factor of 1/cos(latitude) — 1.69 at Leeds, so a "1 km" bar is really about 590 m on the ground — because Mercator's scale varies with latitude. Keeping the axes in EPSG:27700 makes one metre on the axis one metre on the ground, so the bar is honest. The distortion argument is in choosing a map projection for display.

Example 3: diagnosing a basemap that will not appear

import contextily as cx
import geopandas as gpd
import matplotlib.pyplot as plt

def diagnose_basemap(gdf, source=cx.providers.CartoDB.Positron, zoom=None):
    print(f"data CRS      {gdf.crs}")
    print(f"data bounds   {gdf.total_bounds.round(1).tolist()}")

    if gdf.crs is None:
        return "no CRS — set one before adding a basemap"
    if gdf.crs.to_epsg() != 3857:
        print(f"  → not Web Mercator; reproject or pass crs={gdf.crs}")

    g = gdf.to_crs(3857)
    minx, miny, maxx, maxy = g.total_bounds
    LIMIT = 20037508.342789244
    if not (-LIMIT <= minx <= LIMIT and -LIMIT <= maxx <= LIMIT):
        return f"bounds outside the Web Mercator domain: {g.total_bounds}"
    if abs(gdf.to_crs(4326).total_bounds[1]) > 85 or abs(gdf.to_crs(4326).total_bounds[3]) > 85:
        return "data reaches beyond ±85° latitude — Web Mercator is undefined there"

    try:
        img, ext = cx.bounds2img(minx, miny, maxx, maxy,
                                 zoom=zoom or "auto", source=source)
        print(f"  ✓ fetched {img.shape[1]}×{img.shape[0]} px covering {[round(e) for e in ext]}")
    except Exception as exc:
        return f"tile fetch failed: {type(exc).__name__}: {exc}"

    fig, ax = plt.subplots(figsize=(8, 8))
    ax.imshow(img, extent=ext)
    g.boundary.plot(ax=ax, color="red", linewidth=1.5)
    ax.set_title("basemap and data, drawn from the same extent")
    return "ok"

print(diagnose_basemap(gpd.read_file("wards.gpkg")))

data CRS      EPSG:27700
data bounds   [413216.0, 422593.4, 446879.0, 450175.3]
  → not Web Mercator; reproject or pass crs=EPSG:27700
  ✓ fetched 1024×768 px covering [-215247, -136975, 7103140, 7161844]
ok

The checks are ordered by how badly each one breaks things: no CRS makes every later step meaningless, out-of-domain bounds make the tile request nonsensical, and only then is it worth actually fetching. Calling bounds2img directly separates "the tiles cannot be fetched" from "the tiles were fetched and drawn in the wrong place", which are different problems with different fixes.

A basemap is a pre-rendered raster pyramid

A tile basemap is a pre-rendered raster pyramid, and every constraint in this article follows from that.

The tiles were drawn in advance, in Web Mercator, at fixed zoom levels, in 256-pixel squares whose positions are defined by a scheme every provider shares: at zoom $z$ the world is a $2^z \times 2^z$ grid, and tile (x, y, z) covers a known rectangle of Web Mercator space. To place a tile on your axes, contextily needs those axes to be in the same coordinate space. If they are not — if they are in metres from a national grid origin, or in degrees — the tile's coordinates mean nothing there, and it lands wherever those numbers happen to fall.

This is why EPSG:3857 is not negotiable. It is not that contextily prefers it; it is that the tiles were rendered in it and cannot be re-rendered. The crs= argument does the only other thing possible: fetch the tiles, warp the resulting image into your CRS, and draw the warped version. That works, costs a resampling step, and softens the labels slightly — the same trade as any raster reprojection.

Zoom level is a resolution choice, not a scale choice. Because tiles are raster images with labels baked in, their text size is fixed in pixels. Rendering a zoom-10 tile across a large figure enlarges its labels along with everything else, which is why an under-zoomed basemap looks like a blurry poster. Going the other way, an over-zoomed basemap has tiny labels and costs 4× the requests per level. contextily's automatic choice targets a sensible pixel density for the axes; override it when the output medium differs — print at 300 dpi wants one or two levels more than a screen.

The extent has to exist before the tiles are chosen, which is why order matters. add_basemap reads ax.get_xlim() and ax.get_ylim() to decide which tiles to fetch, so calling it on an empty axis asks for the default extent — the whole world, or nothing. Setting the limits explicitly before adding the basemap, as in Example 1, also solves the related annoyance where the basemap stops exactly at the data's bounding box and leaves a white margin.

Finally, there is an ethical constraint that is easy to overlook because the code does not enforce it. Public tile services — OpenStreetMap most of all — are volunteer-funded, and their usage policies prohibit bulk downloading. At zoom 16 a single map of Leeds needs 7,332 tiles (step 4), and a batch job without a cache repeats that for every map it draws. Setting a cache directory, fetching one image and reusing it across a batch, or paying for a commercial provider are the three honest answers, and the first two are one line each.

Cache directories, attribution and zoom

• cx.set_cache_dir() should be one of the first lines in any script that draws maps. Without it, every run re-fetches.

• Plot data before add_basemap, or the axes have no extent to fetch tiles for.

• attribution=False removes a licence-required credit. Do not use it on published output without an alternative credit.

• zoom="auto" is the default; pass an integer for print, and check the tile count first.

• crs= warps the tiles, which is slower and slightly blurry but keeps your axes in your CRS.

• Beyond ±85° latitude Web Mercator is undefined. Clip before reprojecting.

• cx.providers.query_name("CartoDB.Positron") looks a provider up by string, useful for config-driven scripts.

• Some providers need an API key, supplied through the provider dict: cx.providers.Stadia.AlidadeSmooth(api_key="…"). The xyzservices catalogue now lists CARTO's basemaps among them (cx.providers.CartoDB.Positron(apikey="…")).

• A 404 usually means a dead provider or an unsupported zoom, not a bug in your code. Try another provider.

• cx.bounds2img returns an array and an extent, which is the right primitive for batch work and for caching by hand.

Internal links

• How to plot maps in Python with GeoPandas and Matplotlib — the plotting layer underneath
• Choosing a map projection for display — why tiles force Web Mercator
• How to plot multiple layers on one map in GeoPandas — z-order and layer discipline
• How to add legends, labels and a scale bar to a GeoPandas map — finishing the map
• How to batch-generate a map image for every region — reusing one basemap fetch
• How to reproject spatial data in Python (GeoPandas) — the to_crs call this depends on
• How to save a map as an image in Python with Matplotlib — export settings for print
• How to create interactive maps in Python with Folium — when a static basemap is not enough

FAQ

Why does my basemap appear in the wrong place?

The data is not in EPSG:3857. Either call gdf.to_crs(3857) before plotting, or pass crs=gdf.crs to add_basemap so the tiles are warped into your CRS instead.

Should I reproject the data or the tiles?

Reproject the data for display-only maps — the tiles stay sharp. Warp the tiles when the axes themselves must be in your CRS, for instance so a scale bar in metres is correct.

Why is the basemap blurry?

The zoom level is too low for the figure size, so tiles are being stretched. Pass a higher zoom=, and raise it one or two more for print output at 300 dpi.

Which provider should I use for a choropleth?

CartoDB.Positron, or PositronNoLabels if you add your own labels. A detailed basemap competes with the data for attention.

How do I avoid hammering the tile server?

Call cx.set_cache_dir(), keep the zoom no higher than needed, and for batch work fetch one image with cx.bounds2img and reuse it. Public tile services prohibit bulk downloading.

Do I have to show the attribution?

Yes for OpenStreetMap and most free providers — it is a licence condition. add_basemap adds it by default; leave it on unless you credit the source another way.

My basemap hides my data.

The basemap is drawn after and therefore on top. Pass zorder=1 to the basemap and a higher zorder to the data layers, and give the data some transparency.