8

I have downloaded some data (as GeoPackages) from this excellent website: https://data.mfe.govt.nz/layer/52757-petroleum-wells/ and https://data.mfe.govt.nz/layer/52675-location-and-extent-of-nzs-aquifers-2015/ with the ultimate goal of writing a third GeoPackage where the two layers are related according to their spatial relationship (i.e. aquifers with references to intersecting wells).

QGIS is able to read these GeoPackages, and to display their geometries.

However when I use the Spatialite GUI (v4.3.0a) to query these datasets, I am running into odd geometry validation problems.

For example:

SELECT Count(*), GeometryType("geom"), Srid("geom"), CoordDimension("geom")
FROM "petroleum_wells"
GROUP BY 2, 3, 4

returns

1188    POINT    2193    NULL

i.e. that CoordDimension("geom") is NULL. This sensibly flows into literally every ST_* function I've tried: they all return -1 for any function operating on the geometries, which is explained in the documentation as: "–1 for UNKNOWN corresponding to a function invocation on NULL arguments." However, the geometry itself is not NULL. It is returned as the outcome of a simple select, albeit with a rather funny tpyoe: BLOB sz=47445 GeoPackegeGEOMETRY.

Oddly, ST_Is3D(geom) correctly returns 0 (false) when it could return -1 (unknown)... so it seems to have at least an understanding that's it's not in 3D.

ST_IsValidReason(geom) claims that the reason is Invalid: NULL Geometry even when I explicitly include a WHERE "geom" IS NOT NULL clause, and have successfully rendered the geometries in QGIS.

Is the GeoPackage geometry really invalid, or is there some issue with Spatialite, or some issue with the way I am using Spatialite?

9

The issue is with the way I am using Spatialite.

After talking with the company behind the web service that provides these datasets (since I initially assumed it may have been an export bug), they revealed to me that there is a distinction between Spatialite geometries and GeoPackage geometries.

To work with GeoPackages in Spatialite, there are several options:

  • Use GeomFromGPB() or AsGPB() everywhere you're working with a GeoPackage geometry.
  • Call AutoGPKStart() then use the vgpkg_ prefixed names.
  • (In Spatialite 4.3+) start queries with SELECT EnableGpkgMode();.

The last option is the easiest, and all spatial queries after that point "just work" as expected.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.