Re: ST_Area bug?

Hi,

The area is 899078042.2 m2 as verified using:
>
> https://geographiclib.sourceforge.io/cgi-bin/Planimeter?type=polygon&rhumb=geodesic&radius=6378137&flattening=1%2F298.257223563&input=-79.736894381+-0.780387531%0D%0A-79.736894381+170.375361283%0D%0A-79.692508368+170.375361283%0D%0A-79.692508368+-0.780387531%0D%0A-79.736894381+-0.780387531&option=Submit


But there is also an option to switch Edge type to Rhumb line which would
give you the following result:
area (m^2) = 16914914491.4
that's much closer to the expected 16700 km^2.
From Planimeter docs
<https://geographiclib.sourceforge.io/C++/doc/Planimeter.1.html>:

> The algorithm for the area of a rhumb polygon is given in Section 3 of C.
> F. F. Karney, *The area of rhumb polygons*, Stud. Geophys. Geod.
> 68(3--4), 99--120 (2024); DOI https://doi.org/10.1007/s11200-024-0709-z.


The original polygon was drawn in Merkator Web projection, so I expect it
to looks like:

[image: image.png]

instead of

[image: image.png]
which is drawn at
http://www.gcmap.com/mapui?P=79.736894381S+0.780387531W-79.736894381S+170.375361283E-79.692508368S+170.375361283E-79.692508368S+0.780387531W-79.736894381S+0.780387531W&MS=wls&DU=km

Unfortunately, there is no rhumb option in ST_Area or ST_Perimeter.
Also it seems I could workaround this by adding additional points along
latitudes manually. But I can't do this automatically (maybe I didn't find
the proper function?) with ST_Segmentize.
And I couldn't think out any SRID trick to make Postgis process this
polygon in a Merkator-expected way.

If my conclusion is right, then I either need:
1. ST_Segmentatize like method to add points onto rhumb edges in case a
polygon was drawn on the merkator plane.
2. ST_Area and ST_Perimeter with rhumb option
3. Some tricky transformations

Are there any ideas for one of these?


вт, 25 мар. 2025 г. в 23:44, Mike Taves <[email protected]>:

> Note this is a tricky polygon in Antarctica:
>
> http://www.gcmap.com/mapui?P=79.736894381S+0.780387531W-79.736894381S+170.375361283E-79.692508368S+170.375361283E-79.692508368S+0.780387531W-79.736894381S+0.780387531W&MS=wls&DU=km
>
> The area is 899078042.2 m2 as verified using:
>
> https://geographiclib.sourceforge.io/cgi-bin/Planimeter?type=polygon&rhumb=geodesic&radius=6378137&flattening=1%2F298.257223563&input=-79.736894381+-0.780387531%0D%0A-79.736894381+170.375361283%0D%0A-79.692508368+170.375361283%0D%0A-79.692508368+-0.780387531%0D%0A-79.736894381+-0.780387531&option=Submit
>
> PostGIS uses GeoGraphicLib's method
> (https://doi.org/10.1007/s00190-012-0578-z) which has high accuracy.
> The bug is likely with other resources.
>
>
> On Wed, 26 Mar 2025 at 05:58, Alexander Trufanov <[email protected]>
> wrote:
> >
> > Hello,
> >
> > I've a following area calculation:
> >
> > select ST_Area(
> > ST_GeometryFromText('POLYGON ((-0.780387531 -79.736894381,
> > 170.375361283 -79.736894381,
> > 170.375361283 -79.692508368,
> > -0.780387531 -79.692508368,
> > -0.780387531 -79.736894381))', 4326)::geography) /1000000
> >
> > which returns ~900 km^2
> >
> > while ST_Perimeter() for the same geography is ~4600km
> > and If I open ST_AsGeoJson() of this geometry in google earth pro and
> some websites I'm getting ~ 16700 km^2.
> > It seems it's a correct result and ST_Area() result is wrong.
> >
> > Is this a known problem?
> >
> > My versions:
> > PostgreSQL 15.10 (Debian 15.10-astra.se1) on x86_64-pc-linux-gnu,
> compiled by gcc (Astra 12.2.0-14.astra3+b1) 12.2.0, 64-bit
> > POSTGIS="3.3.2 68f4434" [EXTENSION] PGSQL="150"
> GEOS="3.11.1-CAPI-1.17.1" PROJ="9.1.1" GDAL="GDAL 3.6.2, released
> 2023/01/02" LIBXML="2.9.14" LIBJSON="0.16" LIBPROTOBUF="1.4.1" WAGYU="0.5.0
> (Internal)" RASTER
> >
> > --
> > With best regards,
> > Truf
>


-- 
With best regards,
Alexander Trufanov