Styling Vector Tiles with Data-Driven Expressions

A vector tile ships geometry and attributes but no colours, so every visual decision — the choropleth ramp, the category palette, how wide a line grows between zoom 8 and zoom 15 — is made in the style document by an expression, a small JSON array evaluated per feature on the GPU. This guide is for Python practitioners who compute the classification in Python and need to emit those arrays correctly. It sits under MapLibre GL Vector Web Maps in Web Mapping & Interactive Visualization, and picks up where Serving GeoJSON to MapLibre GL JS stops.

Why This Approach / What Goes Wrong

An expression is data, not code. It is a nested JSON array whose first element names an operator and whose remaining elements are the arguments — ["get", "median_rent"] reads a property off the current feature, ["to-number", …] coerces it, ["interpolate", …] blends between stops. Because it is plain JSON, Python can build it with a list comprehension and json.dump it, and because it is evaluated in the renderer, changing a break re-paints a million features without touching the tiles.

Four operators cover almost every real map. interpolate produces a continuous ramp: given ascending input/output stop pairs, it blends linearly (or exponentially) between them, so two features one dollar apart differ by one shade. step produces classes: it emits a fixed output until the input crosses the next boundary, which is what a choropleth actually is. match maps discrete labels to outputs with a mandatory fallback, for zoning codes or land-use classes. And ["zoom"] may appear as the input of a top-level interpolate or step — and only there — which is how a circle radius or line width scales as the user zooms in.

What goes wrong is rarely the syntax. It is that MapLibre has no classification primitive — no quantiles, no Jenks, no equal interval. It will faithfully render whatever boundaries you hand it, which means the statistics belong in Python and the expression is only their serialization. Two failure modes follow from that split. The first is drift: the ramp in the style and the numbers printed in the HTML legend are edited by hand in two places, and after the third data refresh they disagree. The second is typing. Expressions are strictly typed at evaluation time, and a median_rent that arrives from the tiles as the string "1240" rather than the number 1240 makes interpolate fail on every feature — no exception, no crash, just the paint property's default value across the whole layer.

Prerequisites

• geopandas>=1.0 — load the layer, reproject, and export the tiling source
• mapclassify>=2.6 — Jenks/natural breaks, quantiles, and the other classifiers behind the breaks
• numpy>=1.26 — array maths and the verification checks
• pandas>=2.0 — to_numeric for the attribute coercion
• MapLibre GL JS 4.x from a CDN, plus pmtiles@3 if the source is a PMTiles archive (nothing to install in Python)

conda install -c conda-forge "geopandas=1.0.*" "mapclassify=2.6.*" "numpy=1.26.*" "pandas=2.2.*"

Step-by-Step Implementation

The worked example styles a residential parcel layer by median_rent, served as vector tiles built with Generating PMTiles from GeoParquet.

1. Load the layer and fix the attribute types before anything else.

Shapefile DBF fields, CSV joins, and untyped GeoJSON all hand you numbers as object dtype. Whatever type the column has when you write the tiling source is the type the browser will see, so the coercion happens here, not in the style.

import geopandas as gpd
import pandas as pd

parcels = gpd.read_file("parcels.gpkg")[["parcel_id", "median_rent", "zoning", "geometry"]]
print(parcels["median_rent"].dtype)          # object  <- the bug, three steps upstream of the browser

parcels["median_rent"] = pd.to_numeric(parcels["median_rent"], errors="coerce")
parcels = parcels.dropna(subset=["median_rent"])

# Tiles are always WGS 84 lon/lat; any metric work (areas, buffers) must already be done
# in a projected CRS — here NAD83 / UTM 10N — never in Web Mercator.
parcels = parcels.to_crs("EPSG:4326")

2. Compute the classification breaks in Python.

mapclassify returns bins as an array of k upper bounds, the last of which is the maximum of the data. The interior boundaries — everything except that last value — are what the step expression needs.

import mapclassify
import numpy as np

rent = parcels["median_rent"].to_numpy(dtype="float64")
classifier = mapclassify.NaturalBreaks(rent, k=5)

bins = classifier.bins                    # k upper bounds; bins[-1] == rent.max()
print(np.round(bins))                     # [1180. 1620. 2130. 2890. 4350.]

3. Emit the step expression from those breaks.

Five colours need four interior boundaries. The first colour is the default output — everything below the first boundary — and each subsequent pair adds "at or above this value, use that colour". Cast every boundary with float(): json refuses to serialize a numpy.float64.

RAMP = ["#e3f3f7", "#58b7cf", "#0d6f87", "#134a5c", "#1d2d44"]

def step_expression(prop, bins, colors):
    """MapLibre step expression: len(colors) outputs, len(colors) - 1 boundaries."""
    if len(colors) != len(bins):
        raise ValueError(f"{len(colors)} colours cannot express {len(bins)} classes")
    expr = ["step", ["to-number", ["get", prop]], colors[0]]
    for boundary, color in zip(bins[:-1], colors[1:]):
        expr += [float(boundary), color]     # float() — np.float64 is not JSON serializable
    return expr

fill_color = step_expression("median_rent", bins, RAMP)
print(fill_color[:5])
# ['step', ['to-number', ['get', 'median_rent']], '#e3f3f7', 1180.0, '#58b7cf']

4. Use interpolate when the quantity is genuinely continuous, match when it is categorical.

step and interpolate read the same numbers and mean different things: step says "these five classes", interpolate says "this quantity, shaded". Choose by whether a reader will be asked to name the class a feature belongs to.

def interpolate_expression(prop, stops, colors):
    expr = ["interpolate", ["linear"], ["to-number", ["get", prop]]]
    for value, color in zip(stops, colors):
        expr += [float(value), color]        # inputs must strictly ascend
    return expr

# Even-percentile anchors keep the ramp readable on a skewed distribution
ramp_stops = np.quantile(rent, np.linspace(0, 1, len(RAMP)))
continuous_fill = interpolate_expression("median_rent", ramp_stops, RAMP)

ZONING = {"R1": "#e3f3f7", "R2": "#58b7cf", "C1": "#0d6f87", "M1": "#1d2d44"}

def match_expression(prop, mapping, fallback="#c3d0e4"):
    expr = ["match", ["get", prop]]
    for label, color in mapping.items():
        expr += [label, color]               # labels must all be the same type
    expr.append(fallback)                    # mandatory, and it catches nulls
    return expr

zoning_fill = match_expression("zoning", ZONING)

5. Scale sizes with zoom, and nest the data expression inside it.

["zoom"] is only legal as the input of a top-level interpolate or step, but the outputs of that top-level expression may themselves be data-driven. That nesting is how a symbol stays legible at zoom 9 and still encodes its attribute at zoom 15.

def zoom_scaled_radius(prop, low_zoom_range, high_zoom_range, domain):
    lo, hi = domain
    return [
        "interpolate", ["exponential", 1.4], ["zoom"],
        9,  ["interpolate", ["linear"], ["to-number", ["get", prop]],
             lo, low_zoom_range[0], hi, low_zoom_range[1]],
        15, ["interpolate", ["linear"], ["to-number", ["get", prop]],
             lo, high_zoom_range[0], hi, high_zoom_range[1]],
    ]

circle_radius = zoom_scaled_radius("unit_count", (2, 9), (5, 26), domain=(1, 400))

6. Write the style and the legend from the same break array.

This is the whole point of generating expressions in Python: the colours and the labels come out of one loop, so they cannot drift apart. Note the interval convention — a step boundary is the lower, inclusive edge of its class, so the labels read < b0, b0 – b1, …, ≥ b_last.

import json

edges = [float(b) for b in bins[:-1]]      # the four interior boundaries
labels = (
    [f"< {edges[0]:,.0f}"]
    + [f"{lo:,.0f} – {hi:,.0f}" for lo, hi in zip(edges, edges[1:])]
    + [f"≥ {edges[-1]:,.0f}"]
)

style = {
    "version": 8,
    "sources": {"parcels": {"type": "vector", "url": "pmtiles://parcels.pmtiles"}},
    "layers": [{
        "id": "parcel-rent",
        "type": "fill",
        "source": "parcels",
        "source-layer": "parcels",          # must match the layer name inside the tiles
        "paint": {
            "fill-color": fill_color,
            "fill-opacity": 0.85,
            "fill-outline-color": "#1d2d44",
        },
    }],
}
legend = {
    "title": "Median rent (USD/month)",
    "entries": [{"color": c, "label": l} for c, l in zip(RAMP, labels)],
}

with open("public/parcels-style.json", "w", encoding="utf-8") as fh:
    json.dump(style, fh, indent=2)
with open("public/legend.json", "w", encoding="utf-8") as fh:
    json.dump(legend, fh, ensure_ascii=False, indent=2)

7. Point the map at the generated style and build the legend from its sibling file.

page = """<!doctype html>
<html><head><meta charset="utf-8">
<script src="https://unpkg.com/maplibre-gl@4.7.1/dist/maplibre-gl.js"></script>
<script src="https://unpkg.com/pmtiles@3/dist/pmtiles.js"></script>
<link href="https://unpkg.com/maplibre-gl@4.7.1/dist/maplibre-gl.css" rel="stylesheet"/>
<style>html,body,#map{height:100%;margin:0}
#legend{position:absolute;bottom:20px;left:20px;background:#fff;padding:10px;font:12px sans-serif}
#legend i{display:inline-block;width:14px;height:14px;margin-right:6px;vertical-align:-2px}
</style></head>
<body><div id="map"></div><div id="legend"></div><script>
maplibregl.addProtocol("pmtiles", new pmtiles.Protocol().tile);   // before the Map is built
const map = new maplibregl.Map({
  container: "map", style: "parcels-style.json",
  center: [-122.335, 47.61], zoom: 11,           // lon, lat
});
fetch("legend.json").then(r => r.json()).then(legend => {
  document.getElementById("legend").innerHTML =
    "<b>" + legend.title + "</b>" +
    legend.entries.map(e => '<div><i style="background:' + e.color + '"></i>' + e.label + "</div>").join("");
});
</script></body></html>"""

with open("public/index.html", "w", encoding="utf-8") as fh:
    fh.write(page)

Verification

Check the expression's shape and its agreement with the classifier before opening a browser — a wrong stop order or a stray numpy scalar is far cheaper to catch in Python.

import numpy as np

# 1. It must survive a JSON round trip (this is where np.float64 raises TypeError)
assert json.loads(json.dumps(fill_color)) == fill_color

# 2. Structure: one more output than boundary, boundaries strictly ascending
boundaries = np.asarray(fill_color[3::2], dtype="float64")
outputs = fill_color[2::2]
assert len(outputs) == len(boundaries) + 1 == len(RAMP)
assert (np.diff(boundaries) > 0).all(), "step inputs must strictly ascend"

# 3. The renderer's class for each parcel must match mapclassify's
maplibre_class = np.searchsorted(boundaries, rent, side="right")
on_a_break = np.isin(rent, boundaries)
assert (maplibre_class[~on_a_break] == classifier.yb[~on_a_break]).all()
print(f"{on_a_break.sum()} parcels sit exactly on a break")   # 2 parcels sit exactly on a break

# 4. The legend cannot disagree with the ramp
assert [e["color"] for e in legend["entries"]] == outputs
print(f"{len(outputs)} classes, {len(legend['entries'])} legend rows")  # 5 classes, 5 legend rows

Check 3 deliberately excludes values sitting exactly on a boundary, because the two tools disagree there by design: mapclassify uses half-open intervals closed on the upper bound, while step closes on the lower bound. A parcel at exactly 1180 is class 0 in Python and class 1 in the browser. The count printed above tells you whether that matters for your data.

Edge Cases & Debugging

• Every feature is one flat colour. The styled property is a string in the tiles. Confirm in the console with map.querySourceFeatures("parcels", {sourceLayer: "parcels"})[0].properties, then fix the dtype in Python and re-tile — ["to-number", …] rescues most cases, but not text like "1,240" or "N/A".
• Expected value to be of type number, but found string instead on addLayer. A literal in your stop pairs is quoted. This one is caught by style validation at load time, unlike the runtime typing failure above.
• Input/output pairs for "step" expressions must be arranged with input values in strictly ascending order. Duplicate breaks, usually from Quantiles on a distribution with heavy ties. Drop the repeats, use fewer classes, or switch to NaturalBreaks.
• match paints everything the fallback colour. The labels and the attribute have different types — ["match", ["get", "zone_code"], 1, …] never matches the string "1". Make the labels match the attribute exactly; MapLibre also requires all labels in one match to share a type.
• Null attributes land in the lowest class. ["to-number", …] converts null to 0. Wrap the ramp in a presence test to paint no-data grey instead: ["case", ["has", "median_rent"], fill_color, "#c3d0e4"].
• TypeError: Object of type float64 is not JSON serializable. A break went into the array straight from a NumPy array; wrap it in float() as in step 3.

Frequently Asked Questions

Can MapLibre compute the classification itself? No. The style specification has no quantile, Jenks, or equal-interval operator — expressions evaluate against a single feature at a time and have no view of the distribution. Every boundary must be computed upstream, which is why mapclassify is a hard dependency of this workflow rather than a convenience. The upside is reproducibility: the same breaks used by the web map can be reused for a static report or a Folium choropleth without a second implementation.

Do expressions behave the same on a GeoJSON source and on vector tiles? The evaluation is identical, but the data reaching them is not. A GeoJSON source preserves JSON types exactly, so a number stays a number; a tiling step can widen, narrow, or stringify an attribute, and it can drop attributes entirely if you filtered them at build time. Vector-tile layers also require source-layer, which a GeoJSON layer must omit. If a ramp works against GeoJSON and fails against tiles, suspect the tiling step first — tippecanoe -T median_rent:float forces the type explicitly.

How do I change the ramp without rebuilding the tiles? Call map.setPaintProperty("parcel-rent", "fill-color", newExpression) with a fresh array. Styling is fully decoupled from the data, so re-classifying is a JSON swap rather than a re-tile — which is exactly what makes an interactive break selector cheap to build. Regenerate legend.json from the same bins array in the same call, or the two drift apart within one deploy.

Should the ramp use step or interpolate? Use step whenever readers need to name the class a feature belongs to, or when the underlying quantity is ordinal — that is nearly every choropleth. Use interpolate for genuinely continuous fields where the shape of the surface matters more than any threshold, such as elevation, temperature, or a density estimate. interpolate on a skewed variable hides most of the range in one shade; if you reach for it there, use ["exponential", 1.5] or classify instead.