DuckDB Spatial Extension

🚧 WORK IN PROGRESS 🚧

Table of contents

• DuckDB Spatial Extension
• What is this?
• Example Usage
• How do I get it?
  • Through the DuckDB CLI
  • Development builds
  • Building from source
• Limitations and Roadmap
• Internals and technical details
  • Multi-tiered Geometry Type System
  • Per-thread Arena Allocation for Geometry Objects
  • Embedded PROJ Database
  • Embedded GDAL based Input/Output Functions
• Supported Functions

What is this?

This is a prototype of a geospatial extension for DuckDB that adds support for working with spatial data and functions in the form of a GEOMETRY type based on the the "Simple Features" geometry model, as well as non-standard specialized columnar DuckDB native geometry types that provide better compression and faster execution in exchange for flexibility.

Please note that this extension is still in a very early stage of development, and the internal storage format for the geometry types may change indiscriminately between commits. We are actively working on it, and we welcome both contributions and feedback. Please see the function table or the roadmap entries for the current implementation status.

If you or your organization have any interest in sponsoring development of this extension, or have any particular use cases you'd like to see prioritized or supported, please consider sponsoring the DuckDB foundation or contacting DuckDB Labs for commercial support.

Example Usage

The following is a slightly contrived example of how you can use this extension to read and export multiple geospatial data formats, transform geometries between different coordinate reference systems and work with spatial property and predicate functions.

Let's start by loading the spatial extension and the parquet extension so we can import the NYC taxi ride data in parquet format, and the accompanying taxi zone data from a shapefile, using the spatial ST_Read gdal-based table function. We then create a table for the rides and a table for the zones. Note that ST_Read produces a table with a wkb_geometry column that contains the geometry data encoded as a WKB (Well-Known Binary) blob, which we then convert to the GEOMETRY type using the ST_GeomFromWKB function.

LOAD spatial;
LOAD parquet;

CREATE TABLE rides AS SELECT * 
FROM './spatial/test/data/nyc_taxi/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, ST_GeomFromWKB(wkb_geometry) AS geom 
FROM st_read('./spatial/test/data/nyc_taxi/taxi_zones/taxi_zones.shx');

Let's compare the trip distance to the linear distance between the pickup and dropoff points to figure out how efficient the taxi drivers are (or how dirty the data is, since some diffs seem to be negative). We transform the coordinates from WGS84 (EPSG:4326) (lat/lon) to the NAD83 / New York Long Island ftUS (ESRI:102718) projection and calculate the distance using the ST_Distance function, which in this case gives the distance in feet, which we then convert to miles (5280 ft/mile). Trips that are smaller than the aerial distance are likely to be erroneous, so we use this query to filter out some bad data. Although this is not entirely accurate since the distance we use does not take into account the curvature of the earth, but it is a good enough approximation for our purposes.

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;

Now lets join the taxi rides with the taxi zones to get the start and end zone for each ride using the ST_Within function to check if a point is within a polygon. Again we need to transform the coordinates from WGS84 to the NAD83 since the taxi zone data also use that projection.

-- Since we dont have spatial indexes yet, use a smaller dataset for the following example.
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);

We can export the joined table to a GeoJSONSeq file using the GDAL copy function, passing in a GDAL layer creation option. Since GeoJSON only supports a single geometry per feature, we can use the ST_MakeLine function to combine the pickup and dropoff points into a single line geometry. The default coordinate reference system for GeoJSON is WGS84, but the coordinates are expected to be in longitude/latitude, so we need to flip the geometry using the ST_FlipCoordinates function.

COPY (
    SELECT 
        ST_AsWKB(ST_FlipCoordinates(ST_MakeLine(pickup_point, dropoff_point))) 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 ] ] } }

How do I get it?

Through the DuckDB CLI

You can install the extension for DuckDB through the DuckDB CLI like you would do for other first party extensions. Simply execute: INSTALL spatial; LOAD spatial!

Development builds

You can also grab the lastest builds directly from the CI runs or the release page here on GitHub and install manually.

Once you have downloaded the extension for your platform, you need to:

• Unzip the archive
• Start duckdb with the -unsigned flag to allow loading unsigned extensions. (This won't be neccessary in the future)
• Run INSTALL 'absolute/or/relative/path/to/the/unzipped/extension';
• The extension is now installed, you can now load it with LOAD spatial; whenever you want to use it.

You can also build the extension yourself following the instructions below.

Building from source

This extension is based on the DuckDB extension template.

Dependencies

You need a recent version of CMake (3.20) and a C++11 compatible compiler. You also need OpenSSL on your system. On ubuntu you can install it with sudo apt install libssl-dev, on macOS you can install it with brew install openssl. Note that brew installs openssl in a non-standard location, so you may need to set a OPENSSL_ROOT_DIR=$(brew --prefix openssl) environment variable when building.

We bundle all the other required dependencies in the third_party directory, which should be automatically built and statically linked into the extension. This may take some time the first time you build, but subsequent builds should be much faster.

We also highly recommend that you install Ninja which you can select when building by setting the GEN=ninja environment variable.

git clone --recurse-submodules https://github.com/duckdblabs/duckdb_spatial
cd duckdb_spatial
make debug

You can then invoke the built DuckDB (with the extension statically linked)

./build/debug/duckdb

Please see the Makefile for more options, or the extension template documentation for more details.

Limitations and Roadmap

The main limitations of this extension currently are:

• No support for higher-dimensional geometries (XYZ, XYZM, XYM)
• No support for spherical geometry (e.g. lat/lon coordinates)
• No support for spatial indexing.

These are all things that we want to address eventually, have a look at the open issues and roadmap entries for more details. Please feel free to also open an issue if you have a specific use case that you would like to see supported.

Internals and technical details

Multi-tiered Geometry Type System

This extension implements 5 different geometry types. Like almost all geospatial databases we include a GEOMETRY type that (at least strives) to follow the Simple Features geometry model. This includes support for the standard subtypes, such as POINT, LINESTRING, POLYGON, MULTIPOINT, MULTILINESTRING, MULTIPOLYGON, GEOMETRYCOLLECTION that we all know and love, internally represented in a row-wise fashion on top of DuckDB BLOBs. The internal binary format is very similar to the one used by PostGIS - basically double aligned WKB, and we may eventually look into enforcing the format to be properly compatible with PostGIS (which may be useful for the PostGIS scanner extension). Most functions that are implemented for this type uses the GEOS library, which is a battle-tested C++ port of the famous JTS library, to perform the actual operations on the geometries.

While having a flexible and dynamic GEOMETRY type is great to have, it is comparatively rare to work with columns containing mixed-geometries after the initial import and cleanup step. In fact, in most OLAP use cases you will probably only have a single geometry type in a table, and in those cases you're paying the performance cost to de/serialize and branch on the internal geometry format unneccessarily, i.e. you're paying for flexibility you're not using. For those cases we implement a set of non-standard DuckDB "native" geometry types, POINT_2D, LINESTRING_2D, POLYGON_2D, and BOX_2D. These types are built on DuckDBs STRUCT and LIST types, and are stored in a columnar fashion with the coordinate dimensions stored in separate "vectors". This makes it possible to leverage DuckDB's per-column statistics, compress much more efficiently and perform spatial operations on these geometries without having to de/serialize them first. Storing the coordinate dimensions into separate vectors also allows casting and converting between geometries with multiple different dimensions basically for free. And if you truly need to mix a couple of different geometry types, you can always use a DuckDB UNION type.

For now only a small amount of spatial functions are overloaded for these native types, but since they can be implicitly cast to GEOMETRY you can always use any of the functions that are implemented for GEOMETRY on them as well in the meantime while we work on adding more (although with a de/serialization penalty).

This extension also includes a WKB_BLOB type as an alias for BLOB that is used to indicate that the blob contains valid WKB encoded geometry.

Per-thread Arena Allocation for Geometry Objects

When materializing the GEOMETRY type objects from the internal binary format we use per-thread arena allocation backed by DuckDB's buffer manager to amortize the contention and performance cost of performing lots of small heap allocations and frees, which allows us to utilizes DuckDB's multi-threaded vectorized out-of-core execution fully. While most spatial functions are implemented by wrapping GEOS, which requires an extra copy/allocation step anyway, the plan is to incrementally implementat our own versions of the simpler functions that can operate directly on our own GEOMETRY representation in order to greatly accelerate geospatial processing.

Embedded PROJ Database

PROJ is a generic coordinate transformation library that transforms geospatial coordinates from one projected coordinate reference system (CRS) to another. This extension experiments with including an embedded version of the PROJ database inside the extension binary itself so that you don't have to worry about installing the PROJ library separately. This also opens up the possibility to use this functionality in WASM.

Embedded GDAL based Input/Output Functions

GDAL is a translator library for raster and vector geospatial data formats. This extension includes and exposes a subset of the GDAL vector drivers through the ST_Read and COPY ... TO ... WITH (FORMAT GDAL) table and copy functions respectively to read and write geometry data from and to a variety of file formats as if they were DuckDB tables. We currently support the over 50 GDAL formats - check for yourself by running

ST_Read also supports limited support for predicate pushdown and spatial filtering (if the underlying GDAL driver supports it), but column pruning (projection pushdown) while technically feasible is not yet implemented. ST_Read also allows using GDAL's virtual filesystem abstractions to read data from remote sources such as S3, or from compressed archives such as zip files.

Note: This functionality does not make full use of parallelism due to GDAL not being thread-safe, so you should expect this to be slower than using e.g. the DuckDB Parquet extension to read the same GeoParquet or DuckDBs native csv reader to read csv files. Once we implement support for reading more vector formats natively through this extension (e.g. GeoJSON, GeoBuf, ShapeFile) we will probably split this entire GDAL part into a separate extension.

Supported Functions

🧭 - GEOS - functions that are implemented using the GEOS library

🦆 - DuckDB - functions that are implemented natively in this extension that are capable of operating directly on the DuckDB types

🔄 - CAST(GEOMETRY) - functions that are supported by implicitly casting to GEOMETRY and then using the GEOMETRY implementation

We are actively working on implementing more functions, and will update this table as we go. Again, please feel free to open an issue if there is a particular function you would like to see implemented. Contributions are also welcome!