Does geographiclib have a faster code path for pure sphere based geodetic calculations?

[mboeringa]

Recently, PostGIS got a fix to restore the option for sphere based calculations in commands like ST_Area() with geography data type as input (fix implemented in PostGIS 3.3.7+ and 3.4.3+). The command for this is:

ST_Area(geography geog, boolean use_spheroid = true)

with 'use_spheroid' set to false obviously to get a sphere based calculation.

As, for my specific case, sphere based calculations using geography data type are "good enough" (it involves largely generalized global OpenStreetMap data), I was hoping setting 'use_spheroid' to false would result in faster calculations as using a sphere should result in less complicated math compared to a spheroid. However, from a first quick test, this did not appear to result in a performance enhancement. Both using 'use_spheroid = false/true' appeared to show similar timings, with no significant difference (contrary to what the current PostGIS Help suggests).

Now PostGIS depends on geographiclib for the actual calculations as I understood from the developers after starting a discussion on the PostGIS mailing list. The developers suggested the current PostGIS implementation simply sets flattening=0 for the sphere case, and hence does not use a specific different call to geographiclib or a specific class in geographiclib optimized for sphere based calculations.

This might explain the fact I did not see any difference in timings for 'use_spheroid = false/true'.

Now my question is simply: is there an optimized code path in geographiclib for pure sphere based calculations where flattening=0? Or is one supposed to use something like the AlbersEqualArea class in geographiclib in case one desires to calculate approximate areas for global data without the major distortions inherent in e.g. WebMercator? To be honest, I am not entirely sure how PostGIS currently uses geographiclib. I guess it may be using the PolygonAreaT or Geodesic(Exact) classes to perform area and length calculations for geography data type, but IDK.

[cffk]

No. Setting the flattening f = 0 will follow the same code path as f > 0. The considerations for GeographicLib were:

• less code;
• the spherical limit serves as a check to the general case;
• I assumed that the normal use case was an ellipsoid.

If you want to pursue a special purpose code path for the spherical limit, the best formula to use for the spherical excess is Eq. (64) in Algorithms for geodesics. (This is used in the bowels of GeographicLib's area calculation.) The other important thing that GeographicLib's code keeps track of is whether the polygon encircles the pole; see the text preceding Eq. (64).

I should note a pretty good approximation for the area on an ellipsoid: convert the latitudes for the polygon vertices to the authalic latitude and then apply the spherical formula. This is better than using an equal area projection because it will work equally well anywhere on the ellipsoid.

[mosot624]

Not sure if this is the write place to comment but when geod.Inverse(a1,a2,b1,b2) is used on python. Does it use Vincenty's formulae?

[cffk]

The quick answer is no. The specific implementation of the inverse method is described Algorithms for geodesics, Section 4.

The more complete answer is that Vincenty and the GeographicLib python library both use the same underlying formulation given by Helmert. However, Vincenty's method fails to converge in some cases and this is fixed in GeographicLib.

[mboeringa]

No. Setting the flattening f = 0 will follow the same code path as f > 0. The considerations for GeographicLib were:

* less code;

* the spherical limit serves as a check to the general case;

* I assumed that the normal use case was an ellipsoid.

@cffk

Thanks for the detailed response. I understand the considerations, and agree the ellipsoid / spheroid case is the more regular use case. Still, as in my particular workflow, there may be some justification for having an optimized spherical code path as well.

Also thanks for the alternative suggestions. I'd never heard of the authalic latitude before, but from what I now understand from reading up a bit on the internet, it is supposed to give a closer approximation of the spheroidal areas by matching the sphere's total surface area to that of an ellipsoid / spheroid to achieve higher accuracy than a simple spheroid to sphere conversion might, based maybe only by taking into account the diameter at the equator, if I understand it right. Interesting...

[cffk]

Also thanks for the alternative suggestions. I'd never heard of the authalic latitude before, but from what I now understand from reading up a bit on the internet, it is supposed to give a closer approximation of the spheroidal areas by matching the sphere's total surface area to that of an ellipsoid / spheroid to achieve higher accuracy than a simple spheroid to sphere conversion might, based maybe only by taking into account the diameter at the equator, if I understand it right. Interesting...

This procedure gives the exact ellipsoidal area for the polygon whose edges map into great circles on the authalic sphere. (Such edges will be close to geodesics for sufficiently short edges.)

[mboeringa]

If you want to pursue a special purpose code path for the spherical limit, the best formula to use for the spherical excess is Eq. (64) in Algorithms for geodesics. (This is used in the bowels of GeographicLib's area calculation.)

@cffk

One last question: are you saying with the above statement that, although geographiclib does not have special treatment for the flattening=0 case, it internally already uses an approximation of the spheroidal area based based on some sphere based implementation for area calculations in the e.g. PolygonAreaT class to speed up calculations? This all hits a bit on my limits of understanding of these type of geometric property calculations in geodesy...

[cffk]

If you want to pursue a special purpose code path for the spherical limit, the best formula to use for the spherical excess is Eq. (64) in Algorithms for geodesics. (This is used in the bowels of GeographicLib's area calculation.)

@cffk

One last question: are you saying with the above statement that, although geographiclib does not have special treatment for the flattening=0 case, it internally already uses an approximation of the spheroidal area based based on some sphere based implementation for area calculations in the e.g. PolygonAreaT class to speed up calculations? This all hits a bit on my limits of understanding of these type of geometric property calculations in geodesy...

No, this isn't the way to think about it. The expression for the area, Eq. (59) of my paper, contains two terms on the right-hand side. The first gives the spherical contribution and the second, proportional to e^2, is the ellipsoidal correction. No approximations are involved at this stage. (However, the ellipsoidal term is then expressed as a Taylor series, Eq. (62), and enough terms are retained in this series to give full double-precision accuracy.)

[cffk]

Closing for now...