Getting started with modern GIS using DuckDB

Understanding geometries

When working with geospatial functions, you will learn how to work with geometries. In short, these can be points, lines, polygons, or collections of them.

Many databases support spatial function and spatial type to store these geometries, and they are typically prefixed with ST_, which stands for "spatial and temporal". So if you have a coordinate - a point in spatial type - you would convert it with DuckDB by just using the ST_Point from the spatial extension.

An example using the DuckDB CLI :

Copy code
D install spatial;
D load spatial;
D SELECT ST_Point(30, 50) AS location;
┌───────────────┐
│   location    │
│   geometry    │
├───────────────┤
│ POINT (30 50) │
└───────────────┘

File format

The second important point regarding geospatial is the file format. To share geospatial data, there are multiple formats you can work with:

• Vector data represents the discrete features we discussed above, such as points, lines, and polygons (e.g., city locations and roads).
• Raster data is more like a photo and represents continuous information. It consists of a grid of cells (or pixels), and each cell has a value representing something, like temperature, elevation, or colors in a satellite image.

You can find both on the web, but vector data are usually easier to share because of their smaller size. GeoJSON is the easiest one to work with, as you can directly edit, but it's pretty inefficient in terms of size. GeoParquet adoption has been increasing, yet not many tools support it. However, there's no need to worry—DuckDB does! DuckDB has many possibilities to read/write from many of these. You can use the table function FROM ST_Drivers(); to display all supported drivers.

Copy code
D FROM ST_Drivers();
┌────────────────┬──────────────────────────────────────────────────────┬────────────┬──────────┬──────────┬────────────────────────────────────────────────────┐
│   short_name   │                      long_name                       │ can_create │ can_copy │ can_open │                      help_url                      │
│    varchar     │                       varchar                        │  boolean   │ boolean  │ boolean  │                      varchar                       │
├────────────────┼──────────────────────────────────────────────────────┼────────────┼──────────┼──────────┼────────────────────────────────────────────────────┤
│ ESRI Shapefile │ ESRI Shapefile                                       │ true       │ false    │ true     │ https://gdal.org/drivers/vector/shapefile.html     │
│ MapInfo File   │ MapInfo File                                         │ true       │ false    │ true     │ https://gdal.org/drivers/vector/mitab.html         │
│ UK .NTF        │ UK .NTF                                              │ false      │ false    │ true     │ https://gdal.org/drivers/vector/ntf.html           │
│ LVBAG          │ Kadaster LV BAG Extract 2.0                          │ false      │ false    │ true     │ https://gdal.org/drivers/vector/lvbag.html         │
│ S57            │ IHO S-57 (ENC)                                       │ true       │ false    │ true     │ https://gdal.org/drivers/vector/s57.html           │
│ DGN            │ Microstation DGN                                     │ true       │ false    │ true     │ https://gdal.org/drivers/vector/dgn.html           │
│ OGR_VRT        │ VRT - Virtual Datasource                             │ false      │ false    │ true     │ https://gdal.org/drivers/vector/vrt.html           │
│ Memory         │ Memory                                               │ true       │ false    │ true     │                                                    │
│ CSV            │ Comma Separated Value (.csv)                         │ true       │ false    │ true     │ https://gdal.org/drivers/vector/csv.html           │
│ GML            │ Geography Markup Language (GML)                      │ true       │ false    │ true     │ https://gdal.org/drivers/vector/gml.html           │
│ GPX            │ GPX                                                  │ true       │ false    │ true     │ https://gdal.org/drivers/vector/gpx.html           │
│ KML            │ Keyhole Markup Language (KML)                        │ true       │ false    │ true     │ https://gdal.org/drivers/vector/kml.html           │

This makes it super helpful again to convert and join data in a standard format, which is usually a big preparation step for a geospatial project.

As I mentioned in the intro, there are many standards! It's nice to prepare all your geo data with a single tool.

Getting your tools ready

Now that we understand the fundamentals, let's see what we need regarding tooling. You can use SQL and DuckDB, but you still need something to display the data.

A typical stack is to use Python in a notebook environment to render the results directly. We'll use a Google Collab but any notebook environment is good for you, it's just for the simplicity of this tutorial and sharing.

Regarding the data visualization library, again, there are many options. Leafmap is definitely interesting and pretty mature to use. For this blog, however, I'll show you a new kid on the block called Lonboard. It's performant because it doesn't use GeoJSON as an intermediate step to transfer the data to the front end, as many of these tools do. Plus, it supports direct integration with DuckDB.

Now, let's zoom in on the code and the data.

Moving to the cloud

You can create an account in MotherDuck for free. Once sign-up, you can get your access token in the settings menu. Connecting to MotherDuck is simply as changing one line of code, using md: in the duckdb.connect() line - this assumes we set motherduck_token as an environment variable to authenticate to MotherDuck.

Copy code
import duckdb

# Initialize a MotherDuck Connection
con = duckdb.connect('md:')

# Load spatial extension
con.sql('INSTALL spatial;')
con.sql('LOAD spatial;')

# URL filtered with bounding box around France data
poi_url = 'https://api-01.openchargemap.io/v3/poi?client=ocm-data-export&maxresults=100000&compact=true&verbose=false&boundingbox=(51.124,-5.142),(41.342,9.562)'
  
# Ingest the data from the API and create as a table
# Create database if not exist
con.sql("CREATE DATABASE IF NOT EXISTS geo_playground")
con.sql(f"CREATE TABLE IF NOT EXISTS geo_playground.poi_france AS SELECT * FROM read_json_auto('{poi_url}')")

Now, the above query and the rest of the pipeline will leverage cloud computing.

Copy code
# Uploading the current local table to MotherDuck

sql = """ CREATE TABLE IF NOT EXISTS geo_playground.poi_france_display AS
SELECT ID,
       ST_Point(AddressInfo.Longitude, AddressInfo.Latitude) AS geom,
       AddressInfo.Title AS Title
FROM geo_playground.poi_france
WHERE AddressInfo.Latitude IS NOT NULL
  AND AddressInfo.Longitude IS NOT NULL
"""

con.sql(sql)

Finally, as shown below, you can reuse the database and even share it with one line of code .

Copy code
from lonboard import Map, HeatmapLayer

query = con.sql("SELECT * FROM geo_playground.poi_france_display")
layer = HeatmapLayer.from_duckdb(query, con)
m = Map(layer)

Copy code
# Create a MotherDuck Share to share with other MotherDuck users
con.sql("CREATE SHARE IF NOT EXISTS ev_poi_france FROM geo_playground (ACCESS UNRESTRICTED);")

To export into a flat file, for instance GeoJSON, it's a simple as :

Copy code
# Export to geojson
con.sql("COPY geo_playground.poi_france_display TO './ev_poi_france.geojson' WITH (FORMAT GDAL, DRIVER 'GeoJSON');")