2023-04-28

PostGEESE? Die Spatial-Erweiterung von DuckDB

Max Gabrielsson

Geospatiale Daten werden in modernen Anwendungen und Data-Engineering-Workflows immer wichtiger – von standortbasierten Diensten bis zur Umweltüberwachung.

Es gibt viele gute, spezialisierte Werkzeuge für geospatiale Daten. Geospatiale Fähigkeiten direkt in DuckDB zu integrieren hat trotzdem mehrere Vorteile. Zum einen können Sie geospatiale Daten zusammen mit regulären, unstrukturierten oder Zeitreihendaten über DuckDBs reiches Typsystem und Erweiterungen wie JSON und ICU verarbeiten, transformieren und joinen. Zum anderen lassen sich räumliche Abfragen mit geometrischen Prädikaten und Relationen überraschend gut in SQL ausdrücken – SQL dreht sich ja um Relationen. Dazu kommen DuckDBs weitere Stärken: transaktionale Semantik, hochperformante mehrfädige vektorisierte Ausführung und Verarbeitung größerer Datenmengen als der Speicher.

Deshalb freuen wir uns, dass DuckDB jetzt eine Spatial-Erweiterung mit vielen Features hat, die sich einfach aus der DuckDB-CLI und anderen Clients installieren lässt. Einfach ausführen:

INSTALL spatial;
LOAD spatial;

Und schon kann es losgehen!

Nein, GeoDuck nennen wir es auch nicht, das wäre einfach nur eklig.

Was steckt drin?

Kern der Erweiterung ist ein Typ GEOMETRY nach dem Geometriemodell „Simple Features“ und dazugehörige Funktionen wie ST_Area und ST_Intersects. Außerdem gibt es Methoden zum Lesen und Schreiben geospatialer Formate und zur Umrechnung zwischen Koordinatenreferenzsystemen (Details weiter unten). Volle Konformität mit OGC Simple Feature Access und SQL/MM wollen wir noch nicht zusagen; wer geospatiale Funktionen in anderen Datenbanksystemen wie PostGIS oder SpatiaLite kennt, sollte sich aber sofort zurechtfinden.

Die meisten implementierten Funktionen stützen sich auf das Trio grundlegender geospatialer Bibliotheken: GEOS, GDAL und PROJ – für Algorithmen, Formatumwandlungen und Transformationen von Koordinatenreferenzsystemen. Insbesondere nutzen wir GDAL für Tabellen- und Copy-Funktionen, mit denen sich Tabellen aus über 50 geospatialen Formaten importieren und exportieren lassen (bisher!), darunter die gängigsten wie Shapefiles, GeoJSON, GeoPackage, KML, GML, WKT, WKB usw.

Überzeugen Sie sich selbst:

Zunächst haben wir vor allem Breite geschaffen, indem wir bestehende Bibliotheken umhüllen. Künftig wollen wir mehr Kernfunktionen und Algorithmen nativ implementieren, für mehr Leistung und effizientere Speicherverwaltung.

Als ersten Schritt in diese Richtung gibt es eine Reihe nicht-standardisierter, spezialisierter spaltenorientierter DuckDB-nativer Geometrietypen wie POINT_2D, BOX_2D usw. Sie sollen bessere Kompression und schnellere Ausführung bei etwas weniger Flexibilität bieten; die Arbeit daran ist noch stark experimentell.

Beispielnutzung

Das Folgende zeigt, wie Sie mit der Spatial-Erweiterung mehrere geospatiale Formate lesen und exportieren, Geometrien zwischen Koordinatenreferenzsystemen transformieren und mit räumlichen Eigenschafts- und Prädikatfunktionen arbeiten. Das Beispiel ist etwas konstruiert, soll aber die derzeit verfügbaren Möglichkeiten zeigen. Die verwendeten Datensätze liegen im Repository der Spatial-Erweiterung.

Wir importieren die NYC-Taxi-Fahrten im Parquet-Format sowie die dazugehörigen Taxi-Zonen aus einem Shapefile, über die Tabellenfunktion ST_Read der Spatial-Erweiterung. Diese Taxi-Zonen zerlegen NYC in Polygone, die Regionen darstellen, etwa den Newark Airport. Dann legen wir eine Tabelle für die Fahrten und eine für die Zonen an. ST_Read erzeugt eine Tabelle mit einer Spalte wkb_geometry, die die Geometrie als WKB-Blob (Well-Known Binary) enthält; wir wandeln sie mit ST_GeomFromWKB in den Typ GEOMETRY um.

Das mag viel wirken, wenn Sie das geospatiale Ökosystem nicht kennen – zum Einstieg reicht das aber. Kurz: – Shapefile (.shp, .shx, .dbf) ist ein gängiges Format für Geometrie-Vektordaten und Hilfsmetadaten wie Indizes und Attribute. – WKB (Well Known Binary) ist selbst kein Dateiformat, aber eine gängige Binärkodierung von Vektorgeometrie, etwa in GeoParquet. Es gibt mehrere Varianten; uns geht es vorerst um „standard“-WKB. – GEOMETRY ist ein DuckDB-Typ für ein Simple-Features-Geometrieobjekt. Die Standards modellieren Vektorgeometrie als Punkte, Linestrings, Polygone oder Sammlungen davon. Das ist der Kern-Datentyp der Spatial-Erweiterung, den die meisten Funktionen entgegennehmen und zurückgeben.

INSTALL spatial;
LOAD spatial;

CREATE TABLE rides AS
    SELECT *
    FROM 'yellow_tripdata_2010-01-limit1mil.parquet';

-- Load the NYC taxi zone data from a shapefile using the gdal-based ST_Read function
CREATE TABLE zones AS
    SELECT zone, LocationId, borough, geom
    FROM ST_Read('taxi_zones/taxi_zones.shx');

Vergleichen wir die Fahrstrecke mit der linearen Distanz zwischen Abhol- und Absetzpunkt, um zu sehen, wie effizient die Taxifahrer sind (oder wie schmutzig die Daten, denn manche Differenzen wirken negativ). Wir transformieren die Koordinaten von „WGS84“ (Kennung EPSG:4326), auch einfach Breite/Länge genannt, ins Koordinatenreferenzsystem „NAD83 / New York Long Island ftUS“ (ESRI:102718) – eine Projektion mit geringer Verzerrung um New York. Die Distanz berechnen wir mit ST_Distance. Hier erhalten wir Fuß, weil wir nach NAD83 gewandelt haben; in Meilen (5280 ft/Meile) umrechnen – die Einheit des Fahrten-Datensatzes – ist einfach, damit der Vergleich stimmt.

Fahrten, die kürzer sind als die Luftlinie, sind wahrscheinlich fehlerhaft; mit dieser Abfrage filtern wir schlechte Daten. Die Abfrage nutzt DuckDBs Möglichkeit, auf in demselben Select definierte Spaltenaliase zu verweisen. Ein kleines Beispiel, wie DuckDBs reicher SQL-Dialekt geospatiale Analyse vereinfacht.

CREATE TABLE cleaned_rides AS
    SELECT
        ST_Point(pickup_latitude, pickup_longitude) AS pickup_point,
        ST_Point(dropoff_latitude, dropoff_longitude) AS dropoff_point,
        dropoff_datetime::TIMESTAMP - pickup_datetime::TIMESTAMP AS time,
        trip_distance,
        ST_Distance(
            ST_Transform(pickup_point, 'EPSG:4326', 'ESRI:102718'),
            ST_Transform(dropoff_point, 'EPSG:4326', 'ESRI:102718')) / 5280
            AS aerial_distance,
        trip_distance - aerial_distance AS diff
    FROM rides
    WHERE diff > 0
    ORDER BY diff DESC;

Das ist nicht ganz genau, weil die von uns genutzte Funktion ST_Distance die Erdkrümmung nicht berücksichtigt. Für unsere Zwecke reicht die Näherung. Sphärische und geodätische Distanzberechnungen stehen auf der Roadmap!

Jetzt joinen wir die Taxifahrten mit den Taxi-Zonen, um Start- und Endzone jeder Fahrt zu bekommen. Als Join-Bedingung nutzen wir ST_Within, um zu prüfen, ob ein Abhol- oder Absetzpunkt in einem Taxi-Zonen-Polygon liegt. Wieder müssen wir von WGS84 nach NAD83 transformieren, weil die Zonendaten dieselbe Projektion nutzen. Räumliche Joins wie diese sind das tägliche Brot geospatialer Verarbeitung; Optimierungen (etwa räumliche Indizes) haben wir dafür noch nicht, deshalb nutzen wir im nächsten Schritt nur eine Teilmenge der Daten.

-- Since we don't have spatial indexes yet, use a smaller dataset for the join.
DELETE FROM cleaned_rides WHERE rowid > 5000;

CREATE TABLE joined AS
    SELECT
        pickup_point,
        dropoff_point,
        start_zone.zone AS start_zone,
        end_zone.zone AS end_zone,
        trip_distance,
        time,
    FROM cleaned_rides
    JOIN zones AS start_zone
      ON ST_Within(ST_Transform(pickup_point, 'EPSG:4326', 'ESRI:102718'), start_zone.geom)
    JOIN zones AS end_zone
      ON ST_Within(ST_Transform(dropoff_point, 'EPSG:4326', 'ESRI:102718'), end_zone.geom);

Die gejointe Tabelle können wir mit der GDAL-Copy-Funktion in eine GeoJSONSeq-Datei exportieren und eine GDAL-Layer-Creation-Option übergeben. GeoJSON unterstützt nur eine GEOMETRY pro Datensatz; mit ST_MakeLine kombinieren wir Abhol- und Absetzpunkt zu einer Liniengeometrie. Das Standard-CRS für GeoJSON ist WGS84, die Koordinatenpaare werden aber als Länge/Breite erwartet, daher drehen wir die Geometrie mit ST_FlipCoordinates.

COPY (
    SELECT
        ST_MakeLine(pickup_point, dropoff_point)
            .ST_FlipCoordinates()
            .ST_AsWKB()
            AS wkb_geometry,
        start_zone,
        end_zone,
        time::VARCHAR AS trip_time
    FROM joined)
TO 'joined.geojsonseq'
WITH (
    FORMAT gdal,
    DRIVER 'GeoJSONSeq',
    LAYER_CREATION_OPTIONS 'WRITE_BBOX=YES'
);

{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:52:00" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.789923, 40.643515 ], [ -73.97608, 40.680395 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:35:00" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.776445, 40.645422 ], [ -73.98427, 40.670782 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:45:42" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.776878, 40.645065 ], [ -73.992153, 40.662571 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:36:00" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.788028, 40.641508 ], [ -73.97584, 40.670927 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:47:58" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.781855, 40.644749 ], [ -73.980129, 40.663663 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:32:10" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.787494, 40.641559 ], [ -73.974694, 40.673479 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:36:59" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.790138, 40.643342 ], [ -73.982721, 40.662379 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:32:00" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.786952, 40.641248 ], [ -73.97421, 40.676237 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:33:21" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.783892, 40.648514 ], [ -73.979283, 40.669721 ] ] } }
{ "type": "Feature", "properties": { "start_zone": "JFK Airport", "end_zone": "Park Slope", "trip_time": "00:35:45" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.776643, 40.645272 ], [ -73.978873, 40.66723 ] ] } }

Fertig! Wir haben tabellarische Daten aus Parquet geholt, mit geospatialen Daten aus einem Shapefile kombiniert, die kombinierten Daten bereinigt und analysiert und ins menschenlesbare geospatiale Format geschrieben. Die derzeit unterstützten Funktionen und ihr Umsetzungsstand stehen in der Dokumentation in dieser Tabelle.

Wie geht es weiter?

Bis wir den vollen Funktionsumfang von z. B. PostGIS erreichen, wird es dauern. DuckDBs vektorisiertes Ausführungsmodell und spaltenorientiertes Speicherformat sollten aber eine neue Klasse von Optimierungen für geospatiale Verarbeitung ermöglichen, die wir gerade erst erkunden. Die Leistung räumlicher Joins und Prädikate steht daher weit oben auf der Liste.

Es gibt auch Grenzen unseres Typs GEOMETRY, die wir angehen wollen: keine zusätzlichen Z- und M-Dimensionen und nicht die volle Bandbreite der vom OGC-Standard geforderten Geometrie-Subtypen wie Kurven oder polyedrische Flächen.

Sphärische und ellipsoidische Berechnungen wollen wir bald unterstützen, vielleicht als eigenen Typ GEOGRAPHY.

Wasm-Builds stehen ebenfalls kurz bevor!

Die vollständige Roadmap und den aktuellen Stand finden Sie im GitHub-Repository. Wenn Sie beim Aufbau helfen wollen, melden Sie sich auf GitHub!

Fazit

Die DuckDB-Spatial-Erweiterung ist ein weiterer Schritt, DuckDB zum Schweizer Taschenmesser für Data Engineering und Analytics zu machen. Sie liefert einen flexiblen, vertrauten Typ GEOMETRY, reprojizierbar zwischen Tausenden von Koordinatenreferenzsystemen, plus Import und Export geospatialer Daten aus über 50 Quellen. Alles in einer Erweiterung mit minimalen Laufzeitabhängigkeiten. DuckDB fügt sich so in bestehende GIS-Workflows ein, unabhängig von Format oder Projektion.

Wir sind gespannt, was Sie aus der DuckDB-Spatial-Erweiterung machen. Es ist noch früh, aber wir hoffen, bald mehr teilen zu können! Fragen, Vorschläge, Ideen oder Issues gern auf Discord oder GitHub.