Chapter 15: Visualizing Spatial Data¶

Part II — Geographic Data Science

Learning Objectives¶

By the end of this chapter you will be able to:

• Create static maps with matplotlib and contextily
• Build interactive maps with folium and leafmap
• Apply cartographic best practices

In [28]:

# Standard imports — add chapter-specific imports below
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd

# Standard imports — add chapter-specific imports below import numpy as np import pandas as pd import matplotlib.pyplot as plt import geopandas as gpd

Every map in this book so far has been a quick .plot() — a sanity check, not something meant to be read by anyone else. This chapter treats the map as the deliverable: styled well enough to publish, interactive when that adds real value, and honest about what it's showing. Maranhão's municipalities, still the running dataset since Chapter 6, are the canvas throughout.

In [29]:

url = (
    "https://geoftp.ibge.gov.br/organizacao_do_territorio/"
    "malhas_territoriais/malhas_municipais/municipio_2022/"
    "UFs/MA/MA_Municipios_2022.zip"
)
maranhao = gpd.read_file(url)

url = ( "https://geoftp.ibge.gov.br/organizacao_do_territorio/" "malhas_territoriais/malhas_municipais/municipio_2022/" "UFs/MA/MA_Municipios_2022.zip" ) maranhao = gpd.read_file(url)

Static Maps with matplotlib¶

A bare gdf.plot() is a start, not a finished map. Four small additions turn it into something publishable: a meaningful column to color by, a legend, a title, and removing the latitude/longitude tick marks that rarely add anything to a thematic map:

In [30]:

fig, ax = plt.subplots(figsize=(8, 8))
maranhao.plot(
    column="AREA_KM2",
    cmap="YlOrRd",
    edgecolor="black",
    linewidth=0.3,
    legend=True,
    legend_kwds={"label": "Area (km²)", "shrink": 0.6},
    ax=ax,
)
ax.set_title("Maranhão municipalities by area")
ax.set_axis_off()
plt.show()

fig, ax = plt.subplots(figsize=(8, 8)) maranhao.plot( column="AREA_KM2", cmap="YlOrRd", edgecolor="black", linewidth=0.3, legend=True, legend_kwds={"label": "Area (km²)", "shrink": 0.6}, ax=ax, ) ax.set_title("Maranhão municipalities by area") ax.set_axis_off() plt.show()

ax.set_axis_off() is a small habit worth keeping for every thematic map from here on: the raw lat/lon tick values almost never help a reader, and removing them focuses attention on the shapes and colors that do.

Adding Basemaps with contextily¶

A map of municipality boundaries floating on a blank white background tells a reader nothing about where they are relative to roads, cities, or the coastline. contextily solves this by fetching map tiles from an online provider and drawing them underneath your own layer — but it has one hard requirement: tile providers speak Web Mercator (EPSG:3857), so anything plotted on top needs the same reprojection Chapter 6 and Chapter 7 already made a habit of checking first:

In [31]:

import contextily as cx

maranhao_web = maranhao.to_crs(epsg=3857)

fig, ax = plt.subplots(figsize=(8, 8))
maranhao_web.plot(ax=ax, facecolor="none", edgecolor="steelblue", linewidth=0.6)
cx.add_basemap(ax, source=cx.providers.Esri.WorldGrayCanvas)
ax.set_axis_off()
ax.set_title("Maranhão municipalities over a basemap")
plt.show()

import contextily as cx maranhao_web = maranhao.to_crs(epsg=3857) fig, ax = plt.subplots(figsize=(8, 8)) maranhao_web.plot(ax=ax, facecolor="none", edgecolor="steelblue", linewidth=0.6) cx.add_basemap(ax, source=cx.providers.Esri.WorldGrayCanvas) ax.set_axis_off() ax.set_title("Maranhão municipalities over a basemap") plt.show()

cx.providers lists dozens of tile sources beyond Esri.WorldGrayCanvas — satellite imagery, terrain shading, OpenStreetMap's own detailed style, and CARTO's Positron/DarkMatter among them, each with its own visual weight and its own access requirements. WorldGrayCanvas above was chosen deliberately: it's a light, low-contrast basemap that stays out of the way of whatever thematic layer sits on top of it, and Esri's free tier serves it without an API key or any special request headers. The two more obvious choices both have a catch worth knowing about, not worth getting stuck on: CARTO's Positron/DarkMatter now require a free API key (CARTO began enforcing this in late 2026 — sign up at carto.com/basemaps/apikey and pass it via cx.providers.CartoDB.Positron(api_key=...)), and OpenStreetMap.Mapnik — the "official" OSM style — now enforces its own tile usage policy strictly enough that a plain, unconfigured request from contextily can come back blocked, since OSM's volunteer-run servers require a proper Referer header identifying the requesting application. Neither is a dead end, just extra setup a first example shouldn't need.

Interactive Maps with folium¶

A static map is the right tool for a printed page or a fixed figure. When the reader needs to pan, zoom, or hover a specific municipality to read its name, an interactive map earns its extra weight. folium wraps the JavaScript mapping library Leaflet, and reading a GeoDataFrame into it is close to direct:

In [32]:

import folium

center = maranhao.geometry.centroid.union_all().centroid
m = folium.Map(location=[center.y, center.x], zoom_start=6, tiles="OpenStreetMap")

folium.GeoJson(
    maranhao,
    tooltip=folium.GeoJsonTooltip(fields=["NM_MUN", "AREA_KM2"]),
    style_function=lambda feature: {"fillColor": "steelblue", "fillOpacity": 0.2, "color": "black", "weight": 0.5},
).add_to(m)

m

import folium center = maranhao.geometry.centroid.union_all().centroid m = folium.Map(location=[center.y, center.x], zoom_start=6, tiles="OpenStreetMap") folium.GeoJson( maranhao, tooltip=folium.GeoJsonTooltip(fields=["NM_MUN", "AREA_KM2"]), style_function=lambda feature: {"fillColor": "steelblue", "fillOpacity": 0.2, "color": "black", "weight": 0.5}, ).add_to(m) m

/tmp/ipykernel_268850/1449308392.py:3: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  center = maranhao.geometry.centroid.union_all().centroid

Out[32]:

Make this Notebook Trusted to load map: File -> Trust Notebook

folium.GeoJsonTooltip is doing the interactive work here — hovering any municipality now shows its name and area, information a static map would need a hundred tiny labels to convey, most of them illegibly overlapping. m on its own, as the last line of a Jupyter cell, is enough to render the map inline; saving it for use outside a notebook is a one-line m.save("maranhao.html"). If the map doesn't render and Jupyter shows "Make this Notebook Trusted to load map" instead, that's a Jupyter security feature, not a bug: folium embeds raw HTML/JavaScript in the cell output, which Jupyter only renders automatically for notebooks it already trusts — one you wrote yourself in that session, not one opened from a downloaded file. File → Trust Notebook (or just re-running the cell, in most Jupyter versions) fixes it for the rest of the notebook.

Choropleth Maps¶

Every map so far has colored municipalities by a number already sitting in maranhao itself. A choropleth more commonly joins in an outside attribute first — annual population estimates published separately by IBGE are a natural fit here, merged on the same municipality code Chapter 7's Geoprocessing Workflows section already established a pattern for:

population = pd.read_csv(
    "estimativa_dou_2024.csv", sep=";", thousands=",",
    dtype={"COD. UF": str, "COD. MUNIC": str},
)
population["cod_mun"] = population["COD. UF"] + population["COD. MUNIC"]

maranhao_pop = pd.merge(maranhao, population, left_on="CD_MUN", right_on="cod_mun", how="left")

The interesting choice in a choropleth isn't the color palette — it's how the numeric range gets split into bins, which is exactly what mapclassify handles. Four common schemes answer that question differently:

• EqualInterval — splits the full value range into k equal-width bins. Simple, but a handful of outliers can compress every other municipality into one or two bins.
• Quantiles — puts an equal count of municipalities in each bin, regardless of value spread. Always uses every color, but bin boundaries can land at oddly specific numbers.
• FisherJenks (natural breaks) — finds the bin boundaries that best separate genuinely distinct clusters in the data itself, minimizing within-bin variance. Usually the most visually defensible default.

geopandas.plot() accepts a scheme argument directly, so the classification and the plotting happen in one call — this example uses maranhao's own AREA_KM2, already on hand without needing the population merge above:

In [33]:

fig, ax = plt.subplots(figsize=(8, 8))
maranhao.plot(
    column="AREA_KM2",
    scheme="FisherJenks",
    k=6,
    cmap="YlGnBu",
    legend=True,
    legend_kwds={"fmt": "{:.0f}", "title": "Area (km²)"},
    edgecolor="black",
    linewidth=0.2,
    ax=ax,
)
ax.set_axis_off()
ax.set_title("Maranhão municipalities — natural breaks classification")
plt.show()

fig, ax = plt.subplots(figsize=(8, 8)) maranhao.plot( column="AREA_KM2", scheme="FisherJenks", k=6, cmap="YlGnBu", legend=True, legend_kwds={"fmt": "{:.0f}", "title": "Area (km²)"}, edgecolor="black", linewidth=0.2, ax=ax, ) ax.set_axis_off() ax.set_title("Maranhão municipalities — natural breaks classification") plt.show()

Swap scheme="FisherJenks" for "EqualInterval" or "Quantiles" and rerun — the same six colors will cover noticeably different sets of municipalities each time, which is the whole point of comparing schemes before publishing a choropleth rather than accepting whichever one happens to be the default.

Visualizing Simulation Output¶

Every technique above visualizes a static dataset — one value per municipality, fixed in time. A simulation produces something different: the same grid's state changing step by step, which a single map can't show at all. Two approaches cover most needs.

A small multiple — several snapshots side by side — works well for a handful of key time steps:

In [34]:

# Illustrative: a toy state array evolving over a few steps, not a real simulation
np.random.seed(0)
n_steps = 4
grid_side = 15
states = [np.random.rand(grid_side, grid_side)]
for _ in range(n_steps - 1):
    states.append(np.clip(states[-1] + np.random.normal(0, 0.1, states[-1].shape), 0, 1))

fig, axes = plt.subplots(1, n_steps, figsize=(14, 4))
for t, (ax, state) in enumerate(zip(axes, states)):
    ax.imshow(state, cmap="viridis", vmin=0, vmax=1)
    ax.set_title(f"t = {t}")
    ax.axis("off")
plt.show()

# Illustrative: a toy state array evolving over a few steps, not a real simulation np.random.seed(0) n_steps = 4 grid_side = 15 states = [np.random.rand(grid_side, grid_side)] for _ in range(n_steps - 1): states.append(np.clip(states[-1] + np.random.normal(0, 0.1, states[-1].shape), 0, 1)) fig, axes = plt.subplots(1, n_steps, figsize=(14, 4)) for t, (ax, state) in enumerate(zip(axes, states)): ax.imshow(state, cmap="viridis", vmin=0, vmax=1) ax.set_title(f"t = {t}") ax.axis("off") plt.show()

For a full run rather than a handful of snapshots, matplotlib.animation.FuncAnimation turns the same sequence into a genuine animation — one frame per time step, played back at a chosen speed — which is precisely what DisSModel's own Map and Chart visualization classes automate for you during env.run(): attach one to a model, and every step() call updates the figure automatically, without hand-writing a FuncAnimation callback yourself. Chapter 21 shows that mechanism directly, once a real model exists to visualize.

Did you know?

All of this chapter's maps have used at most a few hundred geometries — fine for matplotlib and folium alike. A dataset with millions of points (every building centroid in a country, say) renders both libraries unusably slow. datashader takes a different approach entirely: instead of drawing one shape per point, it rasterizes the whole collection directly into a density grid, the same way a histogram bins numbers instead of drawing one bar per data point — genuinely fast even at tens of millions of points, at the cost of losing the ability to click an individual one.

Exercises¶

1. A different color scheme. Redo the Choropleth Maps example with scheme="Quantiles" and scheme="EqualInterval", keeping k=6 both times. Pick the São Luís municipality specifically — does it land in the same color bin under all three schemes, or does its classification change?
2. A different basemap. Rerun Adding Basemaps with contextily with cx.providers.Esri.WorldImagery (satellite) instead of WorldGrayCanvas. Does the extra visual detail help or compete with the municipality boundaries drawn on top?
3. Tooltip, extended. Add a third field to the folium.GeoJsonTooltip in Interactive Maps with folium — any column already present in maranhao. Confirm it appears on hover.
4. From snapshot to animation. Using the toy states list from Visualizing Simulation Output, look up matplotlib.animation.FuncAnimation and sketch (in comments, no need to fully run it) how you'd turn those four static frames into a playable animation instead.

In [ ]:

# Your code here

# Your code here

Summary¶

Key concepts introduced¶

• Publishable static maps with matplotlib: a meaningful column, a legend, and ax.set_axis_off()
• contextily basemaps, and the hard requirement that both layers share Web Mercator (EPSG:3857)
• Interactive maps with folium, including GeoJsonTooltip for hover-based attribute lookup
• Choropleth classification with mapclassify — EqualInterval, Quantiles, and FisherJenks (natural breaks) each tell a visually different story from the same numbers
• Visualizing simulation output over time: small multiples for a few key steps, animation for a full run, and DisSModel's Map/Chart classes as the automated version of the same idea
• datashader as the answer once a dataset's point count outgrows what matplotlib or folium can render interactively

This closes Part II's tour of the geographic data science toolkit. Chapter 16 picks up the one gap deliberately left open since Chapter 6: what happens when a vector layer and a raster layer need to work together in the same analysis.

Further Reading¶

• contextily documentation, Adding a background map to plots: https://contextily.readthedocs.io/en/latest/intro_guide.html
• folium documentation: https://python-visualization.github.io/folium/latest/
• mapclassify documentation, Choosing a classification scheme: https://pysal.org/mapclassify/
• datashader documentation, Getting Started: https://datashader.org/getting_started/index.html