Folium Choropleth from a GeoDataFrame

This guide builds a binned, value-shaded choropleth directly from a GeoDataFrame, the way an analyst actually has their data — not from the separate GeoJSON-plus-CSV inputs Folium's older API assumed. It is written for anyone who has a polygon layer with a numeric column and wants a clean interactive map. It sits under Interactive Maps with Folium in the Web Mapping & Interactive Visualization domain.

Why This Approach / What Goes Wrong

Folium ships two ways to make a choropleth. The legacy folium.Choropleth class wants a GeoJSON file and a separate GeoPandas DataFrame joined on a key_on string — a fragile setup where a single key mismatch produces a map of uniformly grey polygons with no error. The modern approach binds a single GeoDataFrame through folium.GeoJson with a classified style_function, so the geometry and the value travel together and there is no join to break.

Worth understanding why the legacy join fails silently, because the same shape of bug appears whenever a map is assembled from two inputs. key_on is not a column name — it is a JavaScript-style accessor path evaluated against each GeoJSON feature, typically "feature.properties.tract_id". Folium walks the features, resolves that path, and looks the resulting string up in the first column of the data frame. Anything that fails to match is not an error; it is simply absent from the colour dictionary, so the polygon is painted with nan_fill_color (black at 40 % opacity by default) and the map still renders. Leading zeros stripped by a CSV reader, an integer key on one side and a string on the other, a trailing space from a shapefile's fixed-width attribute — every one of those produces a plausible-looking map that is quietly wrong. Binding a single GeoDataFrame removes the class of bug entirely.

The other common failure is classification. A linear color ramp on skewed data (population, income, emissions) paints almost everything the same shade because a few outliers stretch the scale. Quantile or natural-breaks classification via mapclassify fixes this by binning on the data distribution instead of its raw range.

The third failure is conceptual and no library will catch it: a choropleth must show an intensive quantity. Shading polygons by a raw count — total population, number of permits, incident tallies — draws a map of polygon size, because a large tract contains more of everything. Convert to a rate or a density before classifying: incidents per thousand residents, permits per square kilometre, income as a per-household median. Where the denominator is area, compute that area in an equal-area projection rather than the projected CRS you happen to have, since a UTM zone's area distortion grows with distance from its central meridian and Web Mercator's is severe at high latitude.

Prerequisites

• geopandas>=0.14
• folium>=0.16
• mapclassify>=2.6 (required for binning schemes; Folium uses it silently)
• branca>=0.7 (color maps and the legend; installed with Folium)

conda install -c conda-forge "geopandas=0.14.*" "folium=0.16.*" "mapclassify=2.6.*"

Step-by-Step Implementation

1. Load the layer and reproject for the web. Compute any metric column in a projected CRS upstream, then reproject to WGS84 for display — Folium renders on a Leaflet slippy map and expects geographic EPSG:4326. Never derive areas, densities, or distances from EPSG:4326 itself; do that work in the projected CRS first, as covered in Coordinate Systems with PyProj.

import geopandas as gpd

# census_tracts carries a "median_income" column, analysed in EPSG:32610 upstream
census_tracts = gpd.read_file("census_tracts.gpkg")
census_tracts = census_tracts.to_crs(epsg=4326)  # Folium needs WGS84 (lon/lat)

2. Classify the values with mapclassify. Quantiles give visually balanced bins on skewed data.

import mapclassify

classifier = mapclassify.Quantiles(census_tracts["median_income"], k=5)
census_tracts["bin"] = classifier.yb  # integer bin index 0..4

3. Map bins to a color ramp. Use a sequential branca colormap stepped into the same number of classes.

import branca.colormap as cm

palette = ["#f0ebd8", "#9bb3c9", "#748cab", "#3e5c76", "#1d2d44"]
colormap = cm.StepColormap(
    palette,
    vmin=census_tracts["median_income"].min(),
    vmax=census_tracts["median_income"].max(),
    index=[classifier.bins[i] for i in range(-1, len(classifier.bins))][1:],
    caption="Median income (USD)",
)

4. Build the map with a classified style function. Folium takes coordinates latitude-first, so the map centre reads [centroid.y, centroid.x] — a reversal that silently drops the map in the ocean if you swap it.

import folium

center = census_tracts.geometry.union_all().centroid
income_map = folium.Map(location=[center.y, center.x], zoom_start=11, tiles="CartoDB positron")

folium.GeoJson(
    census_tracts,
    name="Median income",
    style_function=lambda feat: {
        "fillColor": palette[feat["properties"]["bin"]],
        "color": "#1d2d44",
        "weight": 0.6,
        "fillOpacity": 0.78,
    },
    tooltip=folium.GeoJsonTooltip(
        fields=["tract_id", "median_income"],
        aliases=["Tract", "Median income"],
    ),
).add_to(income_map)

colormap.add_to(income_map)        # legend
folium.LayerControl().add_to(income_map)
income_map.save("income_choropleth.html")

5. Give missing data its own colour. Tracts with a null value are not "low" — they are unknown, and painting them the bottom class of the ramp is a factual error the reader cannot detect. mapclassify will not bin nulls at all, so decide their fate explicitly: assign a sentinel bin, shade it a neutral grey that is outside the sequential ramp, and say so in the legend caption.

import numpy as np

NO_DATA = "#c3d0e4"          # deliberately outside the sequential palette

known = census_tracts["median_income"].notna()
census_tracts["bin"] = np.where(known, census_tracts["bin"], -1).astype(int)

def income_style(feature):
    bin_index = feature["properties"]["bin"]
    return {
        "fillColor": NO_DATA if bin_index < 0 else palette[bin_index],
        "color": "#1d2d44",
        "weight": 0.6,
        "fillOpacity": 0.45 if bin_index < 0 else 0.78,
    }

Reclassify on the non-null subset before this step, not after — running mapclassify.Quantiles over a column containing NaN either raises or silently shifts every break, depending on the version, and a shifted break set produces a map that looks fine and reports the wrong quintiles.

6. Make the tooltip carry the number, formatted. A choropleth answers "where", and the tooltip answers "how much". Raw JSON numbers render as 118400.0, which reads badly and buries the units. GeoJsonTooltip accepts aliases for the labels, localize=True to apply the browser's thousands separators, and sticky=True so the box follows the cursor rather than pinning to the polygon centroid.

tooltip = folium.GeoJsonTooltip(
    fields=["tract_id", "median_income", "households"],
    aliases=["Tract", "Median income (USD)", "Households"],
    localize=True,
    sticky=True,
    labels=True,
    style=(
        "background-color: #f7fafe; border: 1px solid #c3d0e4;"
        "border-radius: 4px; padding: 6px; font-size: 12px;"
    ),
)

Any field named here must still exist on the GeoDataFrame at render time. A common self-inflicted wound is trimming columns for payload size after building the tooltip: Folium raises no error for a missing field, it just renders an empty row in the box.

7. Revisit the scheme once you can see the map. Quantiles is a safe default, not a universal answer, and mapclassify makes the comparison cheap. Each classifier exposes .bins (the class upper bounds), .counts (members per class) and .adcm — the absolute deviation around class medians, a goodness-of-fit measure where lower means the classes fit the distribution more tightly. Fit two or three and look at the numbers before committing.

import mapclassify

values = census_tracts.loc[known, "median_income"]

for scheme in ("Quantiles", "NaturalBreaks", "EqualInterval", "HeadTailBreaks"):
    cls = mapclassify.classify(values, scheme, k=5) if scheme != "HeadTailBreaks" \
        else mapclassify.classify(values, scheme)
    print(f"{scheme:<16} k={len(cls.counts)}  adcm={cls.adcm:>12,.0f}  counts={list(cls.counts)}")

# Quantiles        k=5  adcm=  8,214,900  counts=[104, 103, 104, 103, 104]
# NaturalBreaks    k=5  adcm=  6,001,340  counts=[171, 186, 108, 42, 11]
# EqualInterval    k=5  adcm= 11,880,205  counts=[126, 262, 74, 44, 12]
# HeadTailBreaks   k=4  adcm=  7,442,118  counts=[398, 92, 22, 6]

Natural breaks usually wins on adcm because it is explicitly minimising within-class variance, but it buys that fit with wildly uneven class sizes — here the top class holds eleven tracts, so four fifths of the colour ramp describe under a tenth of the data. Quantiles gives every shade equal weight, which is what you want when the map is meant to be read comparatively (which tracts are in the top fifth?). Reach for NaturalBreaks when the distribution has genuine gaps you want the map to expose, UserDefined when the thresholds are policy rather than statistics (a poverty line, a regulatory limit), and HeadTailBreaks for heavy-tailed data where the interesting structure is in the tail. Note that FisherJenks is the optimal-partition classifier and is O(n²k): it is fine for a few thousand features and unusably slow beyond that, where FisherJenksSampled fits the breaks on a random sample instead.

If you want none of this control and just need a fast look, GeoPandas ships the whole pipeline as one call — census_tracts.explore(column="median_income", scheme="quantiles", k=5, legend=True) returns a Folium Map object you can keep adding layers to. It is excellent for exploration and deliberately opinionated about styling, which is exactly why production maps end up back on the explicit style_function written here.

Verification

Confirm the bins are populated and the output is reasonable before sharing it.

# Each quantile bin should hold a similar count of tracts
print(census_tracts["bin"].value_counts().sort_index())
# 0    104
# 1    103
# 2    104
# 3    103
# 4    104

# The saved file should be well under the 5 MB inline-payload ceiling
import os
size_mb = os.path.getsize("income_choropleth.html") / 1e6
print(f"Output: {size_mb:.2f} MB")   # Output: 1.84 MB
assert size_mb < 5, "Too large for inline GeoJSON — switch to vector tiles"

Roughly equal bin counts confirm the quantile classifier worked; a flat-grey map means the bin lookup failed.

Two more assertions are worth keeping in the script, because both failures are invisible in the rendered map. The first checks that the class breaks are strictly increasing — duplicated breaks happen when more than a fifth of the values are identical (a column padded with zeros, a capped survey response), and they collapse two classes into one shade without any warning. The second checks that every feature actually received a bin, catching the case where a row was reintroduced after classification by a later join.

import numpy as np

bins = classifier.bins
assert np.all(np.diff(bins) > 0), f"Duplicate class break — too many tied values: {bins}"

assert census_tracts["bin"].notna().all(), "Some tracts never got a bin"
print(f"{(census_tracts['bin'] < 0).sum()} tracts shaded as no-data")   # 7 tracts shaded as no-data

# Colour is baked at save time, so grep the output to prove the ramp reached the file
html = open("income_choropleth.html", encoding="utf-8").read()
print({c: html.count(c) for c in palette})
# {'#f0ebd8': 104, '#9bb3c9': 103, '#748cab': 104, '#3e5c76': 103, '#1d2d44': 622}

The last check exploits the fact that the style function runs in Python: every fill colour appears literally in the saved HTML, so counting occurrences proves the whole ramp was used. The darkest hex appears far more often than 104 times here because it is also the polygon stroke colour — subtract the stroke count, or use a stroke colour that is not in the palette when you want the check to be exact.

Edge Cases & Debugging

• All polygons one color. The style_function is reading the wrong property name, or you used the legacy Choropleth with a broken key_on. Print census_tracts.columns and confirm bin exists.
• NaN values break classification. mapclassify cannot bin nulls. Drop or impute first: census_tracts = census_tracts.dropna(subset=["median_income"]).
• Map lands in the ocean. You passed [lon, lat] to folium.Map; Folium wants [lat, lon] — reverse the centroid coordinates.
• Legend missing. You forgot colormap.add_to(income_map), or branca is outdated.
• Skewed, unreadable ramp. Switch Quantiles to mapclassify.NaturalBreaks for clustered data, or mapclassify.UserDefined for domain thresholds.
• File too big to share. Simplify geometry (in a projected CRS) and trim precision; past 5 MB, move to Generating PMTiles from GeoParquet.
• AttributeError: 'GeoSeries' object has no attribute 'set_precision'. GeoSeries.set_precision landed in GeoPandas 1.0. On 0.14 call the Shapely 2 function directly: shapely.set_precision(gdf.geometry.to_numpy(), 1e-6).
• The legend's numbers do not match the classes. StepColormap takes an index of n+1 boundaries for n colours. Passing classifier.bins alone gives it n values, so the ramp silently shifts by one class; prepend the data minimum.
• Two classes render as the same shade. Tied values collapsed a break — mapclassify returns duplicate bin edges when more than 1/k of the observations share a value. Drop to a smaller k, or switch to UserDefined breaks.
• The map is a picture of polygon size. You classified a count instead of a rate. Divide by population, households, or an equal-area-projected area before classifying.
• Everything shifts when you re-run on new data. Quantile breaks are recomputed from each dataset, so a monthly refresh silently re-scales the map and destroys comparability. Freeze the breaks with mapclassify.UserDefined(values, bins=FROZEN_BREAKS) for any map published as a time series.

Frequently Asked Questions

Quantiles or natural breaks — which should I default to? Default to quantiles when the map will be read comparatively, which is most of the time. Equal class counts mean each shade carries the same amount of data, so "in the darkest class" is a statement about rank that a reader can trust across the whole map. Natural breaks fits the distribution better by construction and is the right choice when the data has real gaps you want visible — two distinct populations of parcels, a bimodal elevation histogram. The failure mode of natural breaks is a top class holding a handful of features, which makes the extreme shade look far more common than it is; check .counts before you publish, not after.

How many classes should I use? Five is the working default and seven is a practical ceiling for a sequential ramp — beyond that most readers cannot match a polygon to a legend swatch, especially at the light end where adjacent steps differ by a few percent of lightness. Use fewer, not more, if the map is small on the page or destined for a report thumbnail. If the story genuinely needs fine gradation, the answer is a continuous ramp with an interactive tooltip carrying the exact value, not eleven classes.

Should I just use GeoDataFrame.explore() instead of all this? For exploration, yes — it is one line, it wires up mapclassify and a legend for you, and it returns a Folium Map you can add further layers to. Move to the explicit style_function when you need something explore() does not expose: a distinct colour for missing data, frozen class breaks across releases, a bespoke tooltip, or control over what ends up in the payload. The two are not different technologies — explore() builds the same folium.GeoJson object underneath — so migrating is a matter of writing out what the convenience call was doing implicitly.

How do I keep the colours stable when I republish next month? Freeze the breaks and store them next to the data. Quantile and natural-breaks classifiers derive their bin edges from whatever values you hand them, so a dataset refresh re-scales the ramp and a tract can change colour without changing value. Compute the breaks once, save them (a JSON file, a config constant), and classify every subsequent release with mapclassify.UserDefined(values, bins=saved_breaks). The trade-off is that a genuine shift in the distribution will start piling features into the end classes — which is a signal worth seeing, not a bug.

Sequential, diverging, or something else? Sequential (one hue, light to dark) for a quantity with a natural low-to-high reading: income, density, rainfall. Diverging (two hues meeting at a neutral midpoint) only when there is a meaningful centre — change from a baseline, deviation from a target, a net gain or loss — and the midpoint of the ramp must be pinned to that centre, not to the data median, or the map invents a story. Avoid rainbow ramps entirely: they are not perceptually ordered, so readers cannot rank two shades without consulting the legend, and they collapse badly for readers with colour-vision deficiency.

Can I put a choropleth of a million polygons in a Folium map? No, and the constraint is the payload rather than the classification. Every polygon is serialised into the HTML file, so the practical limit is the same few megabytes that bounds any inline-GeoJSON map. Simplify and trim first; if the layer is still too heavy, the value should be baked into vector tiles and shaded client-side, which is the path described in Styling Vector Tiles with Data-Driven Expressions under MapLibre GL Vector Web Maps.