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 A ⟨x⟩ =
φ2 sin(λ2 − λ0)
−
φ1 sin(λ1 − λ0)
2 A ⟨y⟩ =
cos α0 (σ2 − σ1) − (φ2 cos(λ2 − λ0) − φ1 cos(λ1 − λ0))
2 A ⟨z⟩ =
(λ2 − λ1) − sin α0 (σ2 − σ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.)