geometry - Calculating a spherical polygon centroid

I'd like a general way to calculate centroids for polygons on a sphere.

So far, the best online reference appears to be:

Tools for Graphics and Shapes by Jeff Jenness.

The method described there suggests decomposing the polygon into multiple spherical triangles, and computing the average of spherical triangle centroids, weighted by the spherical triangle area.

I know that there are several ways to define a spherical polygon centroid, but I'm looking for something analogous to the following definitions for points and polylines:

• Points: arithmetic mean of the Cartesian vectors representing the points.


• Polylines: weighted mean of Cartesian vectors representing midpoints of each line segment, weighted by the (spherical) length of each segment.

It seems a reasonable continuation to have polygon centroids defined as the weighted mean of the triangular decomposition, weighted by area.

My question is whether the method in the above reference will work regardless of the triangle decomposition used. In particular, it mentions decomposing into triangles relative to an arbitrary point, even external to the polygon, such that some triangles will have negative areas that contribute a negative weight.

Related: How to find the center of geometry of an object?

Answer

It won't work consistently even when you perform all triangulations relative to a single, fixed point. The problem is that spherical and Euclidean calculations are being mixed without any consideration of what they might mean.

One way to make this obvious is to consider a rather extreme triangle, such as almost one-half of a hemisphere. For instance, starting at (lon,lat) = (-179, 0), run along the equator to (0, 0), then up to the north pole at (0, 90), then back to the beginning at (-179, 0). This is a 90-179-90 triangle comprising most of the northern half of the western hemisphere. The problem is that its endpoints (shown as white dots in the figure) lie practically in a plane: one is at the pole and the other two are almost on opposite sides of it. Thus their average, projected back to the sphere (the red dot), is almost at the pole--but that's about as far from any reasonable center as one can get:

As another example, let's triangulate a polygon representing the upper hemisphere relative to its center, the North Pole. We will always divide the Western hemisphere into two equal halves, each of them a 90-90-90 triangle (thereby avoiding any problems with huge, hemisphere-spanning triangles). The Eastern hemisphere, however, will be divided into n equal semi-lunes. The vertices of lune k (k = 1, 2, ..., n) have (lon, lat) coordinates

((k-1) * 180/n, 0),  (k * 180/n, 0),  (k * 180/n, 90).

This figure shows the setup for k=8. The red dots are the individual triangle "centers" computed according to the "Tools for Graphics and Shapes" document, pp 65-67.

Doing the calculations, I find that with k = 2, the area-weighted center indeed is at the North Pole (as would be indicated by symmetry considerations), but as n increases, the result quickly shifts into the Western hemisphere and, in the limit, approaches a latitude of 89.556 degrees along longitude -90 degrees. This is approximately 50 kilometers south of the North Pole itself.

Admittedly, a +/-50 kilometer error for a polygon spanning 20,000 kilometers is small; the total amount of arbitrary variation due to different triangulations in this case is only 0.5%. Evidently the relative errors can be made arbitrarily large by including negative triangles (just add and subtract some really large triangles relative to a small triangle). Regardless, anyone going to the effort of doing spherical computations evidently is trying to avoid projection errors, so they are looking for high accuracy. This triangulation method cannot be recommended.