qgis / QGIS

QGIS is a free, open source, cross platform (lin/win/mac) geographical information system (GIS)
https://qgis.org
GNU General Public License v2.0
10.34k stars 2.98k forks source link

ST_Area(geom,use_ellipsoid) not working in virtual layer using a Spatialite datasource on a system with SpatiaLite 5 #41890

Open pigreco opened 3 years ago

pigreco commented 3 years ago

I tried to calculate the area of polygons using a spatialite vector :

Virtual Layers:

image

the result is this:

image

🔗 https://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html

DBManager:

image

with DBManager it calculates all areas

SO and Version

OSGeo4W64-NG:

🔗 http://download.osgeo.org/osgeo4w/testing/osgeo4w-setup.exe

image

image

here, geojson file for testing

https://gist.github.com/pigreco/f0729f630a7298f13f62223a14d47cf1

jef-n commented 3 years ago

Findings:

Trying following with plain spatialite to verify it's not an qgis issue (with utmzones_area.geojson from above gist):

SELECT load_extension('mod_spatialite');
SELECT InitSpatialMetadata(1);
SELECT 
sqlite_version() as sqlite_version,
spatialite_version() as spatialite_version,
proj_version() as proj_version,
geos_version() as geos_version,
rttopo_version() as rttopo_version,
libxml2_version() as libxml2_version,
spatialite_target_cpu() as spatialite_target_cpu;
SELECT ImportGeoJSON('utmzones_area.geojson', 'utmzones_area', 'geom', 1, 4326);
UPDATE utmzones_area
 SET (area_spatialite_5_true,area_spatialite_5_false) =
(
        select st_area(geom, true) as area_spatialite_5_true,
               st_area(geom, false) as area_spatialite_5_false
        from utmzones_area t
        where t.pk = utmzones_area.pk
);

On osgeo4w testing:

D:\TEMP\osgeo4w-2>set SPATIALITE_SECURITY=relaxed
D:\TEMP\osgeo4w-2>sqlite3 <check.sql

1
3.35.2|5.0.1|Rel. 8.0.0, March 1st, 2021|3.9.1-CAPI-1.14.2|1.1.0-dev|2.9.10|Windows_64bit
1204
RTTOPO error: ptarray_area_spheroid: cannot handle ptarray that crosses equator

[62x in total]

And on debian unstable (older sqlite3, older PROJ, slightly different RTTOPO):

$ SPATIALITE_SECURITY=relaxed sqlite3 <check.sql

1
3.34.1|5.0.1|Rel. 7.2.1, January 1st, 2021|3.9.0-CAPI-1.16.2|1.1.0|2.9.10|x86_64-linux-gnu
1204
RTTOPO error: ptarray_area_spheroid: cannot handle ptarray that crosses equator

[62x in total]

so it could be an upstream issue - in rttopo maybe?

Not sure what I make of the GEOS CAPI: 3.9.0 has 1.16.2 on debian and on 3.9.1 in osgeo4w 1.14.2, so the higher GEOS version has a lower CAPI version!?

agiudiceandrea commented 3 years ago

@pigreco could you test Spatialite 5 st_area (geometry, true) in Virtual Layers using a layer that contains only one simple feature which geometry doesn't cross the Equator?

The following Spatialite layer (zipped) contains only one polygon with and EPSG:4326 extent of 12, 40 : 18, 48 not crossing the Equator: test_spatialite5_spatialite.zip (the same layer as GeoJson: test_spatialite5_geojson.zip)

Anyway, why Spatialite 5 st_area (geometry, true) works in DBManager, but it should not work in Virtual Layers?

pigreco commented 3 years ago

@agiudiceandrea image

OSGeo4W64 testing

image

agiudiceandrea commented 3 years ago

@jef-n so it seems Spatialite 5 st_area (geometry, true) doesn't work in Virtual Layers also with a simple layer which contains only one polygon with an EPSG:4326 extent of 12, 40 : 18, 48 not crossing the Equator.

agiudiceandrea commented 3 years ago

Using test_spatialite5,geojson (from test_spatialite5_geojson.zip), the following query (modified from https://github.com/qgis/QGIS/issues/41890#issuecomment-806497317):

SELECT load_extension('mod_spatialite');
SELECT InitSpatialMetadata(1);
SELECT sqlite_version(), spatialite_version(), proj_version(),
       geos_version(), rttopo_version(), libxml2_version(),
       spatialite_target_cpu();
SELECT ImportGeoJSON('test_spatialite5.geojson', 'test_spatialite5', 'geom', 1, 4326);
SELECT st_area(geom), st_area(geom, true), st_area(geom, false)
FROM test_spatialite5

correctly returns, on OSGeo4W testing (Windows 7 64 bit) shell:

F:\QGIS\qgis>set SPATIALITE_SECURITY=relaxed
F:\QGIS\qgis>sqlite3 <check.sql

1
3.33.0|5.0.1-rc1|Rel. 7.2.1, January 1st, 2021|3.9.1-CAPI-1.14.2|1.1.0-dev|2.9.10|Windows_64bit
1
48.0|427409988645.602|426568461564.972

while I can confirm that

agiudiceandrea commented 2 years ago

See also: https://github.com/qgis/QGIS/issues/50188.

Skipper-is commented 1 year ago

Can confirm that this is still an issue in QGIS 3.28.9 LTR, running Spatialite 5.1.0.

PaulvanderKroft commented 10 months ago

Can confirm that this is still an issue in Qgis 3.34.0-2, running Spatialite 5.1.0.

Some clarification on points made elsewhere: from the spatialite-functions that were handled by modules GEOS / LWGEOM up to spatialite 4.3.0, but are handled by RTTOPO in spatialite 5.x, not a single one is functional in a virtual layer in Qgis. Equally: not a single function acting on a geometry, added since the use of RTTOPO (even fairly simple ondes such as AsTWKB), is actually functional in a virtual layer in Qgis. The only RTTOPO-related functions that are working in virtual layers are, as far as I can tell, HasRtTopo( ) and rttopo_version( ).

On statements made earlier, that RTTOPO-functions work in DB Manager: these statements are only valid if those queries are executed on a sqlite (or geopackage) connection. Queries using RTTOPO functions will also fail in DB Manager when executed on a "Virtual Layers / Project Layers" connection.