[gdal-dev] Spatial view in GeoPackage

Rahkonen Jukka jukka.rahkonen at maanmittauslaitos.fi
Sun May 28 13:00:48 PDT 2023 

Hi,

I found https://github.com/opengeospatial/geopackage/issues/282 from 5 years ago that has also discussion about views and spatial indexes.

-Jukka-

Lähettäjä: Even Rouault <even.rouault at spatialys.com>
Lähetetty: sunnuntai 28. toukokuuta 2023 22.52
Vastaanottaja: Rahkonen Jukka <jukka.rahkonen at maanmittauslaitos.fi>; Bo Victor Thomsen <bo.victor.thomsen at gmail.com>; gdal-dev at lists.osgeo.org
Aihe: Re: [gdal-dev] Spatial view in GeoPackage


Jukka,

currently the driver doesn't try to use the RTree of the base table (it could probably do it, but it would require first to check that both the FID column and geometry column of the same underlying table are available)

Jukka, Bo,

Regarding the use of Spatialite functions, I think the cleanest way is to define an extension to register such computed geometry column. That's what I've done in https://github.com/OSGeo/gdal/pull/7848/commits/8d7468f5342e3e9e54c11c913b8ca53f3effc4b3 . It is the responsibility of the user creating the view to register the extension for its table, as show in https://github.com/OSGeo/gdal/blob/8d7468f5342e3e9e54c11c913b8ca53f3effc4b3/doc/source/drivers/vector/gpkg_spatialite_computed_geom_column.rst#extension-name-or-template

The https://github.com/OSGeo/gdal/blob/master/swig/python/gdal-utils/osgeo_utils/samples/validate_gpkg.py script (once updated by the PR) will use that to avoid erroring out.

Even
Le 28/05/2023 à 19:12, Rahkonen Jukka a écrit :
Hi,

Fixed for this use case, but what if user wants to use  ST_Tranform(...  in the view? Or something that may change the geometry type like ST_Intersection, ST_Buffer etc? In these cases I think that at least the usage of the spatial index of the main table should be disabled. But how could GDAL know when to skip the spatial index and when not? With a special column in gpkg_ogr_contents perhaps?

Somehow I do not feel comfortable. Maybe there should be a warning in
https://gdal.org/drivers/vector/gpkg.html#spatial-views that SpatiaLite functions should be used carefully in the views and such databases should be delivered only for known users, while views with uncomputed geometries should be safe and interoperable.

-Jukka Rahkonen-

Lähettäjä: gdal-dev <gdal-dev-bounces at lists.osgeo.org><mailto:gdal-dev-bounces at lists.osgeo.org> Puolesta Even Rouault
Lähetetty: sunnuntai 28. toukokuuta 2023 0.39
Vastaanottaja: Bo Victor Thomsen <bo.victor.thomsen at gmail.com><mailto:bo.victor.thomsen at gmail.com>; gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
Aihe: Re: [gdal-dev] Spatial view in GeoPackage


Bo Victor,

there was indeed a bug in the OGR GPKG driver preventing you to do what you wanted. Fixed per https://github.com/OSGeo/gdal/pull/7848

The issue was specific to the fact that the geometry column of the view is computed, and not just the simple selection of one of a table.

Note that such GeoPackage will probably only be readable by implementations that install Spatialite (or Spatialite compatible functions). It might also be prudent to wrap the result of ST_Multi() with AsGPB() to force creating a GeoPackage geometry blob instead of a Spatialite one. The OGR GPKG driver is tolerant to receiving both types of geometries, but technically the content should be a GeoPackage one.

Inserting into gpkg_ogr_contents is optional. Its purpose is to accelerate counting the number of features. But for layers < let's say 100 k features, it probably doesn't impact performance much.

Even
Le 27/05/2023 à 22:08, Bo Victor Thomsen a écrit :

Hi list -

I'm working with making spatial views "visible" for QGIS, i.e. so they look like a read-only table from the data source manager.

However, I have one problem:  QGIS doesn't recognize the crs for the view. I've tracked the problem down to that OGR dosn't recognize the crs for the view. Hence this mail...

What I've done:

- creating the view: A union of parcels with the same owner (CVR-number)
CREATE VIEW marker_2019_cvr AS
  SELECT "CVR", SUM(COALESCE(st_area(geom),0.0)) as area_m2, count(*) as antal, st_multi(st_union(geom)) as geom from marker_2019 group by 1;

-- insert or update gpkg_contents
INSERT OR REPLACE INTO gpkg_contents(table_name,data_type,identifier,min_x,min_y,max_x,max_y,srs_id)
  SELECT
    'marker_2019_cvr' AS table_name,
    'features' AS data_type,
    'marker_2019_cvr' AS identifier,
    MIN(ST_MinX(geom)) AS min_x,
    MIN(ST_MinY(geom)) AS min_y,
    MAX(ST_MaxX(geom)) AS max_x,
    MAX(ST_MaxY(geom)) AS max_y,
    25832 AS srs_id -- same srs as the base layer for the view
  FROM marker_2019_cvr;

-- insert or update gpkg_geometry_columns
INSERT OR REPLACE INTO gpkg_geometry_columns ( table_name, column_name , geometry_type_name, srs_id, z, m)
  VALUES ('marker_2019_cvr', 'geom', 'MULTIPOLYGON', 25832, 0, 0);

-- insert or update gpkg_ogr_columns (necessary ?)
INSERT OR REPLACE INTO gpkg_ogr_contents(table_name,feature_count)
  SELECT
    'marker_2019_cvr' AS table_name,
    count(*) AS feature_count
  FROM marker_2019_cvr;

AFAIK, this is the necessary steps to make a spatial view visible in a GeoPackage. And it works; sort of...

The QGIS datasource manager shows the view as a spatial table. I can add it to the QGIS map, but... it can't find the CRS for the layer and defaults to the projects CRS.

Digging a litle bit deeper: Using ogrinfo like this:

ogrinfo -ro -so data_marker_vandloeb_soe.gpkg marker_2019_cvr

shows the layers SRS WKT as (unknown).

Using the same command on the base layer marker_2019 gives the correct WKT.

25832 is the correct EPSG number for the CRS. And the definition exists in the gpkg_spatial_ref_sys table.

What am I doing wrong ? Or is it just "out of scope" for OGR ??

--

Med venlig hilsen / Best regards



Bo Victor Thomsen




_______________________________________________

gdal-dev mailing list

gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/gdal-dev

--

http://www.spatialys.com<http://www.spatialys.com/>

My software is free, but my time generally not.

--

http://www.spatialys.com<http://www.spatialys.com/>

My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20230528/d361db82/attachment.htm>

More information about the gdal-dev mailing list