DuckDB - A game changer for analyzing OpenStreetMap

Detailed installation instructions can be found in the official documentation. I tend to download the binary from GitHub and put it on my $PATH (e.g. into ~/.local/bin). Thus, I can start duckdb from anywhere.

Once installed, DuckDB provides a minimal CLI that is perfect for testing some queries. It can be launched via duckdb.

$ duckdb
v1.1.0 fa5c2fe15f
Enter ".help" for usage hints.
Connected to a transient in-memory database.
Use ".open FILENAME" to reopen on a persistent database.

Inside this prompt, one could craft queries, perform some simple tests for development even perform some analyses. To quit DuckDB, just press Ctrl + D or type .quit.

Many DuckDB client APIs exist, but I am probably going to use it most from Python and R, as it bridges the gap between convenient development and performance, where the languages have their known limitations.

Of special interest for me in geographic applications is the Spatial extension which enables spatial data analyses in DuckDB.

INSTALL spatial;
LOAD spatial;

This extension provides support for spatial geometry types, spatial scalar functions for creating, manipulating, aggregating and analyzing spatial objects, as well as read / write functionality for common spatial formats such as GeoPackages. Also, OSM data can be queried directly from .osm.pbf files.

Converting the data for convenient use in DuckDB is made simple by the tool QuackOSM. It encodes all the logic for converting the OSM-style relational data storage of nodes, ways (open & closed), and relations into the more familiar points, lines, and polygons. Running quackosm on a file converts it to a geoparquet file.

pip install quackosm
quackosm germany-latest.osm.pbf

Although, I have recently been working on a Pull Request to allow for direct import into a native DuckDB database, too, data can easily be imported from the parquet files:

CREATE TABLE osm AS
SELECT * REPLACE (st_geomfromwkb(geometry) as geometry)
FROM read_parquet('files/germany-latest_nofilter_noclip_compact.geoparquet');

The feature that makes OSM especially valuable are the diverse tags, key=value pairs that are used to annotate features in the OSM database. DuckDB makes these tags accessible inside one column as a Map data type. This closely mimics how attribute data are stored in OSM tags.

SELECT tags FROM osm WHERE tags IS NOT NULL LIMIT 10;
/*┌───────────────────────────────────────────────────────────────────────────────────────────────────────┐
  │                                                 tags                                                  │
  │                                         map(varchar, varchar)                                         │
  ├───────────────────────────────────────────────────────────────────────────────────────────────────────┤
  │ {description=Weltkriegsdenkmal, historic=memorial, memorial=war_memorial}                             │
  │ {TMC:cid_58:tabcd_1:Class=Point, TMC:cid_58:tabcd_1:Direction=positive, TMC:cid_58:tabcd_1:LCLversi…  │
  │ {TMC:cid_58:tabcd_1:Class=Point, TMC:cid_58:tabcd_1:Direction=positive, TMC:cid_58:tabcd_1:LCLversi…  │
  │ {bus=yes, name=Michaeliskirche, network=HVV, public_transport=stop_position}                          │
  │ {bicycle=yes, button_operated=yes, crossing=traffic_signals, highway=crossing, tactile_paving=no, t…  │
  │ {TMC:cid_58:tabcd_1:Class=Point, TMC:cid_58:tabcd_1:Direction=negative, TMC:cid_58:tabcd_1:LCLversi…  │
  │ {crossing:activation=local, crossing:barrier=no, crossing:chicane=no, crossing:light=yes, crossing:…  │
  │ {TMC:cid_58:tabcd_1:Class=Point, TMC:cid_58:tabcd_1:Direction=positive, TMC:cid_58:tabcd_1:LCLversi…  │
  │ {crossing=no, highway=traffic_signals, traffic_signals=signal, traffic_signals:direction=forward}     │
  │ {TMC:cid_58:tabcd_1:Class=Point, TMC:cid_58:tabcd_1:Direction=both, TMC:cid_58:tabcd_1:LCLversion=8…  │
  ├───────────────────────────────────────────────────────────────────────────────────────────────────────┤
  │                                                10 rows                                                │
  └───────────────────────────────────────────────────────────────────────────────────────────────────────┘
  Run Time (s): real 0.003 user 0.009330 sys 0.002276*/

In DuckDB they can be accessed with similar syntax to Python's Pandas or R's DataFrames. This makes them very easy to query. The following example shows how to count the number of OSM objects in Germany are tagged with leisure=park.

SELECT count(*) FROM osm WHERE list_contains(tags['leisure'], 'park');
/*┌──────────────┐
  │ count_star() │
  │    int64     │
  ├──────────────┤
  │        44931 │
  └──────────────┘
  Run Time (s): real 5.900 user 44.860880 sys 0.090244*/

Of course, the true power of the combination of OSM and DuckDB in my view lies in the ability to perform spatial analyses. Similar to PostgreSQL with the PostGIS extension, the spatial extension allows for spatial operations on the geometry column in the database. For example, let's find out, what the total area in m² is of all buildings in Germany (at least the ones mapped in OSM). But first, let's check how many buildings are in the dataset.

  SELECT COUNT(*) FROM osm WHERE tags['building'] != [];
  /*┌──────────────┐
    │ count_star() │
    │    int64     │
    ├──────────────┤
    │     37451546 │
    └──────────────┘ */
  SELECT sum(ST_Area(ST_Transform(geom, 'EPSG:4326', 'EPSG:3035')))
  FROM osm
  WHERE tags['building'] != [] AND ST_GeometryType(geometry)::VARCHAR LIKE 'POLYGON';
  /*┌────────────────────────────────────────────────────────────────┐
    │ sum(st_area(st_transform(geometry, 'EPSG:4326', 'EPSG:3035'))) │
    │                             double                             │
    ├────────────────────────────────────────────────────────────────┤
    │                                              9664382739.825466 │
    └────────────────────────────────────────────────────────────────┘*/

What still blows my mind is how fast DuckDB handles these kinds of queries. Afterall, the query first had to find all the buildings using WHERE inside nested tags and filter out all geometries but polygons, then ST_Transform the geometries, then calculate the ST_Area per building, and last, add them up. All that on a 5 year old laptop in under 5 minutes seconds. This just makes me gitty. It enables large scale spatial analyses without the overhead of a server-hosted PostGIS database.

Currently, the spatial extension provides most functionality for doing basic spatial operations and manipulations, but it will surely evolve in the future. The ecosystem around spatial data is promising. One honorable mention here must be the H3 community extension.