Problem with calculating area in sphericalTrigonometry

[theirix]

Hello! In latest versions of pygeodesy (approximate after 20.8.15) there is a subtle bug in calculating polygon area using sphericalGeometry. Area value depends on the order of coordinates provided to areaOf. Elliptical geometry is not affected by this bug. Check the script for reproducing the problem. It prints the ratio of area with forward and reverse coordinate order. It should be exactly 1.

Not all polygons are calculated in a wrong way. Polygon generated with delta=1 is bad:

{"type": "Polygon", "coordinates": [[[10, 22], [10, 20], [11, 20], [11, 21], [10, 22]]]}

Polygon generated with delta=0 is good:

{"type": "Polygon", "coordinates": [[[10, 21], [10, 20], [11, 20], [11, 21], [10, 21]]]}

The script

import pygeodesy
import json

print("PyGeodesy {}".format(pygeodesy.version))

lon1,lon2,lat1,lat2 = 10, 11, 20, 21

# If this value is 1, area is calculated in a wrong way.
# If it is zero, area is good.
delta = 1

geojson = dict(type='Polygon', coordinates=[[
        [lon1,lat2+delta],
        [lon1,lat1],
        [lon2,lat1],
        [lon2,lat2],
        [lon1,lat2+delta],
    ]])
print("Geojson for problem polygon: " + json.dumps(geojson))

b = [
        pygeodesy.sphericalTrigonometry.LatLon(lat2+delta,lon1),
        pygeodesy.sphericalTrigonometry.LatLon(lat1,lon1),
        pygeodesy.sphericalTrigonometry.LatLon(lat1,lon2),
        pygeodesy.sphericalTrigonometry.LatLon(lat2,lon2)
        ]
area = pygeodesy.sphericalTrigonometry.areaOf(b)
rarea = pygeodesy.sphericalTrigonometry.areaOf(reversed(b))
print("Spherical area of CW {}, CCW {}, err {}".format(area, rarea, area/rarea))

b = [
        pygeodesy.ellipsoidalKarney.LatLon(lat2+delta,lon1),
        pygeodesy.ellipsoidalKarney.LatLon(lat1,lon1),
        pygeodesy.ellipsoidalKarney.LatLon(lat1,lon2),
        pygeodesy.ellipsoidalKarney.LatLon(lat2,lon2)
        ]
area = pygeodesy.ellipsoidalKarney.areaOf(b)
rarea = pygeodesy.ellipsoidalKarney.areaOf(reversed(b))
print("Ellipsoidal area of CW {}, CCW {}, match {}".format(area, rarea, area/rarea))

Script output for the latest version:

PyGeodesy 22.8.4
Geojson for problem polygon:
{"type": "Polygon", "coordinates": [[[10, 22], [10, 20], [11, 20], [11, 21], [10, 22]]]}
Spherical area of CW 255015585029405.88, CCW 17353168501.86581, err 14695.62086036257
Ellipsoidal area of CW 17304786956.666412, CCW 17304786956.666412, match 1.0

Old versions are not affected by this problem:

PyGeodesy 20.8.15
Geojson for problem polygon:
{"type": "Polygon", "coordinates": [[[10, 22], [10, 20], [11, 20], [11, 21], [10, 22]]]}
Spherical area of CW 17353168501.86581, CCW 17353168501.86581, err 1.0
Ellipsoidal area of CW 17304786956.666412, CCW 17304786956.666412, match 1.0

[mrJean1]

Investigated and resolved: the initial forward azimuth/bearing was set to 0 and should be left undefined.

The upcoming release PyGeodesy 22.8.22 will include fix in function sphericalTrigonometric.areaOf and another test case for all spherical and ellipsoidal LatLons.

[theirix]

Great, thank you for the quick fix!

[mrJean1]

PyGeodesy 22.8.22 is available and includes the fix for function sphericalTrigonometry.areaOf, among several other things.

Please confirm and if this issue #71 is fixed, please close it.

Thank you for reporting the problem and especially for providing a test case to reproduce it. Greatly appreciated!

[theirix]

Tested the new version – works as expected.

Glad to help!