Simon Willison's Weblog: spatialite

TG: Polygon indexing

2023-09-23T04:32:14+00:00

TG is a brand new geospatial library by Josh Baker, author of the Tile38 in-memory spatial server (kind of a geospatial Redis). TG is written in pure C and delivered as a single C file, reminiscent of the SQLite amalgamation.

TG looks really interesting. It implements almost the exact subset of geospatial functionality that I find most useful: point-in-polygon, intersect, WKT, WKB, and GeoJSON—all with no additional dependencies.

The most interesting thing about it is the way it handles indexing. In this documentation Josh describes two approaches he uses to speeding up point-in-polygon and intersection using a novel approach that goes beyond the usual RTree implementation.

I think this could make the basis of a really useful SQLite extension—a lighter-weight alternative to SpatiaLite.

Via Hacker News

Tags: c, geospatial, spatialite, sqlite, geojson, tg

Weeknotes: AI hacking and a SpatiaLite tutorial

2023-01-15T19:45:39+00:00

Short weeknotes this time because the key things I worked on have already been covered here:

How to implement Q&A against your documentation with GPT3, embeddings and Datasette described my experiment using custom SQL functions in Datasette to provide a Q&A interface against an existing corpus of documents. This resulted in a flurry of software releases:
- datasette-openai - custom SQL functions for querying the OpenAI APIs from Datasette
- datasette-faiss - a plugin that builds a FAISS index for fast vector similarity search
- datasette-cookies-for-magic-parameters - a supporting plugin that helps store secrets such as API keys in cookies so they won't be transmitted in query strings
- openai-to-sqlite - a command-line tool for storing OpenAI API responses in a SQLite database - useful for calculating embeddings
Datasette 0.64, with a warning about SpatiaLite - I released Datasette 0.64 with a new --setting default_allow_sql off setting, due to the discovery that some of the SQL functions added by the SpatiaLite module can be used to crash the Datasette instance.
Building a location to time zone API with SpatiaLite is the latest in the official series of tutorials for Datasette. It demonstrates how SpatiaLite can be used to safely build and deploy a geospatial API.

Releases this week

datasette-openai: 0.2 - (4 releases total) - 2023-01-14
SQL functions for calling OpenAI APIs
openai-to-sqlite: 0.2 - (2 releases total) - 2023-01-13
Save OpenAI API results to a SQLite database
datasette-cookies-for-magic-parameters: 0.1.2 - (3 releases total) - 2023-01-12
UI for setting cookies to populate magic parameters
git-history: 0.7a0 - (10 releases total) - 2023-01-11
Tools for analyzing Git history using SQLite
datasette: 0.64.1 - (124 releases total) - 2023-01-11
An open source multi-tool for exploring and publishing data
datasette-faiss: 0.1a0 - 2023-01-11
Maintain a FAISS index for specified Datasette tables
json-to-files: 0.1 - 2023-01-10
Create separate files on disk based on a JSON object
datasette-auth-passwords: 1.1 - (12 releases total) - 2023-01-09
Datasette plugin for authentication using passwords
datasette-publish-fly: 1.3 - (9 releases total) - 2023-01-09
Datasette plugin for publishing data using Fly
shapefile-to-sqlite: 0.4.2 - (6 releases total) - 2023-01-08
Load shapefiles into a SQLite (optionally SpatiaLite) database

TIL this week

Tags: spatialite, ai, datasette, weeknotes, gpt-3, openai, generative-ai, vector-search, llms

Datasette 0.64, with a warning about SpatiaLite

2023-01-09T21:22:51+00:00

I release Datasette 0.64 this morning. This release is mainly a response to the realization that it's not safe to run Datasette with the SpatiaLite extension loaded if that Datasette instance is configured to enable arbitrary SQL queries from untrusted users.

Here are the release notes quoted in full:

Datasette now strongly recommends against allowing arbitrary SQL queries if you are using SpatiaLite. SpatiaLite includes SQL functions that could cause the Datasette server to crash. See SpatiaLite for more details.

New default_allow_sql setting, providing an easier way to disable all arbitrary SQL execution by end users: datasette --setting default_allow_sql off. See also Controlling the ability to execute arbitrary SQL. (#1409)

Building a location to time zone API with SpatiaLite is a new Datasette tutorial showing how to safely use SpatiaLite to create a location to time zone API.

New documentation about how to debug problems loading SQLite extensions. The error message shown when an extension cannot be loaded has also been improved. (#1979)

Fixed an accessibility issue: the <select> elements in the table filter form now show an outline when they are currently focused. (#1771)

The problem with SpatiaLite

Datasette allows arbitrary SQL execution as a core feature. It takes a bunch of steps to provide this safely: database connections are opened in read-only mode, it imposes a strict time limit on SQL queries and Datasette is designed to be run in containers for a further layer of protection.

SQLite itself is an excellent platform for this feature: it has a set of default functionality that supports this well, protected by a legendarily thorough test suite.

SpatiaLite is a long running third-party extension for SQLite that adds a bewildering array of additional functionality to SQLite - much of it around GIS, but with a whole host of extras as well. It includes debugging routines, XML parsers and even it's own implementation of stored procedures!

Unfortunately, not all of this functionality is safe to expose to untrusted queries - even for databases that have been opened in read-only mode.

After identifying functions which could crash the Datasette instance, I decided that Datasette should make a strong recommendation not to expose SpatiaLite in an unprotected manner.

In addition to the new documentation, I also added a feature I've been planning for a while: a simple setting for disabling arbitrary SQL queries entirely:

datasette --setting default_allow_sql off

Prior to 0.64 you could achieve the same thing by adding the following line to your metadata.json file:

{
    "allow_sql": false
}

Or in metadata.yml:

allow_sql: false

The new setting achieves the same thing, but is more obvious and can be easily applied even for Datasette instances that don't use metadata.

A new SpatiaLite tutorial

The documentation now recommends running SpatiaLite instances with pre-approved SQL implemented using Datasette's canned queries feature.

To help clarify how this works, I decided to publish a new entry in the official series of Datasette tutorials:

Building a location to time zone API with SpatiaLite

This is an updated version of a tutorial I first wrote back in 2017.

The new tutorial now includes material on Chris Amico's datasette-geojson-map plugin, SpatiaLite point-in-polygon queries, polygon intersection queries, spatial indexes and how to use the simplify() function to reduce huge polygons down to a size that is more practical to display on a map.

I'm really happy with this new tutorial. Not only does it show a safe way to run SpatiaLite, but it also illustrates a powerful pattern for using Datasette to create and deploy custom APIs.

The resulting API can be accessed here:

https://timezones.datasette.io/timezones

It's hosted on Fly, using their $1.94/month instance size with 256MB of RAM - easily powerful enough to host this class of application.

I also updated the datasette-publish-fly plugin to make it easier to deploy instances with SQL execution disabled, see the 1.3 release notes.

Tags: security, spatialite, datasette, annotated-release-notes

Spatialite Speed Test

2021-04-04T16:28:42+00:00

Spatialite Speed Test

Part of an excellent series of posts about SpatiaLite from 2012—here John C. Zastrow reports on running polygon intersection queries against a 1.9GB database file in 40 seconds without an index and 0.186 seconds using the SpatialIndex virtual table mechanism.

Tags: geospatial, spatialite, sqlite

Drawing shapes on a map to query a SpatiaLite database (and other weeknotes)

2021-01-24T06:50:29+00:00

This week I built a Datasette plugin that lets you query a database by drawing shapes on a map!

datasette-leaflet-freedraw

SpatiaLite is a SQLite extension that adds a plethora of geospatial querying features. Datasette works well with SpatiaLite already, but every now and then I dip in to a geospatial project to see if there are any neat ways I can extend Datasette's spatial querying support.

datasette-leaflet-freedraw is my new plugin which brings the excellent FreeDraw Leaflet module to Datasette.

It's a pure-JavaScript plugin which looks for Datasette SQL parameters called :freedraw and enhances them with an interactive map. You can then draw on the map to populate the form field with a GeoJSON MultiPolygon representing the shape that you drew. Submit that to the server to excute it as part of a query.

This is best illustrated with a demo. Here's a SQL query that let you draw a shape to find California protected areas (parks, nature reserves and suchlike) that intersect with the polygon that you draw:

select
  AsGeoJSON(geometry), *
from
  CPAD_2020a_SuperUnits
where
  Intersects(GeomFromGeoJSON(:freedraw), geometry) = 1
  and CPAD_2020a_SuperUnits.rowid in (
    select
      rowid
    from
      SpatialIndex
    where
      f_table_name = 'CPAD_2020a_SuperUnits'
      and search_frame = GeomFromGeoJSON(:freedraw)
  )

This uses the California Protected Areas Database by GreenInfo Network, discussed here previously.

Try the SQL query here.

(That from SpatialIndex sub-select at the end is a pattern for taking advantage of SpatiaLite spatial indexes.)

The plugin itself is pretty simple: it scans for freedraw input fields and enhances them with a map interface that writes GeoJSON back to the form field.

You can then use SpatiaLite functions such as GeomFromGeoJSON() to transform that user input into something that can be used in a query.

Another fun query: let's add WHERE PARK_NAME like '%mini%' to the query to filter for just the "mini parks" in San Francisco:

select
  AsGeoJSON(geometry), *
from
  CPAD_2020a_SuperUnits
where
  PARK_NAME like '%mini%' and
  Intersects(GeomFromGeoJSON(:freedraw), geometry) = 1
  and CPAD_2020a_SuperUnits.rowid in (
    select
      rowid
    from
      SpatialIndex
    where
      f_table_name = 'CPAD_2020a_SuperUnits'
      and search_frame = GeomFromGeoJSON(:freedraw)
  )

Here are the mini parks for San Francisco and for Oakland.

This demo uses the user's shapes as part of a SELECT query, but it's also possible to use Datasette's Writable canned queries to create forms that will save the shapes to a database using an INSERT or UPDATE query.

If you want to try this plugin out note that it currently requires the Datasette 0.54a alpha. This is because it takes advantage of the new <script type="module"> support I added last week. Update 25th January 2021: Datasette 0.54 has now been released.

Other releases this week

I recorded my talk for FOSDEM - and in doing so I broke my unproductive streak and pushed out releases for a whole bunch of Datasette projects - partly through fixing bugs I spotted while trying to record the talk!

Here's the full list, with some annotations.

datasette-auth-github: 0.13 - 2021-01-24
Datasette plugin that authenticates users against GitHub. This is a big release: I finally finished upgrading the plugin to work against Datasette's new-ish authentication and permissions mechanism.
datasette-leaflet-freedraw: 0.1.4 - 2021-01-24
Draw polygons on maps in Datasette. See above.
datasette-cluster-map: 0.16 - 2021-01-23
Datasette plugin that shows a map for any data with latitude/longitude columns. Fixed some bugs related to the display of marker popups that included foreign key references to other tables.
asgi-csrf: 0.8 - 2021-01-23
ASGI middleware for protecting against CSRF attacks. I fixed a broken test caused by my use of an undocumented and unstable API in httpx - details here. I also added a new always_set_cookie=True option for cases where you always want to ensure a CSRF cookie has been set, see #16.
datasette-leaflet-geojson: 0.7 - 2021-01-20
Datasette plugin that replaces any GeoJSON column values with a Leaflet map. I fixed a long-standing point of concern: this plugin renders a map for every row on a page containing GeoJSON, which can potentially mean hundreds of maps for queries that return may results. The plugin now renders just the first ten maps and provides a grey clickable "Click to show map" widget for anything past the first ten.
datasette: 0.54a0 - 2021-01-19
Lots of stuff in this one, see the release notes. I pushed out an alpha so I could start using the JavaScript modules stuff in datasette-leaflet-freedraw.
sqlite-transform: 0.3.1 - 2021-01-18
Tool for running transformations on columns in a SQLite database. I improved the --help.
swarm-to-sqlite: 0.3.3 - 2021-01-18
Create a SQLite database containing your checkin history from Foursquare Swarm. I fixed a bug caused by new fields showing up in the Swarm JSON output.
sqlite-utils: 3.3 - 2021-01-18
Python CLI utility and library for manipulating SQLite databases. I added the alter=True argument to the .m2m() method, which I needed to fix the bug in swarm-to-sqlite.

TIL this week

Tags: geospatial, projects, spatialite, sqlite, datasette, weeknotes, leaflet

California Protected Areas Database in Datasette

2020-08-21T23:15:58+00:00

California Protected Areas Database in Datasette

I built this yesterday: it’s a Datasette interface on top of the CPAD 2020 GIS database of protected areas in California maintained by GreenInfo Network. This was a useful excuse to build a GitHub Actions flow that builds a SpatiaLite database using my shapefile-to-sqlite tool, and I fixed a few bugs in my datasette-leaflet-geojson plugin as well.

Via calands-datasette on GitHub

Tags: geospatial, projects, shapefiles, spatialite, california, datasette, github-actions, leaflet

Things I learned about shapefiles building shapefile-to-sqlite

2020-02-19T05:25:58+00:00

The latest in my series of x-to-sqlite tools is shapefile-to-sqlite. I learned a whole bunch of things about the ESRI shapefile format while building it.

Governments really love ESRI shapefiles. There is a huge amount of interesting geospatial data made available in the format - 4,614 on Data.gov!

shapefile-to-sqlite

shapefile-to-sqlite loads the data from these files into a SQLite database, turning geometry properties into database columns and the geometry itself into a blob of GeoJSON. Let's try it out on a shapefile containing the boundaries of US national parks.

$ wget http://nrdata.nps.gov/programs/lands/nps_boundary.zip
...
Saving to: ‘nps_boundary.zip’
nps_boundary.zip                           100%[=====================================================================================>]  12.61M   705KB/s    in 22s     
2020-02-18 19:59:22 (597 KB/s) - ‘nps_boundary.zip’ saved [13227561/13227561]

$ unzip nps_boundary.zip 
Archive:  nps_boundary.zip
inflating: temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.xml  
inflating: temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.dbf  
inflating: temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.prj  
inflating: temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shp  
inflating: temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shx

$ shapefile-to-sqlite nps.db temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shp
temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shp
[####################################]  100%

$ datasette nps.db
Serve! files=('nps.db',) (immutables=()) on port 8003
INFO:     Started server process [33534]
INFO:     Waiting for application startup.
INFO:     Application startup complete.
INFO:     Uvicorn running on http://127.0.0.1:8001 (Press CTRL+C to quit)

I recommend installing the datasette-leaflet-geojson plugin, which will turn any column containing GeoJSON into a Leaflet map.

If you've installed SpatiaLite (installation instructions here) you can use the --spatialite option to instead store the geometry in a SpatiaLite column, unlocking a bewildering array of SQL geometry functions.

$ shapefile-to-sqlite nps.db temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shp --spatialite --table=nps-spatialite
temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.shp
[##################################--]   94%  00:00:00

I deployed a copy of the resulting database using Cloud Run:

$ datasette publish cloudrun nps.db \
    --service national-parks \
    --title "National Parks" \
    --source_url="https://catalog.data.gov/dataset/national-parks" \
    --source="data.gov" \
    --spatialite \
    --install=datasette-leaflet-geojson \
    --install=datasette-render-binary \
    --extra-options="--config max_returned_rows:5"

I used max_returned_rows:5 there because these geometrries are pretty big - without it a page with 100 rows on it can return over 90MB of HTML!

You can browse the GeoJSON version of the table here and the SpatiaLite version here.

The SpatiaLite version defaults to rendering each geometry as an ugly binary blob. You can convert them to GeoJSON for compatibility with datasette-leaflet-geojson using the SpatiaLite AsGeoJSON() function:

select id, UNIT_NAME, AsGeoJSON(geometry)
from [nps-spatialite]

Here's the result of that query running against the demo.

Understanding shapefiles

The most confusing thing about shapefiles is that they aren't a single file. A shapefile comes as a minimum of three files: foo.shp containing geometries, foo.shx containing an index into those geometries (really more of an implementation detail) and foo.dbf contains key/value properties for each geometry.

They often come bundled with other files too. foo.prj is a WKT projection for the data for example. Wikipedia lists a whole bunch of other possibilities.

As a result, shapefiles are usually distributed as a zip file. Some shapefile libraries can even read directly from a zip.

The GeoJSON format was designed as a modern alternative to shapefiles, so understanding GeoJSON really helps in understanding shapefiles. In particular the GeoJSON geometry types: Point, LineString, MultiLineString, Polygon and MultiPolygon match how shapefile geometries work.

An important detail in shapefiles is that data in the .shp and .dbf files is matched by array index - so the first geometry can be considered as having ID=0, the second ID=1 and so on.

You can read the properties from the .dbf file using the dbfread Python module like this:

$ ipython
In [1]: import dbfread
In [2]: db = dbfread.DBF("temp/Current_Shapes/Data_Store/06-06-12_Posting/nps_boundary.dbf")
In [3]: next(iter(db))
Out[3]: 
OrderedDict([('UNIT_TYPE', 'Park'),
            ('STATE', ''),
            ('REGION', 'NC'),
            ('UNIT_CODE', 'NACC'),
            ('UNIT_NAME', 'West Potomac Park'),
            ('DATE_EDIT', None),
            ('GIS_NOTES', ''),
            ('CREATED_BY', 'Legacy'),
            ('METADATA', ''),
            ('PARKNAME', '')])

Reading shapefiles in Python

I'm a big fan of the Shapely Python library, so I was delighted to see that Sean Gillies, creator of Shapely, also created a library for reading and writing shapefiles: Fiona.

GIS with Python, Shapely, and Fiona by Tom MacWright was particularly useful for figuring this out. I like how he wrote that post in 2012 but added a note in 2017 that it's still his recommended way of getting started with GIS in Python.

Projections

The trickiest part of working with any GIS data is always figuring out how to deal with projections.

GeoJSON attempts to standardize on WGS 84, otherwise known as the latitude/longitude model used by GPS. But... shapefiles frequently use something else. The Santa Clara county parks shapefiles for example use EPSG:2227, also known as California zone 3.

(Fun fact: ESPG stands for European Petroleum Survey Group, a now defunct oil industry group that today lives on only as a database of projected coordinate systems.)

I spent quite a while thinking about how to best handle projections. In the end I decided that I'd follow GeoJSON's lead and attempt to convert everything to WGS 84, but allow users to skip that behaviour using --crs=keep or to specify an alternative projection to convert to with --crs=epsg:2227 or similar.

SpatiaLite creates its geometry columns with a baked in SRID (a code which usually maps to the EPSG identifier). You can see which SRID was used for a specific geometry using the srid() function:

select srid(geometry) from "nps-spatialite" limit 1

SpatiaLite can also convert to another projection using the Transform() function:

select ':' || AsGeoJSON(Transform(geometry, 2227)) from "nps-spatialite" limit 1

(I'm using ':' || AsGeoJSON(...) here to disable the datasette-leaflet-geojson plugin, since it can't correctly render data that has been transformed to a non-WGS-84 proection.)

Pulling it all together

I now have two tools for imorting geospatial data into SQLite (or SpatiaLite) databases: shapefile-to-sqlite and geojson-to-sqlite.

I'm excited about Datasette's potential as a tool for GIS. I started exploring this back in 2017 when I used it to build a location to timezone API - but adding easy shapefile imports to the toolchain should unlock all kinds of interesting new geospatial projects.

Tags: geospatial, projects, shapefiles, spatialite, sqlite, geojson, weeknotes, tom-macwright, leaflet

geojson-to-sqlite

2020-01-31T06:40:53+00:00

geojson-to-sqlite

I just put out the first release of geojson-to-sqlite—a CLI tool that can convert GeoJSON files (consisting of a Feature or a set of features in a FeatureCollection) into a table in a SQLite database. If you use the --spatialite option it will initalize the table with SpatiaLite and store the geometries in a spacially indexed geometry field—without that option it stores them as GeoJSON.

Via @simonw

Tags: cli, geospatial, projects, spatialite, sqlite, geojson

SpatiaLite — Datasette documentation

2018-05-30T04:34:06+00:00

SpatiaLite — Datasette documentation

Datasette’s documentation now includes extensive coverage of the SpatiaLite extension for SQLite: how to install it, how to import latitude/longitude points, shapefiles and GeoJSON data into SpatiaLite tables, and how to run SQL queries against it that take advantage of spatial indexes. I’m learning SpatiaLite at the moment and filling out the documentation with each new trick I learn as I go—as Mark Pilgrim once taught me, the best way to learn a new technology is to write about it.

Tags: documentation, mark-pilgrim, spatialite, sqlite, datasette

VirtualKNN for SpatiaLite

2018-05-21T21:23:43+00:00

VirtualKNN for SpatiaLite

This looks amazing: a special virtual table shipped as part of SpatiaLite 4.4.0 which implements a fast, R-Tree backed mechanism for finding the X nearest points against a geospatial database table. There’s just one catch: it’s only available in 4.4.0, but the most recent “stable” release of SpatiaLite is 4.3.0a from September 2015 so the version you get if you install from apt-get or homebrew doesn’t yet have this functionality. I’d love to figure out a neat way to package and distribute this along with Datasette. I’d also like to figure out a clean way to ship a more recent version of SQLite than the one that is currently packaged with Python 3 (3.16.2, where the latest SQLite release is 3.23.1).

Tags: geospatial, spatialite, sqlite

GeoCouch: Geospatial queries with CouchDB

2008-10-27T23:48:46+00:00

GeoCouch: Geospatial queries with CouchDB

Interesting approach: uses “external2”, a branch that allows external services to be called from CouchDB. SQLite’s SpatiaLite extension is then used as an external spacial index.

Tags: couchdb, erlang, external2, geocouch, spatialite, sqlite