How to Style Data-Driven Colours in a Web Map

Problem statement

Data-driven styling is what makes vector tiles worth using: the same tiles serve every style, and the colour is computed in the browser from the feature's attributes.

Two things make it harder than it looks.

The attributes have to be in the tiles. Styling by a field that was dropped during generalisation produces the fallback colour everywhere, with no error. Attributes cost real bytes — two string fields added 21% to a measured tile of 401 buildings, 14% after gzip.

The classification decides the map. Equal-interval, quantile and natural-breaks classifications of the same data produce different maps and support different conclusions, and the style sheet is where that choice gets made without anyone noticing.

Quick answer

paint: {
  // continuous: interpolate between stops
  "fill-color": [
    "interpolate", ["linear"], ["get", "density"],
    0,    "#f7fbff",
    50,   "#c6dbef",
    200,  "#6baed6",
    1000, "#08306b",
  ],
  // categorical: match exact values, with a fallback
  "fill-outline-color": [
    "match", ["get", "class"],
    "residential", "#1a3a6b",
    "commercial",  "#0ea5e9",
    /* fallback */ "#94a3b8",
  ],
  // zoom-dependent
  "fill-opacity": ["interpolate", ["linear"], ["zoom"], 12, 0.2, 16, 0.7],
}

match requires a fallback as its last argument. Without one, features whose value is not listed get no paint and are invisible.

buildings ['area_m2', 'building']

Attributes in the tiles, value types, a classification

1. Make sure the attribute is in the tiles

from pmtiles.reader import Reader, MmapSource
with open("tiles.pmtiles", "rb") as f:
    for layer in Reader(MmapSource(f)).metadata()["vector_layers"]:
        print(layer["id"], sorted(layer.get("fields", {})))

buildings ['area_m2', 'building']

A style referring to a missing field is not an error — ["get", "density"] returns null, the expression falls through to its default, and the map renders in one flat colour.

2. Watch the type of the value

MVT properties are typed, and tile generators frequently write everything as strings.

["interpolate", ["linear"], ["to-number", ["get", "density"]], ...]

interpolate on a string produces no output and no error. ["to-number", ...] is the defensive form, and returns null for values it cannot parse — which the fallback then catches.

3. Choose the classification deliberately

The same values classified three ways give three maps:

• Equal interval — even value ranges. Honest about magnitude, and useless when the distribution is skewed, which spatial data usually is.
• Quantile — even feature counts per class. Always uses the full palette, and can put nearly identical values in different classes.
• Natural breaks — minimises within-class variance. Usually the most readable and the hardest to explain.

Compute the breaks in Python, where you can inspect them, and paste them into the style rather than letting the style imply them.

4. Use step for classes, interpolate for continua

["step", ["get", "value"], "#f7fbff", 50, "#c6dbef", 200, "#6baed6"]

step produces discrete classes with sharp boundaries — the choropleth convention, and easier to read from a legend. interpolate produces a continuous ramp, which suits density surfaces and looks vague on polygons.

5. Handle nulls explicitly

["case",
  ["==", ["get", "value"], null], "#e2e8f0",
  ["interpolate", ["linear"], ["get", "value"], 0, "#f7fbff", 100, "#08306b"]]

"No data" and "zero" are different, and a colour ramp that renders them identically is a factual error on the map.

Generating the style expression from the data

The data used on this page

Every building polygon in OpenStreetMap inside −2.31 to −2.18° E, 53.44 to 53.52° N — central Manchester and Salford, centred on the Town Hall — fetched with OSMnx. The setup keeps the building type, name, building:levels and addr:street tags and a footprint area in square metres (British National Grid), and writes buildings.parquet. It then builds tiles.pmtiles from it with GDAL's PMTiles driver: one layer, buildings, zooms 10 to 15, carrying building and area_m2. The buildings are read as Overpass returned them on 9 October 2026.

import io
import os
import requests
import time
import zipfile

import geopandas as gpd
import osmnx as ox

# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-manchester-buildings-2026-10-09.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.requests_timeout = 600
BOX = (-2.31, 53.44, -2.18, 53.52)                   # west, south, east, north

if not os.path.exists("buildings.parquet"):
    for attempt in range(5):                         # the Overpass mirror is sometimes busy
        try:
            raw = ox.features_from_bbox(BOX, tags={"building": True})
            break
        except Exception:
            time.sleep(30)
    raw = raw[raw.geom_type.isin(["Polygon", "MultiPolygon"])].reset_index()
    for tag in ("name", "building:levels", "addr:street"):
        if tag not in raw:
            raw[tag] = None
    b = raw[["id", "building", "name", "building:levels", "addr:street", "geometry"]].rename(
        columns={"building:levels": "levels", "addr:street": "street"})
    b["area_m2"] = b.to_crs(27700).area.round(1)
    b.to_parquet("buildings.parquet")

if not os.path.exists("tiles.pmtiles"):
    gpd.read_parquet("buildings.parquet")[["building", "area_m2", "geometry"]].to_file(
        "tiles.pmtiles", driver="PMTiles", layer="buildings",
        dataset_options={"MINZOOM": 10, "MAXZOOM": 15})

buildings = gpd.read_parquet("buildings.parquet")
print(f"{len(buildings):,} buildings; tiles.pmtiles {os.path.getsize('tiles.pmtiles') / 1e6:.2f} MB")

63,517 buildings; tiles.pmtiles 5.18 MB

Example 1 — generating the style expression from the data

import json
import numpy as np
import mapclassify


def colour_expression(values, field, scheme="quantiles", k=5,
                      palette=("#f7fbff", "#c6dbef", "#6baed6",
                               "#2171b5", "#08306b"), null_colour="#e2e8f0"):
    """A MapLibre step expression with breaks computed and reported."""
    clean = np.asarray(values, dtype=float)
    clean = clean[np.isfinite(clean)]

    classifier = {
        "quantiles": mapclassify.Quantiles,
        "equalinterval": mapclassify.EqualInterval,
        "naturalbreaks": mapclassify.NaturalBreaks,
    }[scheme](clean, k=k)
    breaks = [round(float(b), 2) for b in classifier.bins[:-1]]

    print(f"  {scheme}, k={k}: breaks {[round(b, 2) for b in breaks]}")
    for i, count in enumerate(classifier.counts):
        print(f"    class {i}: {count:,} features "
              f"({count / len(clean):.1%})")

    expression = ["case",
                  ["==", ["get", field], None], null_colour,
                  ["step", ["to-number", ["get", field]], palette[0]]]
    for break_value, colour in zip(breaks, palette[1:]):
        expression[-1].extend([break_value, colour])

    return expression


areas = buildings.area_m2
print(f"{len(areas):,} footprints, median {areas.median():,.1f} m², largest {areas.max():,.0f} m²")
quantile_expr = colour_expression(areas, "area_m2", "quantiles")
equal_expr = colour_expression(areas, "area_m2", "equalinterval")
print(json.dumps(quantile_expr))

63,517 footprints, median 58.1 m², largest 54,804 m²
  quantiles, k=5: breaks [45.7, 53.1, 66.5, 132.4]
    class 0: 12,933 features (20.4%)
    class 1: 12,592 features (19.8%)
    class 2: 12,592 features (19.8%)
    class 3: 12,703 features (20.0%)
    class 4: 12,697 features (20.0%)
  equalinterval, k=5: breaks [10962.44, 21922.78, 32883.12, 43843.46]
    class 0: 63,473 features (99.9%)
    class 1: 37 features (0.1%)
    class 2: 6 features (0.0%)
    class 3: 0 features (0.0%)
    class 4: 1 features (0.0%)
["case", ["==", ["get", "area_m2"], null], "#e2e8f0", ["step", ["to-number", ["get", "area_m2"]], "#f7fbff", 45.7, "#c6dbef", 53.1, "#6baed6", 66.5, "#2171b5", 132.4, "#08306b"]]

Printing the class counts is what turns the classification into a decision. Equal interval on the same footprints put 99.9% of buildings in the first class, because a handful of warehouses and shopping centres stretch the range to 54,804 m²; the counts make that obvious immediately.

Example 2 — checking a style against the tiles

import gzip
import re
import mapbox_vector_tile as mvt
from pmtiles.reader import Reader, MmapSource


def validate_style(style, archive_path, sample_tile):
    """Do the fields the style uses exist in the tiles?"""
    text = json.dumps(style)
    used = set(re.findall(r'\["get",\s*"([^"]+)"\]', text))
    source_layers = set(re.findall(r'"source-layer":\s*"([^"]+)"', text))

    with open(archive_path, "rb") as handle:
        reader = Reader(MmapSource(handle))
        data = reader.get(*sample_tile)                  # z, x, y
        decoded = mvt.decode(gzip.decompress(data))

    available_layers = set(decoded)
    print(f"  style uses layers {sorted(source_layers)}")
    print(f"  tile contains     {sorted(available_layers)}")
    for missing in source_layers - available_layers:
        print(f"  ! source-layer '{missing}' is not in the tiles")

    fields = set()
    for layer in decoded.values():
        for feature in layer["features"][:50]:
            fields.update(feature.get("properties", {}))
    print(f"  style reads fields {sorted(used)}")
    for missing in used - fields:
        print(f"  ! field '{missing}' is not in the tiles — "
              "the expression will fall through to its default")
    return used - fields


import mercantile

town_hall = mercantile.tile(-2.2445, 53.4794, 14)
style = {"layers": [{"id": "footprints", "type": "fill", "source": "osm",
                     "source-layer": "buildings",
                     "paint": {"fill-color": quantile_expr,
                               "fill-outline-color": ["match", ["get", "class"],
                                                      "residential", "#1a3a6b", "#94a3b8"]}}]}
missing = validate_style(style, "tiles.pmtiles", (town_hall.z, town_hall.x, town_hall.y))

  style uses layers ['buildings']
  tile contains     ['buildings']
  style reads fields ['area_m2', 'class']
  ! field 'class' is not in the tiles — the expression will fall through to its default

Both failures — a missing layer and a missing field — are silent in the browser. Checking them against a real tile before deploying is a minute's work.

Example 3 — a legend that matches the style

def legend_from_expression(expression, field, unit=""):
    """Extract the breaks and colours so the legend cannot drift."""
    step = expression
    while isinstance(step, list) and step[0] != "step":
        step = next((part for part in step
                     if isinstance(part, list) and part and part[0] == "step"),
                    None)
        if step is None:
            raise ValueError("no step expression found")

    default_colour = step[2]
    pairs = step[3:]
    entries = [{"label": f"< {pairs[0]}{unit}", "colour": default_colour}]
    for i in range(0, len(pairs), 2):
        value, colour = pairs[i], pairs[i + 1]
        upper = pairs[i + 2] if i + 2 < len(pairs) else None
        entries.append({
            "label": f"{value}{unit} – {upper}{unit}" if upper
                     else f"≥ {value}{unit}",
            "colour": colour,
        })

    for entry in entries:
        print(f"  {entry['colour']}  {entry['label']}")
    return entries


legend = legend_from_expression(quantile_expr, "area_m2", " m²")

  #f7fbff  < 45.7 m²
  #c6dbef  45.7 m² – 53.1 m²
  #6baed6  53.1 m² – 66.5 m²
  #2171b5  66.5 m² – 132.4 m²
  #08306b  ≥ 132.4 m²

Deriving the legend from the style expression rather than writing it separately is the only way to stop the two drifting apart. A legend that does not match the map is worse than no legend.

Example 4 — what two string attributes cost in a tile

import gzip
import mapbox_vector_tile as mvt
from shapely.ops import transform

tile = mercantile.tile(-2.2445, 53.4794, 15)
west, south, east, north = mercantile.xy_bounds(tile)
in_tile = buildings.to_crs(3857).cx[west:east, south:north].copy()
in_tile["geometry"] = in_tile.clip((west, south, east, north)).geometry


def tile_bytes(props):
    features = [{"geometry": g, "properties": {k: row[k] for k in props if isinstance(row[k], str)}}
                for g, (_, row) in zip(in_tile.geometry, in_tile.iterrows()) if not g.is_empty]
    raw = mvt.encode([{"name": "buildings", "features": features}],
                     default_options={"quantize_bounds": (west, south, east, north)})
    return len(raw), len(gzip.compress(raw))


plain, plain_gz = tile_bytes(())
tagged, tagged_gz = tile_bytes(("building", "street"))
print(f"z{tile.z}/{tile.x}/{tile.y}: {len(in_tile)} buildings")
print(f"  geometry only          {plain / 1000:6.2f} kB, gzipped {plain_gz / 1000:5.2f} kB")
print(f"  + building and street  {tagged / 1000:6.2f} kB, gzipped {tagged_gz / 1000:5.2f} kB  "
      f"(+{tagged / plain - 1:.0%}, +{tagged_gz / plain_gz - 1:.0%} gzipped)")

z15/16179/10601: 401 buildings
  geometry only           12.64 kB, gzipped  9.53 kB
  + building and street   15.25 kB, gzipped 10.82 kB  (+21%, +14% gzipped)

Why a missing field is silent

Why a missing field is silent

MVT features carry only the properties the encoder wrote. ["get", "density"] on a feature without that property returns null.

interpolate and step given null fall through to their default, which is usually the first colour. So the map renders — in one flat colour, uniformly, with no console message.

The same thing happens when generalisation drops attributes at low zoom to save bytes: the map is correctly styled at high zoom and flat at low zoom, which looks like a data problem and is a pipeline one.

Why types matter more than in Python

MVT properties are typed as string, number or boolean, and many generators write numbers as strings — often because the source column was object dtype.

MapLibre's interpolate requires a number. Given a string it produces no output, and again there is no error.

["to-number", ["get", "field"]] coerces and returns null on failure, which a case can catch. Alternatively, fix the types when writing the tiles, which is better and requires controlling the generator.

Why the classification is the map

A choropleth's message comes from its classification. Equal interval on a skewed distribution puts almost everything in the first class and shows a nearly uniform map. Quantiles guarantee an even spread of colour, which can exaggerate differences between nearly identical values.

Neither is wrong; they answer different questions. The problem is that the choice usually happens inside a style sheet, invisibly, and the reader has no way to know which was used.

Computing the breaks in Python and printing the class counts makes it a decision. Recording the scheme in the legend makes it honest.

Why null must be distinguished from zero

A colour ramp maps zero to its lowest colour. If missing data is also rendered as the lowest colour — or as the fallback, which is often the same — the map asserts that unmeasured places have low values.

For a density map that reads as "nobody lives here" where the truth is "we did not measure". The convention is a distinct neutral grey with its own legend entry, and it costs one case branch.

Defaults, string numbers and match fallbacks

• A missing field falls through to the default with no error.
• Numbers can arrive as strings. Wrap in ["to-number", ...].
• match needs a fallback as its final argument.
• step for classes, interpolate for continua.
• Null and zero are different. Give null its own colour and legend entry.
• Compute breaks in Python and paste them in; print the class counts.
• Attributes cost tile bytes — two string fields added 21% to a measured tile, 14% gzipped.
• Derive the legend from the expression so they cannot drift apart.

Internal links

• Choropleth classification explained: quantiles, equal interval and natural breaks — choosing the scheme
• Vector tiles explained: MVT, layers and why they are not GeoJSON — what attributes cost
• How to serve a web map from PMTiles with MapLibre — where the style lives
• My web map is blank or the layer never appears — when styling hides everything
• How to add interactivity and popups to a web map — the other use of attributes
• My GeoPandas choropleth colours look wrong — the same problem in Matplotlib
• How to build vector tiles from a GeoDataFrame in Python — choosing which attributes to include
• How to add legends, labels and a scale bar to a GeoPandas map — legend conventions

FAQ

Why is my whole layer one colour?

The styled field is not in the tiles, or arrives as a string. The expression falls through to its default, silently.

How do I check which fields are in my tiles?

Read vector_layers from the archive metadata, or decode one tile and print the properties of a few features.

Why does interpolate not work on my numeric field?

It is probably a string. Wrap it in ["to-number", ["get", "field"]], or fix the type when generating the tiles.

Should I use step or interpolate?

step for classed choropleths with a discrete legend, interpolate for continuous surfaces such as density.

Where should I compute the class breaks?

In Python, where you can print the class counts and compare schemes. Paste the resulting breaks into the style.

How do I show missing data?

A case branch testing for null, rendering a neutral grey, with its own legend entry. Never let null share a colour with zero.

Does adding attributes make tiles much bigger?

Measurably. Two string attributes, building and street, added 21% to a tile of 401 buildings in central Manchester, and 14% after gzip. Include only what the style and interaction need.