ST_Area bug?

Mon Mar 31 10:37:47 PDT 2025

Regarding:

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.

and

3. Some tricky transformations

At least I found a way to do that:

SELECT ST_Area(
>     ST_Transform(
>         ST_Segmentize(
>             ST_Transform(
>                 ST_GeometryFromText('POLYGON ((-0.780387531 -79.736894381,
>                 170.375361283 -79.736894381,
>                 170.375361283 -79.692508368,
>                 -0.780387531 -79.692508368,
>                 -0.780387531 -79.736894381))'
>                 , 4326) -- set wgs84 srid to geometry
>             , 3857) -- wgs84 to web-mercator
>         , 1000) -- segmentation max 1 km by rhumb edges
>     , 4326) -- web-mercator to wgs84
> ::geography) / 1e+6 -- in km^2

The result is ~16915 km^2

вс, 30 мар. 2025 г. в 13:22, Alexander Trufanov <trufanovan at gmail.com>:

> 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 <mwtoews at gmail.com>:
>
>> 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 <trufanovan at gmail.com>
>> 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
>

-- 
With best regards,
Alexander Trufanov
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20250331/a9c0175d/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 928443 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20250331/a9c0175d/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 112373 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20250331/a9c0175d/attachment-0001.png>