Discrepency between Spatialite RTTOPO based ST_Area(geom, use_ellipsoid) and PostGIS/PROJ/QGIS/now GDAL

[Even Rouault]

Hi,

while implementing https://github.com/OSGeo/gdal/pull/9552 and comparing
to other implementations, I found that Spatialite ST_Area(geom,
use_ellipsoid) version doesn't return the same results as the new GDAL
implementation or the PostGIS one (or QGIS since
https://github.com/qgis/QGIS/pull/41726), all relying on functionnality
provided by the geodesic.h header provided by PROJ (actually an extract
of GeographicLib)

Spatialite:

sqlite> select st_area(setsrid(st_geomfromtext('POLYGON((2 49,3 49,3
48,2 49))'), 4326), true);
4068862927.84766

PostGIS:

autotest=# select st_area(st_setsrid('POLYGON((2 49,3 49,3 48,2
49))'::geography, 4326), true);
 4068384291.8911743

Digging further, I see that librttopo has 2 code paths for the geodesic
area computation in rtspheroid.c, one that uses PROJ's extract of
GeographicLib, and another legacy one. This is similar to PostGIS, since
that file is mostly a copy&paste of PostGIS'
https://github.com/postgis/postgis/blob/master/liblwgeom/lwspheroid.c#L32

The root issue is that librttopo's configure.ac has no provision to
build against PROJ, so the legacy code path is used. So this is an
upstream librttopo issue actually. Just filed at
https://git.osgeo.org/gitea/rttopo/librttopo/issues/44

But just wanted to make sure this is known by Spatialite.

A possibility would be that Spatialite directly implements this
functionality by including <geodesic.h> from PROJ and do something like

geod_geodesic g;
geod_init(&g, dfSemiMajor, dfFlattening);
double signed_area = 0;
double lat_array_in_degree[] = ...
double lon_array_in_degree[] = ...
geod_polygonarea(&g, lat_array_in_degree, lon_array_in_degree,
number_of_points, &signed_area, NULL);

Even

--
http://www.spatialys.com
My software is free, but my time generally not.

[sandro furieri]

Hi Even,

Thanks for reporting this odd anomaly.

SpatiaLite already intensively uses the support of PROJ's <geodesic.h> for all Length 

metric measures,  and I sincerely believed that the same was for Area metric measures too. 

Unfortunately when I went to check the code I discovered that it wasn't like that at all,

and that the metric measures of Area were instead delegated to RtTopo.

A rather extravagant and not at all necessary choice that after so long

I don't really know how to explain.

Conclusion: now that the problem has come to light I've modified the code

as it should have been from the beginning.

Starting from my last commit into the Fossil repository the geodesic

measure of Areas now is directly handled by PROJ ignoring RtTopo,

as you also suggested and which is the cleanest choice.

bye Sandro