Source code: PMTiles vector tiles in Python and a MapLibre GL JS web map
Runnable Python scripts, requirements and console output for the tutorial “PMTiles vector tiles in Python and a MapLibre GL JS web map”.
Last updated
These are the exact scripts used to produce every figure and number in PMTiles vector tiles in Python and a MapLibre GL JS web map. Download them, install the requirements into a fresh virtual environment and run them from this folder.
- pmtiles_vector_tiles.py
- serve.py
- requirements.txt
- console.txt: the console output recorded when we ran it
python -m venv .venv.venv\Scripts\activate # macOS/Linux: source .venv/bin/activatepip install -r requirements.txtpython pmtiles_vector_tiles.pypmtiles_vector_tiles.py
"""Turn Overture GeoParquet into PMTiles vector tiles with GDAL, and build a static MapLibre GL JS map.
TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/vector-tiles-pmtiles-maplibre/Part 2 of "Cloud-native vector GIS with Python". It reads the Bristol GeoParquet files written by Part 1(geoparquet-duckdb-overture-maps). If they are missing, it runs the Part 1 script first.
Tiles are written by GDAL's PMTiles driver (GDAL >= 3.8), which ships inside the pyogrio wheel, so thisruns on Windows without Tippecanoe. Then: python serve.py # open http://localhost:8000/ (a server with HTTP Range support) node screenshot_map.mjs # optional: Playwright screenshots of the running map
Data: Overture Maps Foundation. Buildings and water: ODbL-1.0 (© OpenStreetMap contributors, Microsoft).Places: CDLA-Permissive-2.0 / Apache-2.0 / CC0-1.0, depending on the source of each record."""import jsonimport mathimport subprocessimport sysimport timefrom pathlib import Path
import geopandas as gpdimport matplotlibimport pandas as pdimport pyogrioimport requestsimport shapelyfrom pmtiles.reader import MmapSource, Reader, all_tiles
matplotlib.use("Agg")import matplotlib.pyplot as plt # noqa: E402from matplotlib.colors import LogNorm # noqa: E402from matplotlib.patches import Rectangle # noqa: E402
HERE = Path(__file__).resolve().parentPART1 = HERE.parent / "geoparquet-duckdb-overture-maps"DATA = HERE / "data" # experiment tilesets (not committed)OUT = HERE / "output" # figures, console log, and the finished web map in output/web/WEB = OUT / "web"for d in (DATA, OUT, WEB, WEB / "lib"): d.mkdir(parents=True, exist_ok=True)
RELEASE = "2026-09-23.1" # the Overture release Part 1 downloadedMAPLIBRE = "5.24.0" # maplibre-gl, self-hosted from output/web/lib/PMTILES_JS = "4.5.0" # pmtiles (JavaScript), inlined into the HTML pageATTRIBUTION = ('© <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors, ' '<a href="https://docs.overturemaps.org/attribution/">Overture Maps Foundation</a> ' f'(release {RELEASE})')# The copy of the map embedded in the tutorial loads its tiles and libraries from here on technicalgis.comSITE_BASE = "/demos/vector-tiles-pmtiles-maplibre/"
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
def section(title: str) -> None: print(f"\n=== {title} ===")
def write_pmtiles(gdf: gpd.GeoDataFrame, path: Path, layer: str, **options) -> float: """Write one layer to a new PMTiles file with GDAL; return the time taken in seconds.""" path.unlink(missing_ok=True) t0 = time.perf_counter() # encoding="UTF-8": otherwise pyogrio may encode names with the Windows code page and fail on characters # such as Welsh or Polish letters pyogrio.write_dataframe(gdf, path, driver="PMTiles", layer=layer, encoding="UTF-8", **{k: str(v) for k, v in options.items()}) return time.perf_counter() - t0
def tile_stats(path: Path) -> pd.DataFrame: """One row per tile: z, x, y and compressed size in bytes, read from the PMTiles directory.""" with open(path, "rb") as f: reader = Reader(MmapSource(f)) rows = [(z, x, y, len(data)) for (z, x, y), data in all_tiles(reader.get_bytes)] return pd.DataFrame(rows, columns=["z", "x", "y", "bytes"])
def tile_bounds(z: int, x: int, y: int) -> tuple[float, float, float, float]: """Lon/lat bounds of an XYZ web-mercator tile.""" n = 2 ** z lat = lambda t: math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * t / n)))) # noqa: E731 return x / n * 360 - 180, lat(y + 1), (x + 1) / n * 360 - 180, lat(y)
section("versions")print(f"python {sys.version.split()[0]} | geopandas {gpd.__version__} | pyogrio {pyogrio.__version__} " f"(GDAL {pyogrio.__gdal_version_string__}) | shapely {shapely.__version__}")drivers = pyogrio.list_drivers()print(f"GDAL vector-tile drivers: PMTiles={drivers.get('PMTiles')} MVT={drivers.get('MVT')}")
# ---------------------------------------------------------------------------- 1. load Part 1's GeoParquetsection("load Part 1 GeoParquet")names = ["buildings", "places", "water", "hex"]paths = {n: PART1 / "output" / f"bristol_{n}.parquet" for n in names}if not all(p.exists() for p in paths.values()): print("Part 1 output not found: running geoparquet-duckdb-overture-maps/overture_duckdb.py first") subprocess.run([sys.executable, str(PART1 / "overture_duckdb.py")], check=True, cwd=PART1)layers = {n: gpd.read_parquet(p) for n, p in paths.items()}for n, g in layers.items(): print(f"{n:<9} {len(g):>7,} features {paths[n].stat().st_size / 1e6:5.2f} MB GeoParquet crs {g.crs.to_string()}")
# Keep only what the web map needs: fewer attributes mean smaller tiles (measured in section 3)b = layers["buildings"]buildings = gpd.GeoDataFrame({"subtype": b["subtype"], "height": b["height"].round(1)}, geometry=b.geometry)p = layers["places"]places = gpd.GeoDataFrame({"name": p["name"], "category": p["basic_category"], "group": p["category_group"], "confidence": p["confidence"].round(2), "source": p["source"], "licence": p["licence"]}, geometry=p.geometry)water = layers["water"][["subtype", "geometry"]]hexes = layers["hex"][["h3", "buildings", "places", "food_drink", "built_pct", "geometry"]]
# ---------------------------------------------------------------------------- 2. zoom-level trade-offsection("buildings: maximum zoom vs size (MINZOOM=12, GDAL defaults otherwise)")rows, per_zoom = [], {}for maxzoom in [12, 13, 14, 15, 16]: path = DATA / f"buildings_z12-{maxzoom}.pmtiles" secs = write_pmtiles(buildings, path, "buildings", MINZOOM=12, MAXZOOM=maxzoom) stats = tile_stats(path) per_zoom[maxzoom] = stats rows.append({"zooms": f"12-{maxzoom}", "tiles": len(stats), "size_MB": round(path.stat().st_size / 1e6, 2), "write_s": round(secs, 1), "largest_tile_KB": round(stats["bytes"].max() / 1e3)})zoom_df = pd.DataFrame(rows)print(zoom_df.to_string(index=False))
by_z = per_zoom[16].groupby("z")["bytes"].agg( tiles="size", total_KB=lambda s: round(s.sum() / 1e3), mean_KB=lambda s: round(s.mean() / 1e3, 1), max_KB=lambda s: round(s.max() / 1e3))tile_m = {z: 40_075_016.7 * math.cos(math.radians(51.46)) / 2 ** z for z in by_z.index}by_z["tile_width_m"] = [round(tile_m[z]) for z in by_z.index]by_z["grid_step_m"] = [round(tile_m[z] / 4096, 2) for z in by_z.index] # EXTENT=4096 grid units per tileprint("per zoom level in the z12-16 archive (compressed tile sizes):")print(by_z.to_string())
# What each zoom level holds: GDAL keeps tiles under MAX_SIZE (500 KB) by coarsening or dropping featuresfor z in [12, 13, 14, 16]: g = pyogrio.read_dataframe(DATA / "buildings_z12-16.pmtiles", layer="buildings", ZOOM_LEVEL=str(z)) print(f"z{z}: {len(g):,} feature pieces read back (a building cut by a tile edge counts once per tile)")
# ---------------------------------------------------------------------------- 3. attributes and precisionsection("buildings at z14 only: attributes, EXTENT and SIMPLIFICATION")variants = { "geometry only": (buildings[["geometry"]], {}), "subtype + height (used)": (buildings, {}), "all 10 Part 1 columns": (b, {}), "subtype + height, EXTENT=1024": (buildings, {"EXTENT": 1024}), "subtype + height, SIMPLIFICATION=4": (buildings, {"SIMPLIFICATION": 4}),}rows = []for i, (label, (gdf, opts)) in enumerate(variants.items()): path = DATA / f"buildings_z14_variant{i}.pmtiles" secs = write_pmtiles(gdf, path, "buildings", MINZOOM=14, MAXZOOM=14, **opts) rows.append({"variant": label, "size_MB": round(path.stat().st_size / 1e6, 2), "write_s": round(secs, 1)})attr_df = pd.DataFrame(rows)print(attr_df.to_string(index=False))
# ---------------------------------------------------------------------------- 4. the tilesets for the web mapsection("final tilesets for the web map")final = { # layer: (data, minzoom, maxzoom). MapLibre overzooms the last level, so z14 tiles still draw at z19 "hex": (hexes, 10, 13), "water": (water, 10, 14), "buildings": (buildings, 14, 14), "places": (places, 14, 14),}rows = []for name, (gdf, z0, z1) in final.items(): path = WEB / f"{name}.pmtiles" secs = write_pmtiles(gdf, path, name, MINZOOM=z0, MAXZOOM=z1, NAME=f"Bristol {name}", DESCRIPTION=f"Overture Maps {RELEASE}, central Bristol") rows.append({"file": path.name, "features": len(gdf), "zooms": f"{z0}-{z1}", "tiles": len(tile_stats(path)), "size_MB": round(path.stat().st_size / 1e6, 2), "write_s": round(secs, 1)})final_df = pd.DataFrame(rows)print(final_df.to_string(index=False))print(f"total: {final_df['size_MB'].sum():.2f} MB in {final_df['tiles'].sum()} tiles " f"(the four GeoParquet inputs are {sum(p.stat().st_size for p in paths.values()) / 1e6:.1f} MB)")meta = pyogrio.read_info(WEB / "buildings.pmtiles", layer="buildings")["dataset_metadata"]print("buildings.pmtiles metadata:", {k: meta[k] for k in ["name", "minzoom", "maxzoom", "bounds", "format"]})
# ---------------------------------------------------------------------------- 5. the MapLibre pagesection("MapLibre GL JS page")for url, dest in [(f"https://unpkg.com/maplibre-gl@{MAPLIBRE}/dist/maplibre-gl.js", "maplibre-gl.js"), (f"https://unpkg.com/maplibre-gl@{MAPLIBRE}/dist/maplibre-gl.css", "maplibre-gl.css"), (f"https://unpkg.com/maplibre-gl@{MAPLIBRE}/LICENSE.txt", "LICENSE-maplibre-gl.txt"), (f"https://unpkg.com/pmtiles@{PMTILES_JS}/dist/pmtiles.js", "pmtiles.js")]: target = WEB / "lib" / dest if not target.exists(): response = requests.get(url, timeout=60) response.raise_for_status() target.write_bytes(response.content)pmtiles_js = (WEB / "lib" / "pmtiles.js").read_text(encoding="utf-8")
xmin, ymin, xmax, ymax = layers["hex"].total_boundsGROUP_COLOURS = {"food_and_drink": "#c21f3a", "shopping": "#e8a33c", "services_and_business": "#1543b0", "health_care": "#1b9e77", "arts_and_entertainment": "#7b3294", "education": "#5e3c99"}place_colour = ["match", ["get", "group"], *[v for kv in GROUP_COLOURS.items() for v in kv], "#6b6b6b"]
def style(base: str) -> dict: """MapLibre style JSON: one vector source per PMTiles file, no third-party basemap.""" sources = {n: {"type": "vector", "url": f"pmtiles://{base}{n}.pmtiles", "minzoom": z0, "maxzoom": z1, "attribution": ATTRIBUTION} for n, (_, z0, z1) in final.items()} return { "version": 8, "sources": sources, "layers": [ {"id": "background", "type": "background", "paint": {"background-color": "#f3f1ec"}}, {"id": "water", "type": "fill", "source": "water", "source-layer": "water", "paint": {"fill-color": "#9ecae1"}}, {"id": "hex", "type": "fill", "source": "hex", "source-layer": "hex", "maxzoom": 14.5, "paint": {"fill-color": ["interpolate", ["linear"], ["get", "places"], 0, "#fff7ec", 25, "#fdd49e", 75, "#fc8d59", 200, "#d7301f", 500, "#7f0000"], "fill-outline-color": "#ffffff", "fill-opacity": ["interpolate", ["linear"], ["zoom"], 13.5, 0.85, 14.5, 0]}}, {"id": "buildings", "type": "fill", "source": "buildings", "source-layer": "buildings", "minzoom": 14, "paint": {"fill-color": ["match", ["get", "subtype"], "residential", "#e8c48f", "commercial", "#e49aa6", "education", "#b7a5d6", "medical", "#9fd3c0", "#c9cdca"], "fill-outline-color": "#8a908d"}}, {"id": "places", "type": "circle", "source": "places", "source-layer": "places", "minzoom": 15, "paint": {"circle-color": place_colour, "circle-stroke-color": "#ffffff", "circle-stroke-width": 1, "circle-radius": ["interpolate", ["linear"], ["zoom"], 15, 3, 18, 6]}}, ], }
LEGEND = "".join(f'<span><i style="background:{c}"></i>{g.replace("_", " ")}</span>' for g, c in GROUP_COLOURS.items())PAGE = """<!doctype html><html lang="en-GB"><head><meta charset="utf-8"><meta name="viewport" content="width=device-width, initial-scale=1"><meta name="robots" content="noindex, follow"><title>Bristol buildings and places: PMTiles and MapLibre GL JS</title><link rel="stylesheet" href='__BASE__lib/maplibre-gl.css'><script src='__BASE__lib/maplibre-gl.js'></script><!-- pmtiles __PMTILES_VERSION__ (BSD-3-Clause, https://github.com/protomaps/PMTiles), inlined --><script>__PMTILES_JS__</script><style> html, body, #map { margin: 0; height: 100%; font: 13px/1.4 system-ui, sans-serif; } .panel { position: absolute; top: 10px; left: 10px; max-width: 260px; background: rgba(255,255,255,.93); padding: 8px 10px; border-radius: 6px; box-shadow: 0 1px 4px rgba(0,0,0,.25); } .panel b { display: block; margin-bottom: 2px; } .panel span { display: inline-flex; align-items: center; margin: 2px 8px 0 0; font-size: 12px; } .panel i { width: 10px; height: 10px; border-radius: 50%; margin-right: 4px; display: inline-block; } #zoom { color: #555; font-size: 12px; } @media (max-width: 480px) { .panel { max-width: 180px; font-size: 12px; } .panel span { display: none; } }</style></head><body><div id="map"></div><div class="panel"><b>Central Bristol from Overture Maps</b> <div id="zoom"></div> <div>Zoom 10–14: places per H3 hexagon. Zoom 14+: buildings. Zoom 15+: places (click one).</div> __LEGEND__</div><script>const protocol = new pmtiles.Protocol();maplibregl.addProtocol("pmtiles", protocol.tile);const params = new URLSearchParams(location.search);const map = new maplibregl.Map({ container: "map", style: __STYLE__, center: [+(params.get("lon") || -2.5925), +(params.get("lat") || 51.4585)], zoom: +(params.get("z") || 12.4), minZoom: 10, maxZoom: 19, maxBounds: [[__XMIN__ - 0.05, __YMIN__ - 0.03], [__XMAX__ + 0.05, __YMAX__ + 0.03]], attributionControl: { compact: window.innerWidth < 640 },});map.addControl(new maplibregl.NavigationControl(), "top-right");map.addControl(new maplibregl.ScaleControl({ unit: "metric" }), "bottom-left");// MapLibre opens the compact attribution on load; start it collapsed so it doesn't cover the scale bar on phones.const collapseAttrib = () => document.querySelectorAll(".maplibregl-ctrl-attrib.maplibregl-compact").forEach((a) => { a.classList.remove("maplibregl-compact-show"); a.removeAttribute("open");});map.once("load", collapseAttrib); map.once("idle", collapseAttrib);const showZoom = () => document.getElementById("zoom").textContent = "zoom " + map.getZoom().toFixed(1);map.on("zoom", showZoom); showZoom();map.on("click", "places", (e) => { const p = e.features[0].properties; new maplibregl.Popup().setLngLat(e.lngLat).setHTML( `<b>${p.name || "(no name)"}</b><br>${(p.category || "").replaceAll("_", " ")}<br>` + `source: ${p.source} (${p.licence})<br>confidence: ${p.confidence}`).addTo(map);});map.on("mouseenter", "places", () => map.getCanvas().style.cursor = "pointer");map.on("mouseleave", "places", () => map.getCanvas().style.cursor = "");window.map = map;</script></body></html>"""
def page(base: str) -> str: html = PAGE for key, value in {"__BASE__": base, "__PMTILES_VERSION__": PMTILES_JS, "__LEGEND__": LEGEND, "__STYLE__": json.dumps(style(base), ensure_ascii=False), "__XMIN__": f"{xmin:.4f}", "__YMIN__": f"{ymin:.4f}", "__XMAX__": f"{xmax:.4f}", "__YMAX__": f"{ymax:.4f}", "__PMTILES_JS__": pmtiles_js}.items(): html = html.replace(key, value) return html
(WEB / "index.html").write_text(page(""), encoding="utf-8")(OUT / "map.html").write_text(page(SITE_BASE), encoding="utf-8") # the copy embedded on technicalgis.comweb_files = sorted(f for f in WEB.rglob("*") if f.is_file())print(f"output/web/: {len(web_files)} files, {sum(f.stat().st_size for f in web_files) / 1e6:.2f} MB")for f in web_files: print(f" {f.relative_to(WEB).as_posix():<28} {f.stat().st_size / 1e3:8.1f} KB")print("serve it with: python serve.py then open http://localhost:8000/")
# ---------------------------------------------------------------------------- 6. figuressection("figures")
# Figure 1: zoom ladderfig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 3.8))a1.bar(zoom_df["zooms"], zoom_df["size_MB"], color="#1543b0")for i, (s, t) in enumerate(zip(zoom_df["size_MB"], zoom_df["write_s"])): a1.text(i, s, f"{s:g} MB\n{t:g} s", ha="center", va="bottom", fontsize=8)a1.set_ylim(0, zoom_df["size_MB"].max() * 1.25)a1.set_xlabel("zoom levels in the archive")a1.set_title("Buildings archive: size and GDAL write time", loc="left")a2.bar(by_z.index.astype(str), by_z["tiles"], color="#e8a33c")a2.set_yscale("log")for i, (n, kb) in enumerate(zip(by_z["tiles"], by_z["max_KB"])): a2.text(i, n, f"{n} tiles\nmax {kb} KB", ha="center", va="bottom", fontsize=8)a2.set_ylim(1, by_z["tiles"].max() * 8)a2.set_xlabel("zoom level")a2.set_title("Tiles per zoom level (log scale) and the largest tile", loc="left")for a in (a1, a2): a.spines[["top", "right"]].set_visible(False)fig.tight_layout()fig.savefig(OUT / "01-zoom-size-tradeoff.png")plt.close(fig)
# Figure 2: tile sizes at z14, drawn over the buildingsz14 = tile_stats(WEB / "buildings.pmtiles")fig, ax = plt.subplots(figsize=(9, 7.4))water.to_crs(3857).plot(ax=ax, color="#9ecae1", linewidth=0)norm = LogNorm(vmin=max(z14["bytes"].min() / 1e3, 0.1), vmax=z14["bytes"].max() / 1e3)cmap = plt.get_cmap("YlOrRd")for t in z14.itertuples(): x0, y0, x1, y1 = gpd.GeoSeries([shapely.box(*tile_bounds(t.z, t.x, t.y))], crs=4326).to_crs(3857).total_bounds ax.add_patch(Rectangle((x0, y0), x1 - x0, y1 - y0, facecolor=cmap(norm(t.bytes / 1e3)), edgecolor="#555555", linewidth=0.5, alpha=0.6, zorder=1)) ax.text((x0 + x1) / 2, (y0 + y1) / 2, f"{t.bytes / 1e3:.0f}" if t.bytes >= 1e3 else f"{t.bytes / 1e3:.1f}", ha="center", va="center", fontsize=7, zorder=3, bbox={"boxstyle": "round,pad=0.15", "facecolor": "white", "edgecolor": "none", "alpha": 0.8})buildings.to_crs(3857).plot(ax=ax, color="#1d2526", linewidth=0, alpha=0.55, zorder=2)ax.set_aspect("equal")ax.autoscale_view()ax.set_axis_off()fig.colorbar(plt.cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax, shrink=0.6, label="compressed tile size (KB)")ax.set_title(f"{len(z14)} zoom-14 building tiles, {z14['bytes'].sum() / 1e6:.2f} MB in total " f"(median {z14['bytes'].median() / 1e3:.0f} KB)", loc="left")ax.text(1, -0.01, "© OpenStreetMap contributors, Microsoft, Overture Maps Foundation (ODbL)", transform=ax.transAxes, ha="right", va="top", fontsize=7, color="#4a5553")fig.savefig(OUT / "02-z14-tile-sizes.png")plt.close(fig)
# Figure 3: attributes and precisionfig, ax = plt.subplots(figsize=(8, 2.8))ax.barh(attr_df["variant"], attr_df["size_MB"], color=["#9aa3a8", "#1543b0", "#c21f3a", "#1b9e77", "#e8a33c"])for y, v in enumerate(attr_df["size_MB"]): ax.text(v, y, f" {v:g} MB", va="center", fontsize=8)ax.invert_yaxis()ax.set_xlim(0, attr_df["size_MB"].max() * 1.2)ax.set_title(f"{len(buildings):,} buildings at zoom 14: what attributes and options cost", loc="left")ax.spines[["top", "right"]].set_visible(False)fig.savefig(OUT / "03-attributes-precision.png")plt.close(fig)
print("wrote", ", ".join(sorted(p.name for p in OUT.glob("0*.png"))))serve.py
"""Serve output/web/ on http://localhost:8000/ with HTTP Range support, which PMTiles needs.
Python's built-in `python -m http.server` ignores Range headers and always returns the whole file, so thepmtiles JavaScript library cannot read a tile directory from it. This adds byte-range responses (206).Usage: python serve.py [port]"""import osimport reimport sysfrom functools import partialfrom http.server import SimpleHTTPRequestHandler, ThreadingHTTPServerfrom pathlib import Path
WEB = Path(__file__).resolve().parent / "output" / "web"
class RangeHandler(SimpleHTTPRequestHandler): extensions_map = {**SimpleHTTPRequestHandler.extensions_map, ".pmtiles": "application/octet-stream"}
def send_head(self): match = re.fullmatch(r"bytes=(\d*)-(\d*)", self.headers.get("Range", "")) path = self.translate_path(self.path) if not match or not os.path.isfile(path): return super().send_head() size = os.path.getsize(path) start, end = match.groups() if start == "": # "bytes=-500": the last 500 bytes start, end = max(size - int(end), 0), size - 1 else: start, end = int(start), min(int(end) if end else size - 1, size - 1) if start >= size: self.send_error(416, "Requested range not satisfiable") return None f = open(path, "rb") f.seek(start) self.send_response(206) self.send_header("Content-Type", self.guess_type(path)) self.send_header("Content-Range", f"bytes {start}-{end}/{size}") self.send_header("Content-Length", str(end - start + 1)) self.send_header("Accept-Ranges", "bytes") self.end_headers() self._remaining = end - start + 1 return f
def copyfile(self, source, outputfile): remaining = getattr(self, "_remaining", None) if remaining is None: return super().copyfile(source, outputfile) while remaining > 0 and (chunk := source.read(min(65536, remaining))): outputfile.write(chunk) remaining -= len(chunk) self._remaining = None
if __name__ == "__main__": port = int(sys.argv[1]) if len(sys.argv) > 1 else 8000 server = ThreadingHTTPServer(("127.0.0.1", port), partial(RangeHandler, directory=str(WEB))) print(f"serving {WEB} on http://localhost:{port}/ (Ctrl+C to stop)") server.serve_forever()