Skip to content
CityParquet is v0.1.0-draft — the specification is still changing.
Esc
↑↓navigate↵open⌘Jpreview
On this page

Working with 3DBAG at cloud scale

Query a national-scale 3D city model directly on object storage with DuckDB — bounding-box search over 21.5 million objects, 3D measurement, semantic roof analysis, and format conversion, without downloading anything.

The earlier tutorials each take one tool and one local file. This one puts them together on a real, national-scale dataset that lives on object storage: the 3DBAG reconstruction of the Netherlands, converted to CityParquet as one 16.4 GB file on Cloudflare R2. Nothing is downloaded up front — every query below reads only the columns and row groups it needs, over plain HTTP range requests.

That is the whole argument for storing a city model as CityParquet: because it is columnar, a reader touches a tiny fraction of the bytes, so a single laptop can sweep 21.5 million city objects in seconds and then carve out one neighbourhood to measure in 3D.

This walkthrough shows six things end to end — reading the whole country at once, bounding-box search, GeoParquet interoperability, 3D solid measurement, semantic roof analysis, and conversion back to CityJSONSeq.

The dataset: 3DBAG as one CityParquet package

The package is a directory holding one object table and its STAC Item:

https://cityparquet.open3d.city/data/3dbag/
  building.parquet      -- 21 555 522 rows: every Building and BuildingPart in NL, 16.4 GB
  metadata.json         -- STAC Item (city3d extension)
  collection.json       -- the STAC collection holding that one Item

building.parquet carries the full CityParquet shape: typed attribute columns (the 3DBAG b3_* set), a bbox STRUCT of every object’s 3D extent, and — one per Level of Detail present (0.0, 1.2, 1.3, 2.2) — a geometry_lod* column with its paired geometry_properties_lod*, material_lod*, and texture_lod*. As in raw 3DBAG, geometry lives on BuildingPart rows while the descriptive attributes live on the parent Building; the two are linked by parents / children.

The metadata.json is a STAC Item using the city3d extension, so the package is discoverable without opening the Parquet at all:

{
  "type": "Feature",
  "stac_version": "1.1.0",
  "id": "3dbag",
  "properties": {
    "city3d:version": "2.0",
    "city3d:city_objects": 21555522,
    "city3d:lods": ["0.0", "1.2", "1.3", "2.2"],
    "city3d:co_types": ["Building", "BuildingPart"],
    "city3d:semantic_surfaces": true,
    "cityparquet:version": "0.1.0-draft",
    "proj:code": "EPSG:7415"
  }
}

Connect

https:// access needs the httpfs extension. The file is public, so there is no secret and no credential to configure:

INSTALL httpfs; LOAD httpfs;

Peek at the schema — a DESCRIBE reads only the Parquet footer, not the data:

DESCRIBE SELECT * FROM read_parquet('https://cityparquet.open3d.city/data/3dbag/building.parquet');

Query the whole country in one read

Counting every building in the Netherlands is a single query. Because a COUNT over object_type never touches the geometry columns, DuckDB reads the footer and that one column — a few tens of megabytes out of 16.4 GB:

SELECT count(*) AS buildings
FROM read_parquet('https://cityparquet.open3d.city/data/3dbag/building.parquet')
WHERE object_type = 'Building';
┌──────────┐
│buildings │
├──────────┤
│ 10771547 │   -- 10.8 M buildings out of 21.5 M objects, ~7 s over the network
└──────────┘

The remaining 10.8 M rows are the BuildingParts that carry the geometry.

Bounding-box queries that prune row groups

Every row carries a bbox STRUCT — {xmin, ymin, zmin, xmax, ymax, zmax} — and because those fields feed Parquet’s per-row-group min/max statistics, a bounding-box window skips whole row groups without downloading them. The file has 329 row groups; a window over Delft touches a handful:

SELECT count(*) AS buildings
FROM read_parquet('https://cityparquet.open3d.city/data/3dbag/building.parquet')
WHERE object_type = 'Building'
  AND bbox.xmin < 85000 AND bbox.xmax > 84000     -- overlaps the window in X
  AND bbox.ymin < 447000 AND bbox.ymax > 446000   -- … and in Y
┌──────────┐
│buildings │
├──────────┤
│   2135   │   -- ~2 s, against ~7 s for the national count
└──────────┘

The window is more selective than the national scan even though it reads a wider bbox STRUCT, because the statistics eliminate most row groups before any of their bytes are fetched. That matters far more once geometry is involved: the queries below read geometry_lod2_2, and without the window they would pull gigabytes. The spatial-access design decision argues why the extent is a first-class column, and the benchmark page measures it.

Writing the package with tighter row groups, or ordering rows along a Hilbert curve so spatially close objects land together, makes this sharper still. The cityparquet CLI orders rows along a Hilbert curve by default, and the row-group size is a flag.

The footprint is GeoParquet

At LoD0 a CityParquet building is just a footprint, and a footprint is legal GeoParquet — so DuckDB’s spatial extension decodes it with no help from any city model tooling. The geo footer names the column:

SELECT decode(value)::VARCHAR AS geo
FROM parquet_kv_metadata('https://cityparquet.open3d.city/data/3dbag/building.parquet')
WHERE key::VARCHAR = 'geo';
-- {"version":"1.1.0","primary_column":"geometry_lod0_0",
--  "columns":{"geometry_lod0_0":{"encoding":"WKB","geometry_types":["MultiPolygon Z"],…}}}

Because the footer declares only that column, DuckDB promotes exactly it to a GEOMETRY and leaves the solid LoDs as opaque BLOBs — visible in a DESCRIBE, and reversible with SET enable_geoparquet_conversion=false, which returns all four to BLOB:

geometry_lod0_0   GEOMETRY(…EPSG:7415…)
geometry_lod1_2   BLOB
geometry_lod1_3   BLOB
geometry_lod2_2   BLOB

So the footprint is used like any other geometry column — no ST_GeomFromWKB in sight:

INSTALL spatial; LOAD spatial;

SELECT id,
       ST_GeometryType(geometry_lod0_0)         AS geom_type,
       ROUND(ST_Area(geometry_lod0_0), 1)       AS footprint_m2,
       ST_AsText(ST_Centroid(geometry_lod0_0))  AS centroid
FROM read_parquet('https://cityparquet.open3d.city/data/3dbag/building.parquet')
WHERE object_type = 'Building'
  AND bbox.xmin < 85000 AND bbox.xmax > 84000
  AND bbox.ymin < 447000 AND bbox.ymax > 446000
ORDER BY id LIMIT 4;
┌────────────────────────────────┬──────────────┬──────────────┬─────────────────────────────────┐
│               id               │  geom_type   │ footprint_m2 │            centroid             │
├────────────────────────────────┼──────────────┼──────────────┼─────────────────────────────────┤
│ NL.IMBAG.Pand.0503100000000010 │ MULTIPOLYGON │      48078.9 │ POINT Z (84681.568 445972.369 0)│
│ NL.IMBAG.Pand.0503100000000030 │ MULTIPOLYGON │      15341.2 │ POINT Z (84615.045 446230.551 0)│
│ NL.IMBAG.Pand.0503100000000031 │ MULTIPOLYGON │      10231.9 │ POINT Z (84626.183 446147.505 0)│
│ NL.IMBAG.Pand.0503100000000032 │ MULTIPOLYGON │       7604.8 │ POINT Z (84530.164 446289.164 0)│
└────────────────────────────────┴──────────────┴──────────────┴─────────────────────────────────┘

Any GeoParquet-aware tool — GeoPandas, GDAL/OGR, QGIS — reads that same geometry_lod0_0 column with the correct CRS. The higher-LoD Solid columns are a different matter: a PolyhedralSurface Z is not legal GeoParquet, so the geo footer deliberately does not declare it and strict readers ignore it. This is exactly the boundary where CityParquet stops being a GeoParquet superset — conformant wherever the geometry stays in GeoParquet’s WKB vocabulary, and beyond it only for the 3D solids, which the next section measures.

Measure the 3D solids with three_d

To turn the Solid columns into measurable values, use the three_d extension, available from the DuckDB community registry.

Carve the working set out with the window first. Materialising geometry for the whole country would download gigabytes; the same window as above brings the neighbourhood into memory in a few seconds, and everything after it is local:

INSTALL three_d FROM community;
LOAD three_d;
INSTALL httpfs; LOAD httpfs;

CREATE TEMP TABLE t AS
  SELECT id, parents, object_type, b3_dak_type, b3_volume_lod22,
         geometry_lod2_2 AS g, geometry_properties_lod2_2 AS gp   -- three_d reads the STRUCT directly
  FROM read_parquet('https://cityparquet.open3d.city/data/3dbag/building.parquet')
  WHERE bbox.xmin < 85000 AND bbox.xmax > 84000
    AND bbox.ymin < 447000 AND bbox.ymax > 446000;                -- 4272 rows, ~4 s

ST_3DFromWKB takes the WKB geometry and its properties — in a package the geometry_properties_lod2_2 struct, which three_d reads directly, with no to_json(...) conversion. Join each BuildingPart (which holds the geometry) to its parent Building (which holds the ground-truth attributes), and use ST_3DTryFromWKB so an unsupported geometry yields NULL instead of aborting the query:

CREATE TEMP TABLE parts AS
  SELECT p.id, ST_3DTryFromWKB(p.g, p.gp) AS solid, b.b3_dak_type, b.b3_volume_lod22
  FROM t p JOIN t b ON b.id = p.parents[1]
  WHERE p.g IS NOT NULL;

ST_3DVolume and ST_3DSurfaceArea raise on an open or non-manifold solid, and real reconstructions always contain a few — so check the working set’s health first:

SELECT count(*)                                                     AS parts,
       count(*) FILTER (WHERE ST_3DValidationReport(solid).is_valid) AS valid,
       count(*) FILTER (WHERE ST_3DIsClosed(solid))                 AS closed
FROM parts;
-- parts 2137 · valid 2118 · closed 2134

With measurement guarded to valid solids, the metrics that matter for buildings are short expressions — and they check out against 3DBAG’s own numbers:

SELECT id,
       ROUND(ST_3DVolume(solid), 1)            AS volume_m3,
       ROUND(b3_volume_lod22, 1)               AS volume_3dbag,   -- ground truth
       ROUND(ST_3DFootprintArea(solid), 1)     AS footprint_m2,
       ROUND(ST_3DZMax(solid) - ST_3DZMin(solid), 1) AS height_m
FROM parts
WHERE ST_3DValidationReport(solid).is_valid
ORDER BY id LIMIT 4;
┌──────────────────────────────────┬───────────┬──────────────┬──────────────┬──────────┐
│                id                │ volume_m3 │ volume_3dbag │ footprint_m2 │ height_m │
├──────────────────────────────────┼───────────┼──────────────┼──────────────┼──────────┤
│ NL.IMBAG.Pand.0503100000000010-0 │  617669.9 │     617766.4 │      48078.9 │     12.8 │
│ NL.IMBAG.Pand.0503100000000030-0 │  136462.7 │     136456.7 │      15279.7 │     15.0 │
│ NL.IMBAG.Pand.0503100000000031-0 │   79895.2 │      79910.1 │      10231.7 │     12.1 │
│ NL.IMBAG.Pand.0503100000000032-0 │   43405.7 │      43401.1 │       7604.8 │      6.2 │
└──────────────────────────────────┴───────────┴──────────────┴──────────────┴──────────┘

The volume the kernel computes from the stored WKB tracks b3_volume_lod22 to a fraction of a per mille — a live check that the encoding preserved the geometry losslessly. three_d validates, it never silently repairs; the full function surface (distance queries, transforms, serialisation) is in the three_d documentation.

Roof-type analysis over semantics and geometry

CityParquet keeps CityGML’s semantics queryable, which lets you ask domain questions that need both attributes and geometry. 3DBAG classifies every building’s roof, and that lands as ordinary columns — so a roof-type breakdown, with the flat and sloped roof areas 3DBAG measured, is plain SQL over the windowed Building rows:

SELECT b3_dak_type,
       count(*)                        AS n,
       ROUND(avg(b3_opp_dak_plat), 1)   AS avg_flat_m2,
       ROUND(avg(b3_opp_dak_schuin), 1) AS avg_sloped_m2
FROM read_parquet('https://cityparquet.open3d.city/data/3dbag/building.parquet')
WHERE object_type = 'Building'
  AND bbox.xmin < 85000 AND bbox.xmax > 84000
  AND bbox.ymin < 447000 AND bbox.ymax > 446000
GROUP BY b3_dak_type ORDER BY n DESC;
┌─────────────────────┬──────┬─────────────┬───────────────┐
│     b3_dak_type     │  n   │ avg_flat_m2 │ avg_sloped_m2 │
├─────────────────────┼──────┼─────────────┼───────────────┤
│ slanted             │ 1287 │       163.7 │          91.6 │
│ horizontal          │  737 │        20.7 │           0.7 │
│ multiple horizontal │  110 │       122.0 │           3.0 │
│ unknown             │    1 │     48078.9 │           0.0 │
└─────────────────────┴──────┴─────────────┴───────────────┘

The finer semantics are one level down, in geometry_properties. Its surfaces field is the CityGML surface list verbatim as JSON text (a RoofSurface / WallSurface / GroundSurface type plus 3DBAG’s per-surface b3_azimut, b3_hellingshoek, …), and face_semantics maps each WKB face to an entry in it. Unnest surfaces to profile the semantic surfaces directly — here, the average azimuth and slope of every roof surface in the window:

SELECT s->>'type'                                          AS surface_type,
       count(*)                                            AS n,
       ROUND(avg(TRY_CAST(s->>'b3_azimut'      AS DOUBLE)), 1) AS avg_azimuth,
       ROUND(avg(TRY_CAST(s->>'b3_hellingshoek' AS DOUBLE)), 1) AS avg_slope_deg
FROM (
  SELECT unnest(geometry_properties_lod2_2.surfaces::JSON[]) AS s
  FROM read_parquet('https://cityparquet.open3d.city/data/3dbag/building.parquet')
  WHERE geometry_lod2_2 IS NOT NULL
    AND bbox.xmin < 85000 AND bbox.xmax > 84000
    AND bbox.ymin < 447000 AND bbox.ymax > 446000
)
GROUP BY surface_type ORDER BY n DESC;
┌───────────────┬──────┬─────────────┬───────────────┐
│ surface_type  │  n   │ avg_azimuth │ avg_slope_deg │
├───────────────┼──────┼─────────────┼───────────────┤
│ RoofSurface   │ 6929 │       192.8 │          18.2 │
│ WallSurface   │ 4274 │        NULL │          NULL │  -- walls carry no azimuth/slope
│ GroundSurface │ 2137 │        NULL │          NULL │
└───────────────┴──────┴─────────────┴───────────────┘

Finally, cross the semantic roof label with geometry the kernel measures. Back in the three_d session, group the computed volumes and heights by roof type — the payoff of carrying attributes, semantics, and 3D geometry in one columnar row:

SELECT b3_dak_type,
       count(*)                                     AS n,
       ROUND(SUM(ST_3DVolume(solid)), 0)            AS total_volume_m3,
       ROUND(avg(ST_3DZMax(solid) - ST_3DZMin(solid)), 1) AS avg_height_m
FROM parts
WHERE ST_3DValidationReport(solid).is_valid
GROUP BY b3_dak_type ORDER BY total_volume_m3 DESC;
┌─────────────────────┬──────┬─────────────────┬──────────────┐
│     b3_dak_type     │  n   │ total_volume_m3 │ avg_height_m │
├─────────────────────┼──────┼─────────────────┼──────────────┤
│ slanted             │ 1270 │       2441835   │      9.6     │
│ unknown             │    1 │        617670   │     12.8     │
│ multiple horizontal │  110 │        143469   │      9.5     │
│ horizontal          │  737 │         82818   │      3.4     │
└─────────────────────┴──────┴─────────────────┴──────────────┘

Slanted-roof buildings hold the bulk of the window’s built volume and stand nearly three times as tall on average as the flat-roofed ones — a result that needed the roof semantics, the 3DBAG attributes, and the measured 3D geometry together, all in one scan of one columnar file. The single unknown row is one 48 000 m² industrial complex, a reminder that a national dataset always has outliers.

Convert to CityJSONSeq (and other formats)

Because the cityjson extension also writes city models, COPY … TO turns CityParquet back into CityJSON, CityJSONSeq, or FlatCityBuf — DuckDB as a format converter. It consumes the wide layout directly: a plain SELECT * becomes one multi-LoD CityObject per row, and the STRUCT geometry_properties is reconstructed into CityJSON semantics — no per-LoD projection, no to_json. Keep the window on, so the extract is a neighbourhood rather than a country:

INSTALL cityjson FROM community;
LOAD cityjson;

COPY (
  SELECT * FROM read_parquet('https://cityparquet.open3d.city/data/3dbag/building.parquet')
  WHERE object_type IN ('Building', 'BuildingPart')
    AND bbox.xmin < 85000 AND bbox.xmax > 84000
    AND bbox.ymin < 447000 AND bbox.ymax > 446000
)
TO 'delft_window.city.jsonl'
  (FORMAT cityjsonseq, crs 'https://www.opengis.net/def/crs/EPSG/0/7415');

That writes a streaming CityJSONSeq file — one metadata header line, then one CityJSONFeature per line, each carrying its LoD 1.2 / 1.3 / 2.2 geometries with the roof/wall/ground semantics intact — which reads straight back in to prove the round trip:

SELECT count(*) AS city_objects FROM read_cityjsonseq('delft_window.city.jsonl');
-- 4272, the same 2135 Buildings and 2137 BuildingParts the window selected

CityJSONSeq is the right target for large or streamed data (per-feature vertex pools, low memory); FORMAT cityjson writes a single document with one global vertex pool instead.

Why this is fast

Every query on this page ran against one 16.4 GB file on object storage, from a single machine, in seconds — because CityParquet is columnar and DuckDB exploits it at three levels:

  • Projection pushdown — a query reads only the columns it names. Counting 10.8 M buildings never touches a geometry byte; measuring geometry_lod2_2 never reads geometry_lod1_2.
  • Row-group pruning — the bbox STRUCT and the scalar columns feed Parquet’s footer statistics, so a bounding-box window skips most of the 329 row groups before any of their bytes are fetched. This is what makes the geometry queries possible at all.
  • Read-in-place — httpfs fetches byte ranges over HTTP, so nothing is downloaded whole; the STAC metadata.json describes the package before DuckDB is even involved, and a catalogue of such Items narrows the search to a package first.

The same three properties are what make CityParquet a delivery format for 3D city models, not just a storage one. To reproduce these numbers on your own data, convert a CityJSON source with the reference cityparquet CLI or the format-conversion walkthrough, and point the queries at your package.

Last updated on October 23, 2026

Was this page helpful?