Generating PMTiles from GeoParquet

This guide turns a GeoParquet dataset into a single .pmtiles file you can drop on any static bucket and render in MapLibre GL — no tile server, no database. It is for anyone whose GeoJSON has grown too large to inline and who wants cloud-native, server-less vector tiles. It sits under Vector Tile Pipelines with PMTiles in Web Mapping & Interactive Visualization; if you instead need a running endpoint for auth or on-the-fly filtering, see the sibling recipe Serving MBTiles with Python.

Why This Approach / What Goes Wrong

PMTiles solves the "I have a million features and no budget for a tile server" problem. A naive GeoJSON of that size hangs the browser; a traditional tile server adds a process, a database, and operational cost. PMTiles instead packs the whole tile pyramid into one flat file with an internal directory, and clients fetch only the byte ranges for the tiles currently in view over HTTP Range requests. Storage plus egress is the entire bill.

The build has two sharp edges that break most first attempts. The first is a format gap: tippecanoe, the tiler, ingests GeoJSON, line-delimited GeoJSON, or FlatGeobuf — not GeoParquet directly. You bridge that gap with ogr2ogr from GDAL, and if your GDAL build lacks the Arrow/Parquet driver the conversion fails before tiling even starts. The second edge is projection: tippecanoe expects EPSG:4326 (lon/lat) and performs the Web Mercator projection itself. Feed it a projected CRS — a UTM zone, a national grid, or Web Mercator EPSG:3857 — and every feature lands in the wrong place because tippecanoe reads your metres as degrees. Reprojecting to 4326 as the last step, after all metric work is done, is the single rule that keeps tiles aligned.

GeoParquet is the right source format precisely because it preserves CRS metadata, compresses well, and is the analytical store described in Cloud-Native Geospatial Formats. Treat tiles as a derived rendering view of that store — rebuild them from GeoParquet whenever the data changes rather than editing tiles in place.

Prerequisites

• gdal>=3.8 for ogr2ogr (the conda-forge build ships the Parquet driver)
• tippecanoe>=2.51 (system binary; no native Windows build — use WSL2 or a container)
• pmtiles>=3.2 Python package (MBTiles→PMTiles conversion plus a reader for verification)
• geopandas>=0.14 for source preparation and reprojection
• duckdb>=1.0 (optional) when the tiling source is a filtered or joined subset of a larger store

# Python-side stack from conda-forge keeps GDAL/GEOS/PROJ coherent:
conda install -c conda-forge "gdal=3.8.*" "pmtiles=3.2.*" "geopandas=0.14.*"
# tippecanoe is a native binary, installed separately:
#   macOS:  brew install tippecanoe
#   Linux:  git clone https://github.com/felt/tippecanoe && cd tippecanoe && make -j && sudo make install

Step-by-Step Implementation

1. Prepare the GeoParquet source in EPSG:4326. Do every metric computation — areas, lengths, buffers — before this reprojection, while the data is still in a projected CRS, then convert to 4326 as the final export. Reprojection uses the transformer stack behind Coordinate Systems with PyProj; the to_crs call here is the delivery hand-off, not an analysis step.

import geopandas as gpd

# building_footprints were analysed in EPSG:25832 (UTM 32N, metric) upstream
building_footprints = gpd.read_file("building_footprints.gpkg")
assert building_footprints.crs is not None, "unknown CRS — set it before reprojecting"

building_footprints["area_m2"] = building_footprints.geometry.area  # metric attribute FIRST, in UTM
building_footprints = building_footprints.to_crs(epsg=4326)          # reproject LAST, for delivery

assert building_footprints.crs.to_epsg() == 4326, "tippecanoe requires EPSG:4326"
building_footprints[["building_id", "area_m2", "geometry"]].to_parquet("buildings_4326.parquet")

2. Convert GeoParquet to line-delimited GeoJSON with ogr2ogr. GeoJSONSeq writes one feature per line, which tippecanoe streams instead of loading whole — the difference between a bounded and an unbounded memory footprint on large inputs.

import subprocess

subprocess.run([
    "ogr2ogr",
    "-f", "GeoJSONSeq",
    "buildings.geojsonl",       # newline-delimited GeoJSON, streamed by tippecanoe
    "buildings_4326.parquet",   # GeoParquet source
], check=True)

If the tiling source is a subset of a larger store — one city out of a national dataset, or a join against an attribute table — doing that work in ogr2ogr means materialising the whole thing first. DuckDB reads GeoParquet natively and writes GeoJSONSeq through the same GDAL drivers, so the filter, the join and the format conversion collapse into one query. The engine and its spatial extension are covered in Querying GeoParquet with DuckDB Spatial:

import duckdb

con = duckdb.connect()
con.execute("INSTALL spatial; LOAD spatial;")
con.execute("""
    COPY (
        SELECT b.building_id, b.area_m2, b.geometry
        FROM read_parquet('buildings/*.parquet') AS b
        JOIN read_parquet('permits.parquet') AS p USING (building_id)
        WHERE b.area_m2 > 40 AND p.status = 'issued'
    ) TO 'buildings.geojsonl' WITH (FORMAT GDAL, DRIVER 'GeoJSONSeq')
""")

The glob reads a partitioned GeoParquet directory without concatenating it, and only the rows that survive the WHERE clause ever become text. On a hive-partitioned store this is routinely an order of magnitude less I/O than converting the whole dataset and letting tippecanoe discard most of it.

3. Tile with tippecanoe. Let it pick the maximum zoom from feature density, and coalesce (merge) rather than drop features where a tile overflows its size budget — the right choice for contiguous polygons like building footprints.

subprocess.run([
    "tippecanoe",
    "-o", "buildings.mbtiles",
    "-zg",                           # guess the max zoom from feature density
    "--coalesce-densest-as-needed",  # merge features in crowded tiles, don't drop them
    "--extend-zooms-if-still-dropping",
    "-y", "building_id",             # whitelist the fields the style/tooltip need
    "-y", "area_m2",
    "-l", "buildings",               # layer name MapLibre's source-layer must match
    "--force",                       # overwrite an existing output
    "buildings.geojsonl",
], check=True)

Two flags change the build's cost rather than its output. -P reads the input in parallel across cores, and it works only on line-delimited input — the reason to prefer GeoJSONSeq over plain GeoJSON even when memory is not a concern. -t relocates the scratch space, which matters because tippecanoe writes intermediate per-zoom files that routinely exceed the input several times over; on a machine whose /tmp is a small tmpfs, a national dataset will fail hours into the build with a disk-full error and no partial output:

subprocess.run([
    "tippecanoe", "-o", "buildings.mbtiles", "-zg",
    "-P",                              # parallel read — requires line-delimited input
    "-t", "/mnt/scratch/tippecanoe",   # temp space; the build can need several × the input
    "--coalesce-densest-as-needed",
    "-y", "building_id", "-y", "area_m2",
    "-l", "buildings", "--force", "buildings.geojsonl",
], check=True)

As a rough envelope on a modern eight-core laptop with an SSD: a few hundred thousand building footprints to zoom 14 takes single-digit minutes and a few gigabytes of scratch; ten million takes an hour or more and tens of gigabytes. Wall-clock time scales with the deepest zoom far more steeply than with feature count, because each extra level roughly quadruples the tiles written — going from --maximum-zoom 14 to 16 is a sixteen-fold increase in tiles, not a 15% one. Decide the real maximum zoom from how far users actually zoom, and cap it explicitly once -zg has told you what it would have guessed.

4. Convert MBTiles to PMTiles. The pmtiles package rewrites the SQLite tile store into the single-file, range-readable archive.

subprocess.run(["pmtiles", "convert", "buildings.mbtiles", "buildings.pmtiles"], check=True)

If the map needs several layers — footprints, parcels, and a street centreline, each with its own zoom range and its own shedding flags — build one MBTiles per layer and merge them with tile-join before converting. That keeps each layer's tiling decisions independent and still produces a single archive for the client to range-read:

subprocess.run([
    "tile-join", "-o", "city.mbtiles", "--force",
    "buildings.mbtiles", "parcels.mbtiles", "streets.mbtiles",
], check=True)
subprocess.run(["pmtiles", "convert", "city.mbtiles", "city.pmtiles"], check=True)

tile-join is also the incremental-rebuild tool: re-tile only the layer whose source changed and re-merge, rather than re-running a multi-hour build for a dataset that was already correct. Each layer keeps the name it was given by -l, so the MapLibre style references them as three source-layer values against one source.

5. Reference it from MapLibre. Register the pmtiles:// protocol shim before constructing the map, then add a vector source whose source-layer equals tippecanoe's -l value. The full client-side pattern lives in MapLibre GL Vector Web Maps.

# In the page JS (not Python):
#   const protocol = new pmtiles.Protocol();
#   maplibregl.addProtocol("pmtiles", protocol.tile);   // register BEFORE new Map()
#   map.addSource("buildings", {
#     type: "vector",
#     url: "pmtiles://https://cdn.example.com/buildings.pmtiles",
#   });
#   map.addLayer({ id: "bld", type: "fill", source: "buildings",
#                  "source-layer": "buildings",          // must equal tippecanoe -l
#                  paint: { "fill-color": "#3e5c76", "fill-opacity": 0.7 } });

Verification

Inspect the PMTiles header and metadata to confirm the zoom range, bounds, and layer name before deploying — a five-line assertion turns a manual click-through into a CI guardrail.

from pmtiles.reader import Reader, MmapSource

with open("buildings.pmtiles", "rb") as fh:
    reader = Reader(MmapSource(fh))
    header = reader.header()
    meta = reader.metadata()

print("min/max zoom:", header["min_zoom"], header["max_zoom"])       # 0 14
print("bounds (lon):", header["min_lon_e7"] / 1e7, header["max_lon_e7"] / 1e7)
print("layers:", [layer["id"] for layer in meta["vector_layers"]])   # ['buildings']

assert header["max_zoom"] >= 12, "max zoom too low for street-level detail"
assert meta["vector_layers"][0]["id"] == "buildings", "layer name must match the style source-layer"

Bounds are stored as E7 integers (degrees × 10⁷); divide by 1e7 and confirm the box matches your dataset's real extent. A wildly wrong box — or coordinates near Null Island (0, 0) — means the source was not in EPSG:4326 before tiling.

Header checks prove the archive is well-formed; they do not prove any tile contains features. Pull one tile from the middle of the extent and assert it is non-empty, which catches the build that succeeded while dropping everything:

import math

def lonlat_to_tile(lon, lat, zoom):
    """Web Mercator XYZ tile containing a lon/lat, at a given zoom."""
    n = 2 ** zoom
    x = int((lon + 180.0) / 360.0 * n)
    lat_rad = math.radians(lat)
    y = int((1 - math.asinh(math.tan(lat_rad)) / math.pi) / 2 * n)
    return x, y

centre_lon = (header["min_lon_e7"] + header["max_lon_e7"]) / 2e7
centre_lat = (header["min_lat_e7"] + header["max_lat_e7"]) / 2e7
z = header["max_zoom"] - 2
x, y = lonlat_to_tile(centre_lon, centre_lat, z)

tile = reader.get(z, x, y)
assert tile, f"tile {z}/{x}/{y} is empty at the centre of the extent"
assert tile[:2] == b"\x1f\x8b", "tiles should be gzip-compressed inside the archive"
print(f"tile {z}/{x}/{y}: {len(tile):,} bytes")   # tile 12/2200/1343: 41,208 bytes

The gzip magic bytes are worth asserting because they decide how the archive must be served. Tiles are compressed inside the file, so the object store must hand them over untouched; a host that adds its own Content-Encoding: gzip on top produces tiles the browser decompresses once and the PMTiles shim then fails to decompress again. A tile size in the tens of kilobytes is healthy; one that consistently brushes 500 KB means a shedding flag is doing nothing and the renderer will pay for it on every pan.

Edge Cases & Debugging

• ogr2ogr can't read Parquet. Your GDAL build lacks the Arrow/Parquet driver; install gdal from conda-forge, which bundles it, and check with ogrinfo --formats | grep -i parquet.
• Everything is in the wrong place / near Null Island. The source was not reprojected to EPSG:4326, or its CRS was None. Assert crs.to_epsg() == 4326 immediately before invoking tippecanoe; see Coordinate Reference System Transformations.
• Nothing renders even though the file is valid. The MapLibre layer's source-layer does not equal tippecanoe's -l value (buildings here) — read the real name back from the metadata and fix the style, not the source name.
• Attributes missing in dense tiles. Tippecanoe trims attributes past the per-tile size budget in an order you did not choose; whitelist the fields you need with -y building_id -y area_m2 so they always survive.
• PMTiles returns 403 or won't load. The bucket must honour HTTP Range requests and send CORS headers; enable both on the object store, since the client fetches partial byte ranges cross-origin.
• The .pmtiles file is enormous. You are over-detailed at high zoom — cap --maximum-zoom to what users actually reach and add --drop-smallest-as-needed. Each extra zoom level roughly quadruples the tile count.
• The build dies with "no space left on device" after an hour. Tippecanoe's scratch files outgrew /tmp; point -t at a real disk with several times the input's size free.
• -P makes no difference to build time. It only parallelises line-delimited input. Convert with -f GeoJSONSeq (or FlatGeobuf) rather than plain GeoJSON.
• A numeric attribute arrives in the browser as a string. The source column was object dtype, or the intermediate text file lost the type. Force it at tile time with -T area_m2:float, and fix the dtype in Python so the next rebuild does not need the flag.
• Tiles render as garbage or the shim throws on decompress. The host is applying Content-Encoding: gzip to already-gzipped tiles. Serve the .pmtiles object with no transfer encoding at all.
• A second layer never appears. One tippecanoe run produces one layer; build each layer separately and merge with tile-join before converting to PMTiles.

Frequently Asked Questions

Why convert to GeoJSON at all — can't tippecanoe read GeoParquet? No. Tippecanoe accepts GeoJSON, line-delimited GeoJSON, CSV, and FlatGeobuf, but not GeoParquet or GeoPackage. ogr2ogr is the bridge. For very large inputs, prefer FlatGeobuf (-f FlatGeobuf) over GeoJSONSeq — it streams into tippecanoe without materializing a multi-gigabyte text file.

Should I keep the intermediate MBTiles file? Only if you also serve tiles from a process. PMTiles and MBTiles hold the same pyramid; MBTiles is a SQLite database you query per request (see Serving MBTiles with Python), while PMTiles is the static-hosting form. For a public, read-only map, delete the MBTiles after conversion and keep only the .pmtiles.

How do I bust the CDN cache after rebuilding tiles? Set a long, immutable cache header (Cache-Control: public, max-age=31536000, immutable) and change the filename on each rebuild — buildings.v3.pmtiles — rather than overwriting in place. That way a client mid-session never mixes tiles from two builds.

How do I update one region without rebuilding the whole archive? Tile the changed region on its own, then merge with tile-join, which replaces overlapping tiles from the later input. That is minutes instead of hours for a continental dataset where one metropolitan area changed. The alternative — regenerating from the GeoParquet store every time — stays the right default while the full build is under about ten minutes, because a single deterministic build has no merge order to reason about.

Should the archive be one layer or many? One tippecanoe invocation produces exactly one layer, so the question is really how many builds you run before tile-join. Split when the layers want different treatment: a street network needs deeper zooms than an administrative boundary, and a point layer wants dropping where a polygon layer wants coalescing. Keep them together when they share a zoom range and shedding strategy, since one build is simpler to reason about and produces slightly smaller output through shared tile deduplication.

How much does the whole pipeline cost to run at scale? The build is a one-off CPU and scratch-disk expense — an hour of an eight-core machine and tens of gigabytes of temporary space for a ten-million-feature dataset. Serving costs are storage for one object plus egress for the byte ranges clients actually read, which for a typical session is a few hundred kilobytes across a handful of tiles. There is no per-request compute and no database, which is the entire economic argument for this container over a running tile service.

Do I need Web Mercator (EPSG:3857) anywhere in this pipeline? Never as an input and never for measurement. Tippecanoe projects to the Web Mercator tile grid internally; your job is to hand it EPSG:4326 and to compute any areas or distances in an appropriate projected CRS (a UTM zone, a national grid) beforehand. Measuring in 3857 or 4326 gives distorted numbers.