How to publish a 3D building layer to a web map

Problem statement

There are three ways to get buildings into a browser in 3D, and the difference between them is two orders of magnitude of effort. Extruding footprints client-side is a fill-extrusion layer and a GeoJSON. A single glTF is a file and a model layer. 3D Tiles is a pipeline.

Picking the heaviest option for a neighbourhood is the common mistake, and it is expensive: you build a tiling pipeline for a few thousand buildings that a browser would have extruded from a GeoJSON without noticing.

This guide covers all three, what each costs, and the thresholds where you move up. The numbers come from a real district: the 25,225 buildings in a 3 × 3 km square of The Hague, with roof and ground heights from the 3D BAG, the Dutch national 3D building model.

Quick answer

For anything up to a few tens of thousands of buildings, ship footprints with a height property and let the GPU extrude:

import os
import geopandas as gpd

web = b[["id", "height_m", "year", "geometry"]].to_crs(4326)
web["base_m"] = 0.0
web.to_file("buildings.geojson", driver="GeoJSON", COORDINATE_PRECISION=6)
print(f"{len(web):,} buildings, {os.path.getsize('buildings.geojson') / 1e6:.1f} MB, "
      f"height a {web['height_m'].dtype}")

25,225 buildings, 10.4 MB, height a float64

map.addLayer({
  id: 'buildings-3d', type: 'fill-extrusion', source: 'buildings',
  paint: {
    'fill-extrusion-height': ['get', 'height_m'],
    'fill-extrusion-base': ['get', 'base_m'],
    'fill-extrusion-color': [
      'interpolate', ['linear'], ['get', 'height_m'],
      0, '#e8f4fd', 10, '#7dd3fc', 25, '#0ea5e9', 50, '#1a3a6b'
    ],
    'fill-extrusion-opacity': 0.9
  }
});

The payload stays 2D, the browser does the extrusion, and the data is still queryable — a click returns a feature with its properties, which a mesh does not.

Route by scale, prepare vectors, set base and height

1. Choose the route by scale

• Up to a few tens of thousands of buildings: GeoJSON or vector tiles with fill-extrusion.
• One block, with real roofs: a single glTF placed as a model layer.
• A whole city, with LoD and streaming: 3D Tiles.

2. Prepare the vector data properly

Reproject to EPSG:4326, drop every column the map does not use, round coordinates, and simplify the footprints to something appropriate for the zoom levels you serve. The 3D BAG publishes 66 attributes per building; the map needs three:

sizes = {}
b.to_crs(4326).to_file("all_columns.geojson", driver="GeoJSON", COORDINATE_PRECISION=6)
sizes["all 3D BAG columns"] = os.path.getsize("all_columns.geojson")
web.to_file("full_precision.geojson", driver="GeoJSON")
sizes["3 columns, full precision"] = os.path.getsize("full_precision.geojson")
sizes["3 columns, 6 decimals"] = os.path.getsize("buildings.geojson")

simplified = web.assign(geometry=web.geometry.simplify(0.000015, preserve_topology=True))
simplified.to_file("buildings_simplified.geojson", driver="GeoJSON", COORDINATE_PRECISION=6)
sizes["+ simplified to ~1 m"] = os.path.getsize("buildings_simplified.geojson")

import shapely
for name, size in sizes.items():
    print(f"{name:28} {size / 1e6:6.1f} MB  {size / len(web):5.0f} bytes per building")
print("vertices:", shapely.get_num_coordinates(web.geometry.values).sum(), "->",
      shapely.get_num_coordinates(simplified.geometry.values).sum())

all 3D BAG columns             61.8 MB   2449 bytes per building
3 columns, full precision      14.6 MB    578 bytes per building
3 columns, 6 decimals          10.4 MB    413 bytes per building
+ simplified to ~1 m            8.4 MB    335 bytes per building
vertices: 231504 -> 151735

Six decimal places of a degree is about 10 cm — finer than any footprint is surveyed. A tolerance of 0.000015° is about 1.7 m north–south and 1 m east–west at The Hague's latitude.

3. Move to vector tiles when the GeoJSON gets big

Above a few megabytes, a single GeoJSON stalls the browser on load, because it has to be downloaded and parsed before anything draws. PMTiles gives you tiled delivery from static hosting with no server (Example 2).

4. Set base as well as height

fill-extrusion-base lets a building start above zero, which is how you stack parts of a building and how you place buildings on their ground height when there is no terrain. The 3D BAG gives both heights relative to NAP, the Dutch datum: ground (b3_h_maaiveld) and a representative roof height (b3_h_70p), so height_m here is their difference.

5. Decide about terrain

MapLibre and Mapbox can drape extrusions over a terrain source. With terrain enabled, the renderer places extrusions on the terrain surface, so fill-extrusion-base becomes an offset above the ground and a base of 0 is the right answer — the opposite of flat-map code that put each building's ground height into the base.

6. Style by data, not by height alone

Colouring by use, age or energy rating is what makes a 3D map informative rather than decorative. Height is already encoded by the geometry. The 3D BAG carries each building's construction year from the national address register:

decades = (b["year"] // 10 * 10).value_counts().sort_index()
print(decades[decades.index >= 1850].to_string())

year
1850     259
1860     234
1870     273
1880    1191
1890    2653
1900    5565
1910    5098
1920    3591
1930     932
1940     160
1950     661
1960     226
1970     167
1980    1002
1990    1378
2000     783
2010     601
2020     113

7. Budget the payload

A browser will hold a few tens of megabytes of geometry comfortably. Count your features, measure the bytes per building on your own data, and test on a mid-range phone rather than a workstation. The whole budget for this district, every route measured in step 2 and Example 2, is in the chart below.

3D BAG buildings as GeoJSON, PMTiles and glTF

The data used on this page

Every 3D BAG building whose footprint intersects a 3 × 3 km square of central The Hague (78,000–81,000 E, 454,000–457,000 N in the Dutch RD grid), from the 3D BAG's WFS (CC BY 4.0), paged 10,000 at a time. b keeps all 66 published attributes and adds id, height_m (the 70th percentile of the roof points minus the ground height), ground_z (ground height above NAP) and year (the construction year).

import io
import os

import geopandas as gpd
import pandas as pd
import requests

if not os.path.exists("buildings_rd.gpkg"):
    pages = []
    while True:
        reply = requests.get("https://data.3dbag.nl/api/BAG3D/wfs", timeout=300, params={
            "service": "WFS", "version": "2.0.0", "request": "GetFeature", "typeNames": "BAG3D:lod12",
            "outputFormat": "application/json", "bbox": "78000,454000,81000,457000,EPSG:28992",
            "count": 10_000, "startIndex": 10_000 * len(pages)})
        reply.raise_for_status()
        page = gpd.read_file(io.BytesIO(reply.content))
        if page.empty:
            break
        pages.append(page)
    pd.concat(pages).set_crs(28992, allow_override=True).to_file("buildings_rd.gpkg")

b = gpd.read_file("buildings_rd.gpkg")
b["id"] = b["identificatie"]
b["height_m"] = (b["b3_h_70p"] - b["b3_h_maaiveld"]).round(1)
b["ground_z"] = b["b3_h_maaiveld"]
b["year"] = b["oorspronkelijkbouwjaar"]
print(f"{len(b):,} buildings, {b['id'].duplicated().sum()} ids repeated (buildings the 3D BAG splits into parts); "
      f"height {b['height_m'].min()}–{b['height_m'].max()} m, median {b['height_m'].median()} m")

25,225 buildings, 75 ids repeated (buildings the 3D BAG splits into parts); height 1.0–71.8 m, median 9.9 m

Example 1 — the vector route, end to end

import geopandas as gpd, pathlib

def prepare_web_buildings(b, out="web_buildings.geojson", simplify_deg=0.000015,
                          keep=("id", "height_m", "year")):
    web = b.to_crs(4326).copy()
    web["geometry"] = web.geometry.simplify(simplify_deg, preserve_topology=True)
    web["height_m"] = pd.to_numeric(web["height_m"], errors="coerce")     # a string height extrudes to 0
    web = web[[*keep, "geometry"]]
    web = web[web.geometry.notna() & ~web.geometry.is_empty & web["height_m"].notna()]
    web.to_file(out, driver="GeoJSON", COORDINATE_PRECISION=6)

    size_kb = pathlib.Path(out).stat().st_size / 1024
    print(f"{len(web):,} buildings, {size_kb:,.0f} KB "
          f"({size_kb * 1024 / len(web):.0f} bytes each)")
    return out

prepare_web_buildings(b)

25,225 buildings, 7,871 KB (320 bytes each)

Bytes per building is the number to watch. At about 320 bytes each, this district is an 8 MB file — past the point where one GeoJSON is comfortable, which is why the next example tiles it.

Example 2 — the same data as PMTiles

GDAL (3.8 or later) writes PMTiles directly, so no separate tiling tool is needed:

import subprocess

for minzoom in (10, 13):
    out = f"buildings_z{minzoom}.pmtiles"
    pathlib.Path(out).unlink(missing_ok=True)
    subprocess.run(["ogr2ogr", "-f", "PMTiles", out, "web_buildings.geojson", "-nln", "buildings",
                    "-dsco", f"MINZOOM={minzoom}", "-dsco", "MAXZOOM=16"], check=True)
    sizes[f"PMTiles, zooms {minzoom}-16"] = os.path.getsize(out)
    print(f"{out}: {os.path.getsize(out) / 1e6:.1f} MB")

buildings_z10.pmtiles: 4.1 MB
buildings_z13.pmtiles: 2.7 MB

A minimum zoom of 13 is deliberate: at zoom 12 a pixel at this latitude covers about 23 m, so a house is sub-pixel and tiles below 13 carry geometry nobody can see. Extrusions are only meaningful when the map is pitched and zoomed in anyway. And unlike the GeoJSON, a viewer only ever downloads the tiles on screen.

For comparison, the mesh route — every footprint extruded into a triangle mesh and written as a binary glTF:

import numpy as np
import trimesh

origin = np.array([78_000.0, 454_000.0, 0.0])
meshes = [trimesh.creation.extrude_polygon(part, h).apply_translation(-origin)
          for geom, h in zip(b.geometry, b["height_m"]) for part in getattr(geom, "geoms", [geom])]
district = trimesh.util.concatenate(meshes)
district.export("district.glb")
sizes["one glTF mesh"] = os.path.getsize("district.glb")

footprint_vertices = int(shapely.get_num_coordinates(b.geometry.values).sum())
flat = district.copy()
flat.unmerge_vertices()                                  # one vertex per triangle corner, as flat shading needs
print(f"footprints: {footprint_vertices:,} vertices; mesh: {len(district.faces):,} triangles, "
      f"{len(flat.vertices):,} vertices flat-shaded ({len(flat.vertices) / footprint_vertices:.1f}x)")
print(f"district.glb: {sizes['one glTF mesh'] / 1e6:.1f} MB")

footprints: 231,504 vertices; mesh: 724,216 triangles, 2,172,648 vertices flat-shaded (9.4x)
district.glb: 13.6 MB

Example 3 — a single glTF for one block, placed on the map

import json
import pyproj
from pyproj import Transformer

pyproj.network.set_network_enabled(True)                      # fetch the Dutch geoid grid on demand

origin_rd = (78_642.0, 457_940.0, 0.0)                 # RD x, y and NAP height of the model origin
to_wgs84 = Transformer.from_crs("EPSG:7415", "EPSG:4979", always_xy=True)   # RD + NAP -> WGS 84 3D
lon, lat, ellipsoidal = to_wgs84.transform(*origin_rd)

placement = {
    "model": "block.glb",
    "origin_crs": "EPSG:7415",
    "origin": list(origin_rd),
    "origin_wgs84": [round(lon, 7), round(lat, 7), round(ellipsoidal, 2)],
    "axis_convention": "Y-up, rotated -90° about X from Z-up",
    "vertical_datum": "NAP; add the geoid separation for ellipsoidal placement",
}
pathlib.Path("block.placement.json").write_text(json.dumps(placement, indent=2))
print(placement["origin_wgs84"])

[4.2726524, 52.1044255, 43.4]

The vertical datum note is the one that saves an afternoon. NAP 0 at this point is 43.4 m above the WGS 84 ellipsoid, so a model in NAP heights placed with those numbers as ellipsoidal heights sinks by 43 m. The transformation needs PROJ's Dutch geoid grid, which PROJ fetches from its CDN once network access is on, as above (or with PROJ_NETWORK=ON). Without the grid PROJ silently returns the NAP height unchanged — 0.0 here — which is the sinking model again.

Example 4 — check what the renderer does with a bad height

The renderer's own expression engine answers the "why is nothing extruded" question without opening a browser. MapLibre publishes it as @maplibre/maplibre-gl-style-spec (npm install @maplibre/maplibre-gl-style-spec):

%%writefile check_height.mjs
import * as spec from '@maplibre/maplibre-gl-style-spec';

const property = spec.latest['paint_fill-extrusion']['fill-extrusion-height'];
const height = spec.createPropertyExpression(['get', 'height_m'], property).value;
const warnings = [];
console.warn = (message) => warnings.push(message);
for (const properties of [{height_m: 12.3}, {height_m: '12.3'}, {height_m: null}, {}]) {
  const value = height.evaluate({zoom: 15}, {type: 'Polygon', properties});
  console.log(JSON.stringify(properties).padEnd(20), '->', value);
}
console.log('warnings:', warnings);
const cast = spec.createPropertyExpression(['to-number', ['get', 'height_m'], 0], property).value;
console.log("['to-number', ...] on '12.3' ->", cast.evaluate({zoom: 15}, {type: 'Polygon', properties: {height_m: '12.3'}}));

print(subprocess.run(["node", "check_height.mjs"], capture_output=True, text=True).stdout)

{"height_m":12.3}    -> 12.3
{"height_m":"12.3"}  -> 0
{"height_m":null}    -> 0
{}                   -> 0
warnings: [
  'Expected value to be of type number, but found string instead.',
  'Expected value to be of type number, but found null instead.'
]
['to-number', ...] on '12.3' -> 12.3

A string height and a null both extrude to 0, with a line in the browser console that nobody reads; a missing height extrudes to 0 with no warning at all. Cast before export, or wrap the expression in to-number.

Client-side extrusion and zoom limits

Why client-side extrusion is usually right

The GPU extrudes a polygon into a prism essentially for free; the work is uploading the polygon. Sending the prism instead means sending about ten times the vertices for the same information once each wall face gets its own normals, as Example 2 measured, plus losing the ability to restyle, filter or query without regenerating the data. Binary glTF claws some of that back — the district's mesh was not ten times the size of its GeoJSON — but it is still a single file that must arrive whole. The exception is when the geometry is not a prism — real roofs — and then you are sending a mesh because you need one.

Why zoom limits matter more than simplification

An extruded building is only visible on a pitched, zoomed-in map. Serving building tiles at zoom 10 costs bandwidth for geometry that occupies a fraction of a pixel and cannot be seen at any pitch. Setting a minimum zoom is a bigger saving than any amount of coordinate rounding.

Why fill-extrusion-base behaves differently with terrain

Without a terrain source, the base is metres above the map's zero plane. With terrain enabled, the renderer places extrusions on the terrain surface, so a base of 0 sits on the ground and a non-zero base floats above it. Code written for one case is wrong in the other, and the symptom is buildings buried in a hill or hovering over it.

Why 3D Tiles is a different kind of commitment

3D Tiles is a streaming format with a tiling pipeline, a level-of-detail tree, a georeferencing transform and a client that must support it. It is the right answer for a city-scale model with real roofs and per-feature metadata, and it is a lot of machinery for a neighbourhood. The honest test is whether the payload fits in a browser; if it does, do not tile it.

Opacity, shadows and numeric heights

• Opacity is expensive. Semi-transparent extrusions have to be drawn in the right order to look right.
• Extrusions do not cast shadows in MapLibre; the light model is ambient plus directional.
• height must be a number. A string property extrudes to zero, with only a console warning.
• Missing heights render as zero — with no warning at all. Filter them out or substitute a default.
• Pitch limits vary by renderer. Test the pitch your users will actually use.
• Mobile GPUs choke sooner. Test on a phone.
• Keep the id. A click has to resolve back to the data.
• Colour by data, not by height. The geometry already shows height.

Internal links

• 3D Tiles and the web delivery of 3D city models explained — the heavyweight route
• How to export a 3D city model to glTF in Python — producing the middle route
• A 3D viewer shows nothing or a black screen — when it loads and shows nothing
• How to extrude building footprints into 3D in Python — why you often should not
• How to prepare a GeoDataFrame for the web — trimming the payload
• How to serve PMTiles with MapLibre — tiled delivery without a server
• Vector tiles explained — what tippecanoe is producing
• A web map layer is not appearing — the 2D version of the same debugging

FAQ

What is the simplest way to show 3D buildings on a web map?

A GeoJSON of footprints with a numeric height property and a MapLibre fill-extrusion layer. The browser extrudes on the GPU and the payload stays 2D.

When should I move to vector tiles?

When the GeoJSON exceeds a few megabytes — for the 25,225 buildings of a 3 × 3 km square of The Hague it was 8–10 MB even with three columns. PMTiles gives tiled delivery from static hosting with no server, and GDAL writes it directly.

When do I need 3D Tiles?

At city scale with real roof geometry and per-feature metadata. For a neighbourhood it is a great deal of machinery for no benefit.

Why are my buildings buried in the terrain?

fill-extrusion-base is relative to the terrain when a terrain source is enabled and to the map plane when it is not. A base copied from the flat case sinks buildings into a hill.

Why is nothing extruded?

Usually because the height property is a string rather than a number, or is missing. MapLibre's expression engine evaluates both to zero: the string with a warning in the browser console, the missing value with none.

What zoom levels should I serve buildings at?

Thirteen and above. At zoom 12 a pixel at The Hague's latitude is about 23 m, so a house is sub-pixel and cannot be seen at any pitch; the zoom-13 PMTiles archive was also a third smaller than one starting at zoom 10.