Geospatial SQL queries in SQLite using TG, sqlite-tg and datasette-sqlite-tg#

Installation#

Since Alex released it as a Datasette plugin, you can run this:

1
datasette install datasette-sqlite-tg

Confirm installation with:

1
datasette plugins

This should just work. In my case I ran into a problem because I’d previously installed an earlier alpha and needed to upgrade the underlying sqlite-tg dependency, which I did like this:

1
datasette install sqlite-tg==0.0.1a6

This won’t be necessary for you if you are installing it for the first time.

Trying it out#

sqlite-tg can work with three formats: WKT, WKB and GeoJSON.

These queries can be used to convert one to the other. Here’s how to convert a WKT bounding box around San Francisco to GeoJSON:

1
select tg_to_geojson('POLYGON((
2
  -122.51610563264538 37.81424532146113,
3
  -122.51610563264538 37.69618409220847,
4
  -122.35290547288255 37.69618409220847,
5
  -122.35290547288255 37.81424532146113,
6
  -122.51610563264538 37.81424532146113
7
))')

This outputs (after pretty-printing):

1
{
2
  "type": "Polygon",
3
  "coordinates": [
4
    [
5
      [-122.51610563264538, 37.81424532146113],
6
      [-122.51610563264538, 37.69618409220847],
7
      [-122.35290547288255, 37.69618409220847],
8
      [-122.35290547288255, 37.81424532146113],
9
      [-122.51610563264538, 37.81424532146113]
10
    ]
11
  ]
12
}

This is already useful: an easily installed mechanism for converting between WKT and GeoJSON with a SQL query is a great thing to have.

Let’s convert a point within San Francisco to WKB binary format:

1
select hex(tg_to_wkb('POINT(-122.4075 37.787994)'))

Outputs:

1
0101000000AE47E17A149A5EC0DCB8C5FCDCE44240

We can convert that back to GeoJSON like this:

1
select tg_to_geojson(x'0101000000AE47E17A149A5EC0DCB8C5FCDCE44240')

This is using SQLite’s x'...' syntax to treat a hexadecimal string as a blob literal.

And here’s how to confirm that our point exists within our polygon bounding box for San Francisco:

1
select tg_intersects('{
2
  "type": "Polygon",
3
  "coordinates": [
4
    [
5
      [-122.51610563264538, 37.81424532146113],
6
      [-122.51610563264538, 37.69618409220847],
7
      [-122.35290547288255, 37.69618409220847],
8
      [-122.35290547288255, 37.81424532146113],
9
      [-122.51610563264538, 37.81424532146113]
10
    ]
11
  ]
12
}', 'POINT(-122.4075 37.787994)')

I’m mixing GeoJSON and WKT here. I get back:

1
1

Because the point and the polygon intersect.

Instead of using the WKT POINT(...) string here we can use the tg_point(longitude, latitude) function. We’ll use that in the next example.

Let’s try with Times Square in New York:

1
select tg_intersects('{
2
  "type": "Polygon",
3
  "coordinates": [
4
    [
5
      [-122.51610563264538, 37.81424532146113],
6
      [-122.51610563264538, 37.69618409220847],
7
      [-122.35290547288255, 37.69618409220847],
8
      [-122.35290547288255, 37.81424532146113],
9
      [-122.51610563264538, 37.81424532146113]
10
    ]
11
  ]
12
}', tg_point(-73.985130, 40.758896))

And we get back:

1
0

Let’s try something a bit more challenging.

Timezone lookups#

I’ve worked with timezone lookups using SpatiaLite and Datasette before - I even wrote a tutorial about it.

I decided to try implementing a version of that on top of sqlite-tg.

I grabbed the latest release of timezones.geojson.zip from evansiroky/timezone-boundary-builder/ and unzipped it to get a combined.geojson file that starts like this:

1
{
2
  "type": "FeatureCollection",
3
  "features": [
4
    {
5
      "type": "Feature",
6
      "properties": {
7
        "tzid": "Africa/Abidjan"
8
      },
9
      "geometry": {
10
        "type": "Polygon",
11
        "coordinates": [
12
          [
13
            [
14
              -5.440683,
15
              4.896553
16
            ],
17
            [
18
              -5.303699,
19
              4.912035
20
            ],

Then I loaded that into a SQLite database table using this sqlite-utils mechanism (explained in What’s new in sqlite-utils 3.20 and 3.21: —lines, —text, —convert).

1
sqlite-utils insert timezones.db timezones combined.json \
2
  --text --convert '
3
import json
4

5
def convert(s):
6
    data = json.loads(s)
7
    for feature in data["features"]:
8
        tzid = feature["properties"]["tzid"]
9
        yield {"tzid": tzid, "geometry": feature["geometry"]}
10
' --pk tzid

This Bash one-liner loads the entire combined.json file (140MB of it) into memory and then runs a snippet of Python to loop through all of those features and yield {"tzid": "...", "geometry": {...}} dictionaries for each timezone.

sqlite-utils then creates a table with this schema and inserts those records as JSON strings:

1
CREATE TABLE [timezones] (
2
   [tzid] TEXT PRIMARY KEY,
3
   [geometry] TEXT
4
);

The resulting database file is 160MB. You can download a copy from here.

Here’s a query that can be used to find the timezone for a latitude/longitude pair:

1
select
2
  tzid
3
from
4
  timezones
5
where
6
  tg_intersects(
7
    'POINT (' || :longitude || ' ' || :latitude || ')',
8
    geometry
9
  )

I started Datasette like this:

1
datasette timezones.db

And ran that query using the http://localhost:8001/timezones page.

And it worked! It returned the correct timezone for the different points I tried.

The query takes around 700ms on my M2 MacBook Pro. That’s not a great number - usable, but pretty taxing - but it’s also not a huge surprise since this is the most naive implementation of this possible - it’s doing a brute force geometry check against 160MB of JSON strings in the table.

sqlite-tg doesn’t yet include support for indexing. TG provides some very interesting indexing mechanisms as part of the library.

And for added fun… I installed the datasette-geojson-map plugin and added select * to the query and I got this!