# 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.

-

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.

-
You've demonstrated that errors can accumulate for large n, but it's not clear that the approach is necessarily flawed. What value of n did you use to achieve the limiting value? – Jason Davies Dec 7 '12 at 21:13
Also, thanks very much for doing the calculations and looking at this in-depth. I just want a little more clarification before I can put the issue to rest. :) – Jason Davies Dec 7 '12 at 21:19
Jason, I added a preliminary example to give you some intuition. The limit itself is rapidly approached; a few dozen lunes will get you several significant digits. But the new example ought to lay to rest any lingering doubts that this weighted triangulation does anything reasonable--except for small triangles, where you're much better off doing the calculations in projected coordinates in the first place. The only reason for doing spherical calculations is when your area of analysis is truly global, for then all projections introduce a lot of distortion. – whuber Dec 7 '12 at 21:24
Fantastic, thanks. So if I understand correctly, simply averaging the Cartesian vectors doesn't result in a reasonable centroid for a spherical triangle (particularly large ones like your first example). I will investigate better methods e.g. finding the intersection of great-circle medians. – Jason Davies Dec 7 '12 at 23:16
BTW, I'm still hopeful that a spherical-area-weighted centroid similar to the above will work. Imagine each polygon being given a 3D volume by adding a vertex at the sphere's origin. Then suspend the sphere by an invisible string connected to its origin and find a stable equilibrium. The centroid is the bottommost point (it's the projection of the centre of mass onto the spherical surface). This should work aside from a few ambiguous cases e.g. a strip going around the equator, where I can just pick a sensible point. Happy to discuss in a new question if you think it's worth it. – Jason Davies Dec 7 '12 at 23:27

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 R3.

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 (φ1, λ1) and (φ2, λ2). Let the azimuths of the edge and the endpoints by α1 and α2. Assuming the radius of the sphere is 1, the area of the quadrilateral is

A = α2 − α1
= 2 tan−1 [tan ½(λ2 − λ1) sin ½(φ2 + φ1) / cos ½(φ2 + φ1)]

(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 Ax⟩ = φ2 sin(λ2 − λ0) − φ1 sin(λ1 − λ0)
2 Ay⟩ = cos α02 − σ1) − (φ2 cos(λ2 − λ0) − φ1 cos(λ1 − λ0))
2 Az⟩ = (λ2 − λ1) − sin α02 − σ1)

where σ2 − σ1 is the length of the edge, and λ0 and α0 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.)

-
Can you explain why the position of Jenness' centroids will depend on how a polygon is divided into triangles? I know from @whuber's example that Jenness' centroid calculation is wrong for spherical triangles, but what if a centroid based on spherical triangle medians is used instead? Will this still fail? – Jason Davies Dec 27 '12 at 18:24
Jenness effectively replaces the spherical polygon by a set of planar triangles and computes their centroid. Clearly(?), the result will depend on the partitioning. Doing the calculation I outlined using the centroids of spherical triangles is fine. You can find the formula for the centroid in J. E. Brock, The Inertia Tensor for a Spherical Triangle, J. Applied Mechanics 42, 239 (1975) dx.doi.org/10.1115/1.3423535 – cffk Dec 28 '12 at 11:44
I took another look at Brock's paper. His formula for the center of mass of a spherical triangle involves a sum over the edges of the triangle. It can therefore be trivially generalized to apply to a polygon (without the need to break it up into triangles). – cffk Dec 28 '12 at 12:14
Great, thanks for the reference. – Jason Davies Dec 28 '12 at 13:59
Do you mind providing a reference for the area calculation due to Bessel, too? I can't seem to find it anywhere, and I'm interested in writing a fast (and accurate) spherical polygon area routine. Thanks! – Jason Davies Jan 10 '13 at 21:28