[Bug]: Polygon areas around a pole: 0 km² for a cap written as a box, the rest of the Earth for a loop around the North Pole #1353

Before you start

• I searched the existing issues and did not find a duplicate.
• This is a bug report, not a usage question (questions belong in Discussions).

What happened?

A polygon around a pole gets a different area, and in some tools a wrong one, depending on which tool measures it. Measured on main (8ebc4ce) for a cap above 80°N, whose true area on the WGS 84 ellipsoid (the zone formula) is 3.909 million km²:

This corrects my #1033 comment, which put the 0 km² on the attribute table's Calculate Geometry and called a cap with vertices along the parallel right: the 0 km² is the Calculate Geometry Attributes tool's, and the attribute table is 0.9% low either way.

Tabulate Intersection's area_km2 uses the same routine as the Calculate Geometry Attributes tool, on each clipped part. Spatial Join's Largest overlap uses a third, the sphere's geodesic_ring_area_m2.

The causes:

• 0 km²: ground_area_perimeter (crates/analysis/src/tools.rs) unwraps each step of longitude the short way, which makes the -180°→180° edge 0° long and the ring a line.
• The sphere's area, 0.9% low: wgs84_ring_perimeter_area_m2 (crates/geom/src/lib.rs) measures any ring with an edge spanning 360° on the sphere instead (added in Measure on the WGS 84 ellipsoid in the attribute table, measure tool and network tools (#939) #998, for [Bug]: The measure tool, the attribute table's Calculate Geometry and network edge lengths are still on the sphere, and disagree with the ellipsoidal geoprocessing tools #939). Every cap closed along the pole has one, the edge along the pole from 180° back to -180°. The sphere routine, geodesic_ring_area_m2, takes geodetic latitudes on a sphere of the mean radius, which is what makes it 0.9% low there.
• The rest of the Earth: geodesic_ring_area_m2 integrates from the South Pole, so a ring winding a full turn of longitude gets the South-Pole side whichever way it runs, and wgs84_ring_perimeter_area_m2's hemisphere test inherits it. The same loops at 80°S give the cap. (I noted this in my review of #998.)
• Perimeters miss an edge along a parallel written as one edge: a geodesic from (-180°, 80°N) to (180°, 80°N) has no length, where the parallel is 6,982 km.

Steps to reproduce

1. Save boxcap.geojson: {"type":"FeatureCollection","features":[{"type":"Feature","properties":{"name":"box"},"geometry":{"type":"Polygon","coordinates":[[[-180,80],[180,80],[180,90],[-180,90],[-180,80]]]}}]}
2. Run both measures and export the result:

   rustgis-cli run --data boxcap.geojson \
     --cmd 'gp.run={"tool":"calculateGeometry","params":{"input":"boxcap"}}' \
     --cmd 'edit.calculateGeometry={"layer":"boxcap","field":"table_km2","create":true,"overwrite":true,"property":"area","unit":"km2","scope":"all"}' \
     --cmd 'share.exportLayer={"layer":"boxcap","path":"out_box.geojson"}'
   

   out_box.geojson has area_km2 0 (the Calculate Geometry Attributes tool) and table_km2 3,874,523 (the attribute table's Calculate Geometry).
3. The same with the 12-point loop [[0,80],[30,80],[60,80],[90,80],[120,80],[150,80],[180,80],[210,80],[240,80],[270,80],[300,80],[330,80],[0,80]]: area_km2 3,736,196 and table_km2 506,329,425.

Expected behavior

#1033 settled the form for polygons around a pole: cut at ±180° and closed along the pole ("Pole-enclosing form: yes, keep it as the canonical form"). A cap written that way should measure the same in every tool: 3.909 M km² with its 80°N edge read as the parallel, and within 0.001 M km² of that with a vertex every degree, whichever line its edges follow. Two fixes need no wider decision, and the third only the choice of side for a loop around a pole (#1354):

• a 360° edge isn't unwrapped to nothing;
• an edge along the pole doesn't send the ring to the sphere;
• a ring winding a full turn gets the same side at either pole.

A way to make the non-geodesic measures exact: add the authalic latitude ξ and radius R_q to rustgis_geom. proj4rs has them only inside the crate. Then any ring whose edges are paths in longitude and latitude has the exact area R_q² ∮ sin ξ dλ: in closed form along a parallel, and by quadrature along an edge straight in longitude and latitude or on the map. Geodesic edges keep GeographicLib's PolygonArea. Which line an edge follows is a wider question (issue to follow). Tabulate Intersection's shares must still add up, so they need measuring along the lines the overlay cuts along.

How are you running RustGIS?

Command line (rustgis-cli)

Version

main at 8ebc4ce, built from source

Operating system, browser and graphics

macOS (Apple silicon)

Sample data

The files in the steps above. I can take the three fixes and the authalic helpers after #1347–#1349 land.