Viewing GeoPackage data with SpatiaLite and Datasette

I managed to display polygons fom a GeoPackage database file today, using SpatiaLite and Datasette.

GeoPackage

GeoPackage is a standard format for distributing geospatial data as SQLite databases.

I'm still not entirely clear which binary format it uses to store geometries, but I've figured out how to view them by jumping through quite a few hoops.

I started with an example file from here: https://gisdata.mn.gov/dataset/struc-dnr-wld-mgmt-areas-pub

I could browse this in Datasette like so:

datasette struc_dnr_wld_mgmt_areas_pub.gpkg

But this showed the geom columns as pure binary data.

Using SpatiaLite with Datasette

Datasette can load the SpatiaLite extension, if it is available.

Unfortunately it's very easy to end up with a Python 3 version on OS X which entirely disables loading SQLite extensions.

The version you can install from Homebrew does not have this limitation. You can get an install of Datasette that uses that extension-enabled Python using brew install datasette.

I found I had to use the full path to my Homebrew installation to get that to work, like so:

/usr/local/Cellar/datasette/0.61.1/bin/datasette \
  struc_dnr_wld_mgmt_areas_pub.gpkg --load-extension spatialite

To convert the GeoPackage column to a SpatiaLite column, you can use either of these functions:

select GeomFromGPB(geom) from dnr_wildlife_management_areas_public_facilities_lines
-- or
select CastAutomagic(geom) from dnr_wildlife_management_areas_public_facilities_lines

Once it's a SpatiaLite geometry you can convert to GeoJSON like this:

select AsGeoJSON(GeomFromGPB(geom)) from dnr_wildlife_management_areas_public_facilities_lines

This gave me back GeoJSON! But it looked like this (truncated):

{
    "type": "LineString",
    "coordinates": [
        [
            363204.8576250001,
            5155661.741125
        ],
        [
            363227.3801250001,
            5155745.02675
        ],
        [
            363215.0977499997,
            5155760.594499999
        ]
    ]
}

Those don't look like regular latitude/longitude values. They're not - they are in a different projection.

The gpkg_geometry_columns table has a srs_id column, and it showed that for this geometry column the SRID is 26915 - https://epsg.io/26915 says that's "NAD83 / UTM zone 15N" for "North America - between 96°W and 90°W - onshore and offshore".

So we need to convert it to WGS84, also known as SRID 4326, which is preferred by the GeoJSON ecosystem.

After much digging around I figured out this SQL query:

select AsGeoJSON(ST_Transform(SetSRID(GeomFromGPB(geom), 26915), 4326))
from dnr_wildlife_management_areas_public_facilities_lines

The geometry object we get back from GeomFromGPB() doesn't have an SRID associated with it, so we need to specify that it's in 26915 using SetSRID(). Having done that, we can use ST_Transform() to transform it to WGS84.

InitSpatialMetadata() to fix ST_Transform()

The first time I tried this I got this error:

ST_Transform exception - unable to find the origin SRID

After a LOT of frustration, I figured out this was because the database (here the struc_dnr_wld_mgmt_areas_pub.gpkg file) was missing the spatial_ref_sys table it needs to look up different SRIDs.

This isn't a surprise, because it's not a SpatiaLite databsae - it's a regular SQLite database.

You can create those tables in a database by calling this SpatiaLite function:

select InitSpatialMetadata();

Datasette can't call this function because it writes to the database file, so I used sqlite-utils instead - again using a version installed by Homebrew so it could load the extension:

/usr/local/Cellar/sqlite-utils/3.28/bin/sqlite-utils \
  query struc_dnr_wld_mgmt_areas_pub.gpkg \
  'select InitSpatialMetadata()' \
  --load-extension spatialite

This did the job! Having run this, I could open Datasette against the database file (now upgraded to a SpatiaLite database) and view GeoJSON with:

select AsGeoJSON(ST_Transform(SetSRID(GeomFromGPB(geom), 26915), 4326))
from dnr_wildlife_management_areas_public_facilities_lines

Viewing it on a map

Install the datasette-geojson-map by Chris Amico to see that data on a map.

The plugin needs you to return a geometry column, like this:

select ST_Transform(SetSRID(GeomFromGPB(geom), 26915), 4326) as geometry
from dnr_wildlife_management_areas_public_facilities_lines

This plugin works against SpatiaLite geometries as well as GeoJSON, so I dropped the AsGeoJSON() call here.

The binary data here is being rendered using datasette-render-binary.

Created 2022-12-11T19:06:14-08:00, updated 2022-12-11T23:16:54-08:00 · History · Edit

Viewing GeoPackage data with SpatiaLite and Datasette

GeoPackage

Using SpatiaLite with Datasette

InitSpatialMetadata() to fix ST_Transform()

Viewing it on a map

Related