Opened 9 months ago

Closed 8 months ago

#5671 closed defect (fixed)

Bug in ST_Area function with use_spheroid=false

Reported by: azakh Owned by: pramsey
Priority: medium Milestone: PostGIS 3.3.7
Component: postgis Version: 3.3.x
Keywords: ST_Area, use_spheroid Cc: azakh

Description (last modified by azakh)

ST_Area fails to calculate the area of the geography polygon with one inner ring, with use_spheroid=false:

SELECT ST_Area(
'POLYGON(
  (
    -111.772 33.307,
    -111.772 33.304,
    -111.783 33.304,
    -111.783 33.307,
    -111.772 33.307
  ),
  (
    -111.781442 33.304825,
    -111.776964 33.3048,
    -111.773047 33.304778,
    -111.773177 33.306408,
    -111.781442 33.304825
  )
)
'::geography, false);

The result is:

ERROR:  lwgeom_area_spher(oid) returned area < 0.0

use_spheroid=true version works and gives 270218.45392724127, but I would expect the on-the-sphere version also to work.

The polygon seems to be correct, though it has almost flat angle at the (-111.776964 33.3048) vertex. Let's call the vertex X.

If X is removed then ST_Area is 270337.25718276866.

If we retain X but make a cyclic shift of the inner ring: AXBCA => XBCAX, then the function succeeds but gives 91737.85992092389 which is quite different from the correct value.

Let's give small amends to the X latitude with AXBCA vertex order:
X => (-111.776964 (33.3048 + amend)).

  Amend    |        ST_Area
-------------------------------
-1.0e-05   |   269904.2853637372
-1.0e-06   |   270298.26808849035
-1.0e-07   |   270628.6930016738
-1.0e-08   |   269541.18027120724
-1.0e-09   |   315257.7589488026
-1.0e-10   |   ERROR:  lwgeom_area_spher(oid) returned area < 0.0
-1.0e-11   |   262120.12981526955
-1.0e-12   |    41142.8857341336
-1.0e-13   |   144916.06439020386
-1.0e-14   |    41142.8857341336
-1.0e-15   |   ERROR:  lwgeom_area_spher(oid) returned area < 0.0
 0.0e000   |   ERROR:  lwgeom_area_spher(oid) returned area < 0.0
 1.0e-15   |   ERROR:  lwgeom_area_spher(oid) returned area < 0.0
 1.0e-14   |   159848.6351422007
 1.5e-14   |   113982.3188750676
 5.0e-14   |   243067.1526449226
 7.5e-14   |    95027.49039558775
 1.0e-13   |    41142.8857341336
 5.0e-13   |   159848.59909125822
 7.5e-13   |   320241.96397413884
 1.0e-12   |    91867.01242231275
 2.5e-12   |   165036.9786292576
 5.0e-12   |   204984.361034764
 7.5e-12   |    35520.7773073109
 1.0e-11   |   116531.9857301203
 2.5e-11   |    64257.524312080175
 5.0e-11   |   253059.59064874944
 1.0e-10   |   131857.65595877985
 5.0e-10   |   261657.88463095468
 1.0e-09   |   265687.2263686325
 1.0e-08   |   270029.18385391496
 1.0e-07   |   270240.35224941676
 1.0e-06   |   270325.2882698695
 1.0e-05   |   270771.5448612002

Monotonous area growth is expected but oscillating behaviour alternating with ERRORs is observed.

Change History (14)

comment:1 by azakh, 9 months ago

Description: modified (diff)

comment:2 by robe, 8 months ago

Milestone: PostGIS 3.3.6

comment:3 by robe, 8 months ago

Confirmed, odd it works with use_sphere=true and not with use_sphere=false. I haven't checked the difference between the two if maybe use_sphere=true is using geographiclib and false is not.

comment:4 by pramsey, 8 months ago

Fails different for me, but still fails. (Conflicting answers between use_sphere true and false.) Will have to get into the code, I thought the area stuff was all the same code line, but maybe not?

comment:5 by pramsey, 8 months ago

Confirmed, the implementation of ptarray_area_sphere is native and ptarray_area_spheroid is delegating to GeographicLib (in proj). I think the easiest fix is to throw away the native implementation and just use GeographicLib with a spherical major/semimajor parameterization, which should be more robust.

comment:6 by Paul Ramsey <pramsey@…>, 8 months ago

In 0d0461b4/git:

Calculate sphere area using the spheroid code, but with a spherical initialization. References #5671

comment:7 by Paul Ramsey <pramsey@…>, 8 months ago

In 7bd482a/git:

Revert "Calculate sphere area using the spheroid code, but with a spherical initialization. References #5671"

This reverts commit 0d0461b431256ee0f949bcfd0064acea18deab84.

comment:8 by Paul Ramsey <pramsey@…>, 8 months ago

In b7990ec/git:

Calculate sphere area using the spheroid code,
but with a spherical initialization. Remove
no longer used code. References #5671

comment:9 by Regina Obe <lr@…>, 8 months ago

In 4e82dc50/git:

Get rid of no longer
used function sphere_signed_area

Triggering errors on some bots that treat warnings as errors

References #5671 for PostGIS 3.5.0

comment:10 by Regina Obe <lr@…>, 8 months ago

In 4057a91/git:

Remove sphere_angle no longer used
References #5671 for PostGIS 3.5.0

comment:11 by Regina Obe <lr@…>, 8 months ago

Resolution: fixed
Status: newclosed

In 0976b90/git:

Use geographiclib also for sphere calcs. Closes #5671 for PostGIS 3.4.3

comment:12 by robe, 8 months ago

Resolution: fixed
Status: closedreopened

comment:13 by robe, 8 months ago

Milestone: PostGIS 3.3.6PostGIS 3.3.7

comment:14 by Regina Obe <lr@…>, 8 months ago

Resolution: fixed
Status: reopenedclosed

In 96145b11/git:

Use geographiclib also for sphere calcs.
Closes #5671 for PostGIS 3.3.7

Note: See TracTickets for help on using tickets.