Area function ringArea() yields inconsistent results when circle/polygon contains one of the poles

Hello. I am a PhD student working on a geolocation project that involves computing the intersection area of arbitrarily many spherical circles on the Earth’s surface. While I was using Turf.js to perform these calculations, I noticed something strange. If I took two circles that had the same radius but were located at different latitudes, they would have two different areas. In particular, if one of the circles overlapped the North Pole, its area would be much smaller than a circle of the same size that was closer to the equator.

I attached the Node.js program I used to run these tests. The first command line argument is the latitude in degrees, the second argument is the longitude, and the third argument is the haversine radius in meters. The output is the perceived area in square meters. Here is some sample input/output illustrating the bug for one circle:

% node intersections.js 15 0 300000
283210262181.6476
% node intersections.js 45 0 300000
283210262181.64764
% node intersections.js 60 0 300000
283210262181.6447
% node intersections.js 70 0 300000
283210262181.64355
% node intersections.js 80 0 300000
283210262181.65
% node intersections.js 85 0 300000
283210262181.65076
% node intersections.js 86 0 300000
283210262181.6218
% node intersections.js 89 0 300000
170981367514.91577
% node intersections.js 90 0 300000
0.00903291476967305

% node intersections.js 15 0 3000000
27806712578022.703
% node intersections.js 45 0 3000000
27806712578022.715
% node intersections.js 60 0 3000000
27806712539885.29
% node intersections.js 70 0 3000000
25905348153895.25
% node intersections.js 80 0 3000000
16665046131685.965
% node intersections.js 85 0 3000000
9234108320930.45
% node intersections.js 86 0 3000000
7529485958388.105
% node intersections.js 89 0 3000000
1987895606332.9822
% node intersections.js 90 0 3000000
0.00903291476967305

I looked through the Turf.js GitHub repository and found that the spherical polygon area function ringArea() at https://github.com/Turfjs/turf/blob/master/packages/turf-area/index.ts#L99 seemed to be based on an approximation algorithm on page 7 of the paper "Some algorithms for polygons on a sphere" ( https://trs.jpl.nasa.gov/handle/2014/40409 ). Earlier in the article on pages 3 and 5, it mentioned that the algorithm would only apply if the polygon in question did not contain either of the poles. I believe that this issue can result in significant problems when computing the areas of certain intersection regions, and I couldn’t find it mentioned in the Turf.js documentation.

I understand that computing arbitrary global intersection areas is a nontrivial task and that there are few alternatives to the ringArea() approach. However, I think it might be helpful to mention such a limitation more prominently in the documentation. This might prevent future researchers from obtaining inaccurate results when using the library to compute the areas of polar intersection regions.

If any alternative functions exist in the library to compute areas involving polar circles, I would greatly appreciate any suggestions. Please let me know also if you have any questions or if I should submit a pull request for this issue.

Thank you for your time,

William

Object.defineProperty(exports, "__esModule", { value: true });
//import * as turf from "@turf/turf"
var turf = require("turf");
//import * as fs from 'fs'
var process = require("process"); //For command line arguments
//import * as turfIntersect from '@turf/intersect'
//import * as turfArea from '@turf/area'
var turfCircle = require("@turf/circle");
var commandLineArgs = process.argv.slice(2);
function getArea(circles) {
    if (circles.length == 0) {
        return 0;
    }
    var intersection = circles[0];
    for (var i = 1; i < circles.length; i++) {
        intersection = turf.intersect(intersection, circles[i]);
    }
    return turf.area(intersection);
}
function circlesFromArgs() {
    var circleList = [];
    var index = 0;
    while (index < commandLineArgs.length) {
        var lat = parseFloat(commandLineArgs[index]);
        var long = parseFloat(commandLineArgs[index + 1]);
        var center = [long, lat];
        var radius = parseFloat(commandLineArgs[index + 2]);
        var circle = turfCircle.circle(center, radius, { steps: 128, units: "meters" });
        circleList.push(circle);
        index += 3;
    }
    return circleList;
}
if (commandLineArgs.length % 3 != 0) {
    console.log("%f", 0);
}
else {
    var circles = circlesFromArgs();
    try {
        console.log("%f", getArea(circles));
    }
    catch (_a) {
        console.log("%f", -1);
    }
}

Start with packages/turf-area/index.ts and the ringArea() entry point, then reproduce the polar cases using the attached Node.js program. Compare the behavior with the cited spherical-polygon algorithm and determine whether the polar limitation should be fixed or documented; done means the agreed behavior is covered for polar regions.