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 — no tile server. It is for anyone whose GeoJSON has grown too large to inline and who wants cloud-native, server-less tiles. It sits under Vector Tile Pipelines with PMTiles in Web Mapping & Interactive Visualization.

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 operational cost. PMTiles packs the whole tile pyramid into one file with an internal index, and clients fetch only the byte ranges for the tiles in view. The build, though, has a sharp edge: tippecanoe ingests GeoJSON or FlatGeobuf, not GeoParquet directly, so you convert formats with ogr2ogr first. Skipping that step, or feeding tippecanoe a projected CRS, is where most pipelines fail.

GeoParquet is the right source format because it preserves CRS metadata, compresses well, and is the analytical store used in Cloud-Native Geospatial Formats — tiles are just a derived rendering view of it.

Prerequisites

conda install -c conda-forge "gdal=3.8.*" "pmtiles=3.2.*" "geopandas=0.14.*"
# tippecanoe: brew install tippecanoe (macOS) or build felt/tippecanoe (Linux)

Step-by-Step Implementation

1. Prepare the GeoParquet source in EPSG:4326. Do metric work before this step.

import geopandas as gpd

# building_footprints analysed in EPSG:25832 (metric) upstream
building_footprints = gpd.read_file("building_footprints.gpkg")
building_footprints["area_m2"] = building_footprints.geometry.area
building_footprints = building_footprints.to_crs(epsg=4326)
building_footprints[["building_id", "area_m2", "geometry"]].to_parquet(
    "buildings_4326.parquet"
)

2. Convert GeoParquet to line-delimited GeoJSON with ogr2ogr. Tippecanoe streams this efficiently.

import subprocess

subprocess.run([
    "ogr2ogr", "-f", "GeoJSONSeq", "buildings.geojsonl", "buildings_4326.parquet",
], check=True)

3. Tile with tippecanoe. Let it pick max zoom and shed density where tiles overflow.

subprocess.run([
    "tippecanoe",
    "-o", "buildings.mbtiles",
    "-zg",                          # auto max zoom
    "--coalesce-densest-as-needed", # merge features in crowded tiles
    "--extend-zooms-if-still-dropping",
    "-l", "buildings",              # layer name MapLibre will reference
    "buildings.geojsonl",
], check=True)

4. Convert MBTiles to PMTiles.

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

5. Reference it from MapLibre. Register the protocol, then add a vector source.

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

Verification

Inspect the PMTiles header to confirm zoom range, bounds, and the layer name before deploying.

from pmtiles.reader import Reader, MmapSource

with open("buildings.pmtiles", "rb") as fh:
    reader = Reader(MmapSource(fh))
    header = reader.header()
    print("min/max zoom:", header["min_zoom"], header["max_zoom"])  # 0 14
    print("bounds (E7):", header["min_lon_e7"], header["min_lat_e7"])
    meta = reader.metadata()
    print("layers:", [layer["id"] for layer in meta["vector_layers"]])  # ['buildings']

assert header["max_zoom"] >= 12, "Max zoom too low for street-level detail"

The reported bounds, converted from E7 integers, should match your dataset's extent — a wildly wrong box means the source wasn't in EPSG:4326.

Edge Cases & Debugging

Related