Calculating a spherical polygon centroid

Question

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 1

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.

Answer 2

It's a good idea to enumerate properties that the centroid of a polygon should have. Here are my criteria:

(a) It's a property of the polygon interior (instead of the vertices or edges). Thus, splitting an edge in two by inserting an additional vertex should not change the position of the centroid. Note that Jenness' definition of the centroid fails on this criterion, since the position of the centroid will depend on how a polygon is divided into triangles.

(b) Perturbing the shape of the polygon by a little should move the centroid by a little. It's necessary here to impose a restriction on the overall extent of the polygon (e.g., to a single hemisphere). Without this restriction, it's easy to construct cases where the centroid will suddenly swing to the opposite side of the earth with a slight movement of a vertex. This condition excludes methods which require that the centroid lie inside the polygon.

(c) It should reduce to the planar definition of centroid for small polygons.

Here are two approaches satisfy these criteria:

(1) Compute the centroid for ellipsoidal polygon in three dimensions and project back to ellipsoid surface (along a normal to the ellipsoid). Big advantage: the centroid can be computed by breaking polygon into simpler shapes.

(2) The centroid is the point with a minimum RMS geodesic distance to all the points in the interior of the polygon. See Buss and Fillmore, "Spherical Averages and Applications to Spherical Splines and Interpolation", ACM Transactions on Graphics 20, 95–126 (2001). Big advantage: the resulting point doesn't depend on how surface is embedded in R³.

Unfortunately, neither of these definitions are straightforward to put into practice. However, the first method can be carried out simply for a sphere. The best "elementary" area to use is the quadrilateral bounded by an edge of the polygon, two meridians through the end-points of the edge, and the equator. The result for the whole polygon entails summing the contributions over the edges. (Additional steps need to be taken if the polygon encircles a pole.)

Suppose the end-points of the edge are (φ₁, λ₁) and (φ₂, λ₂). Let the azimuths of the edge and the endpoints by α₁ and α₂. Assuming the radius of the sphere is 1, the area of the quadrilateral is

  A = α₂ − α₁
      = 2 tan⁻¹ [tan ½(λ₂ − λ₁) sin ½(φ₂ + φ₁) / cos ½(φ₂ + φ₁)]

(This formula for the area, due to Bessel, is substantially better behaved numerically than the commonly used L'Huilier's formula of the area of a triangle.)

The components of the centroid for this quadrilateral are given by

  2 A ⟨x⟩ = φ₂ sin(λ₂ − λ₀) − φ₁ sin(λ₁ − λ₀)
  2 A ⟨y⟩ = cos α₀ (σ₂ − σ₁) − (φ₂ cos(λ₂ − λ₀) − φ₁ cos(λ₁ − λ₀))
  2 A ⟨z⟩ = (λ₂ − λ₁) − sin α₀ (σ₂ − σ₁)

where σ₂ − σ₁ is the length of the edge, and λ₀ and α₀ are the longitude and azimuth of the edge where it crosses the equator, and the x and y axes are oriented so that the equator crossing is at x = 1, y = 0. (z is the axis through the pole, of course.)

Answer 3

For anyone passing by all this time later, there's a good solution posted on https://stackoverflow.com/questions/19897187#answer-38201499 giving the centre as Σ[a×b / |a×b| ⋅ θab/2] for each pair of consecutive vertices a, b – there’s a JavaScript implentation at https://github.com/chrisveness/geodesy/blob/8f4ef33d/latlon-nvector-spherical.js#L783