Locating the centroid (center of mass) of spherical polygons

Question

I'm trying to work out how best to locate the centroid of an arbitrary shape draped over a unit sphere, with the input being ordered (clockwise or anti-cw) vertices for the shape boundary. The density of vertices is irregular along the boundary, so the arc-lengths between them are not generally equal. Because the shapes may be very large (half a hemisphere) it is generally not possible to simply project the vertices to a plane and use planar methods, as detailed on Wikipedia (sorry I'm not allowed more than 2 hyperlinks as a newcomer). A slightly better approach involves the use of planar geometry manipulated in spherical coordinates, but again, with large polygons this method fails, as nicely illustrated here. On that same page, 'Cffk' highlighted this paper which describes a method for calculating the centroid of spherical triangles. I've tried to implement this method, but without success, and I'm hoping someone can spot the problem?

I have kept the variable definitions similar to those in the paper to make it easier to compare. The input (data) is a list of longitude/latitude coordinates, converted to [x,y,z] coordinates by the code. For each of the triangles I have arbitrarily fixed one point to be the +z-pole, the other two vertices being composed of a pair of neighboring points along the polygon boundary. The code steps along the boundary (starting at an arbitrary point), using each boundary segment of the polygon as a triangle side in turn. A sub-centroid is determined for each of these individual spherical triangles and they are weighted according to triangle area and added to calculate the total polygon centroid. I don't get any errors when running the code, but the total centroids returned are clearly wrong (I have run some very basic shapes where the centroid location is unambiguous). I haven't found any sensible pattern in the location of the centroids returned...so at the moment I'm not sure what is going wrong, either in the math or code (although, the suspicion is the math).

The code below should work copy-paste as is if you would like to try it. If you have matplotlib and numpy installed, it will plot the results (it will ignore plotting if you don't). You just have to put the longitude/latitude data below the code into a text file called example.txt.

from math import * try:     import matplotlib as mpl     import matplotlib.pyplot     from mpl_toolkits.mplot3d import Axes3D     import numpy     plotting_enabled = True except ImportError:     plotting_enabled = False  def sph_car(point):     if len(point) == 2:         point.append(1.0)     rlon = radians(float(point[0]))     rlat = radians(float(point[1]))     x = cos(rlat) * cos(rlon) * point[2]     y = cos(rlat) * sin(rlon) * point[2]     z = sin(rlat) * point[2]     return [x, y, z]  def xprod(v1, v2):     x = v1[1] * v2[2] - v1[2] * v2[1]     y = v1[2] * v2[0] - v1[0] * v2[2]     z = v1[0] * v2[1] - v1[1] * v2[0]     return [x, y, z]  def dprod(v1, v2):     dot = 0     for i in range(3):         dot += v1[i] * v2[i]     return dot  def plot(poly_xyz, g_xyz):     fig = mpl.pyplot.figure()     ax = fig.add_subplot(111, projection='3d')     # plot the unit sphere     u = numpy.linspace(0, 2 * numpy.pi, 100)     v = numpy.linspace(-1 * numpy.pi / 2, numpy.pi / 2, 100)     x = numpy.outer(numpy.cos(u), numpy.sin(v))     y = numpy.outer(numpy.sin(u), numpy.sin(v))     z = numpy.outer(numpy.ones(numpy.size(u)), numpy.cos(v))     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', linewidth=0,                     alpha=0.3)     # plot 3d and flattened polygon     x, y, z = zip(*poly_xyz)     ax.plot(x, y, z)     ax.plot(x, y, zs=0)     # plot the alleged 3d and flattened centroid     x, y, z = g_xyz     ax.scatter(x, y, z, c='r')     ax.scatter(x, y, 0, c='r')     # display     ax.set_xlim3d(-1, 1)     ax.set_ylim3d(-1, 1)     ax.set_zlim3d(0, 1)     mpl.pyplot.show()   lons, lats, v = list(), list(), list() # put the two-column data at the bottom of the question into a file called # example.txt in the same directory as this script with open('example.txt') as f:     for line in f.readlines():         sep = line.split()         lons.append(float(sep[0]))         lats.append(float(sep[1])) # convert spherical coordinates to cartesian for lon, lat in zip(lons, lats):     v.append(sph_car([lon, lat, 1.0]))  # z unit vector/pole ('north pole'). This is an arbitrary point selected to act as one #(fixed) vertex of the summed spherical triangles. The other two vertices of any #triangle are composed of neighboring vertices from the polygon boundary. np = [0.0, 0.0, 1.0] # Gx,Gy,Gz are the cartesian coordinates of the calculated centroid Gx, Gy, Gz = 0.0, 0.0, 0.0 for i in range(-1, len(v) - 1):     # cycle through the boundary vertices of the polygon, from 0 to n     if all((v[i][0] != v[i+1][0],             v[i][1] != v[i+1][1],             v[i][2] != v[i+1][2])):         # this just ignores redundant points which are common in my larger input files         # A,B,C are the internal angles in the triangle: 'np-v[i]-v[i+1]-np'         A = asin(sqrt((dprod(np, xprod(v[i], v[i+1])))**2                       / ((1 - (dprod(v[i+1], np))**2) * (1 - (dprod(np, v[i]))**2))))         B = asin(sqrt((dprod(v[i], xprod(v[i+1], np)))**2                       / ((1 - (dprod(np , v[i]))**2) * (1 - (dprod(v[i], v[i+1]))**2))))         C = asin(sqrt((dprod(v[i + 1], xprod(np, v[i])))**2                       / ((1 - (dprod(v[i], v[i+1]))**2) * (1 - (dprod(v[i+1], np))**2))))         # A/B/Cbar are the vertex angles, such that if 'O' is the sphere center, Abar         # is the angle (v[i]-O-v[i+1])          Abar = acos(dprod(v[i], v[i+1]))         Bbar = acos(dprod(v[i+1], np))         Cbar = acos(dprod(np, v[i]))         # e is the 'spherical excess', as defined on wikipedia         e = A + B + C - pi         # mag1/2/3 are the magnitudes of vectors np,v[i] and v[i+1].         mag1 = 1.0         mag2 = float(sqrt(v[i][0]**2 + v[i][1]**2 + v[i][2]**2))         mag3 = float(sqrt(v[i+1][0]**2 + v[i+1][1]**2 + v[i+1][2]**2))         # vec1/2/3 are cross products, defined here to simplify the equation below.         vec1 = xprod(np, v[i])         vec2 = xprod(v[i], v[i+1])         vec3 = xprod(v[i+1], np)         # multiplying vec1/2/3 by e and respective internal angles, according to the          #posted paper         for x in range(3):             vec1[x] *= Cbar / (2 * e * mag1 * mag2                                * sqrt(1 - (dprod(np, v[i])**2)))             vec2[x] *= Abar / (2 * e * mag2 * mag3                                * sqrt(1 - (dprod(v[i], v[i+1])**2)))             vec3[x] *= Bbar / (2 * e * mag3 * mag1                                * sqrt(1 - (dprod(v[i+1], np)**2)))         Gx += vec1[0] + vec2[0] + vec3[0]         Gy += vec1[1] + vec2[1] + vec3[1]         Gz += vec1[2] + vec2[2] + vec3[2]  approx_expected_Gxyz = (0.78, -0.56, 0.27) print('Approximate Expected Gxyz: {0}\n'       '              Actual Gxyz: {1}'       ''.format(approx_expected_Gxyz, (Gx, Gy, Gz))) if plotting_enabled:     plot(v, (Gx, Gy, Gz)) 

Thanks in advance for any suggestions or insight.

EDIT: Here is a figure that shows a projection of the unit sphere with a polygon and the resulting centroid I calculate from the code. Clearly, the centroid is wrong as the polygon is rather small and convex but yet the centroid falls outside its perimeter.

EDIT: Here is a highly-similar set of coordinates to those above, but in the original [lon,lat] format I normally use (which is now converted to [x,y,z] by the updated code).

  -39.366295      -1.633460   -47.282630      -0.740433   -53.912136       0.741380   -59.004217       2.759183   -63.489005       5.426812   -68.566001       8.712068   -71.394853      11.659135   -66.629580      15.362600   -67.632276      16.827507   -66.459524      19.069327   -63.819523      21.446736   -61.672712      23.532143   -57.538431      25.947815   -52.519889      28.691766   -48.606227      30.646295   -45.000447      31.089437   -41.549866      32.139873   -36.605156      32.956277   -32.010080      34.156692   -29.730629      33.756566   -26.158767      33.714080   -25.821513      34.179648   -23.614658      36.173719   -20.896869      36.977645   -17.991994      35.600074   -13.375742      32.581447   -9.554027      28.675497   -7.825604      26.535234   -7.825604      26.535234   -9.094304      23.363132   -9.564002      22.527385   -9.713885      22.217165   -9.948596      20.367878   -10.496531      16.486580   -11.151919      12.666850   -12.350144       8.800367   -15.446347       4.993373   -20.366139       1.132118   -24.784805      -0.927448   -31.532135      -1.910227   -39.366295      -1.633460 

EDIT: A couple more examples...with 4 vertices defining a perfect square centered at [1,0,0] I get the expected result: However, from a non-symmetric triangle I get a centroid that is nowhere close...the centroid actually falls on the far side of the sphere (here projected onto the front side as the antipode): Interestingly, the centroid estimation appears 'stable' in the sense that if I invert the list (go from clockwise to counterclockwise order or vice-versa) the centroid correspondingly inverts exactly.

Answer 1

Anybody finding this, make sure to check Don Hatch's answer which is probably better.

---

I think this will do it. You should be able to reproduce this result by just copy-pasting the code below.

• You will need to have the latitude and longitude data in a file called longitude and latitude.txt. You can copy-paste the original sample data which is included below the code.
• If you have mplotlib it will additionally produce the plot below
• For non-obvious calculations, I included a link that explains what is going on
• In the graph below, the reference vector is very short (r = 1/10) so that the 3d-centroids are easier to see. You can easily remove the scaling to maximize accuracy.
• Note to op: I rewrote almost everything so I'm not sure exactly where the original code was not working. However, at least I think it was not taking into consideration the need to handle clockwise / counterclockwise triangle vertices.

---

Legend:

• (black line) reference vector
• (small red dots) spherical triangle 3d-centroids
• (large red / blue / green dot) 3d-centroid / projected to the surface / projected to the xy plane
• (blue / green lines) the spherical polygon and the projection onto the xy plane

---

from math import * try:     import matplotlib as mpl     import matplotlib.pyplot     from mpl_toolkits.mplot3d import Axes3D     import numpy     plotting_enabled = True except ImportError:     plotting_enabled = False   def main():     # get base polygon data based on unit sphere     r = 1.0     polygon = get_cartesian_polygon_data(r)     point_count = len(polygon)     reference = ok_reference_for_polygon(polygon)     # decompose the polygon into triangles and record each area and 3d centroid     areas, subcentroids = list(), list()     for ia, a in enumerate(polygon):         # build an a-b-c point set         ib = (ia + 1) % point_count         b, c = polygon[ib], reference         if points_are_equivalent(a, b, 0.001):             continue  # skip nearly identical points         # store the area and 3d centroid         areas.append(area_of_spherical_triangle(r, a, b, c))         tx, ty, tz = zip(a, b, c)         subcentroids.append((sum(tx)/3.0,                              sum(ty)/3.0,                              sum(tz)/3.0))     # combine all the centroids, weighted by their areas     total_area = sum(areas)     subxs, subys, subzs = zip(*subcentroids)     _3d_centroid = (sum(a*subx for a, subx in zip(areas, subxs))/total_area,                     sum(a*suby for a, suby in zip(areas, subys))/total_area,                     sum(a*subz for a, subz in zip(areas, subzs))/total_area)     # shift the final centroid to the surface     surface_centroid = scale_v(1.0 / mag(_3d_centroid), _3d_centroid)     plot(polygon, reference, _3d_centroid, surface_centroid, subcentroids)   def get_cartesian_polygon_data(fixed_radius):     cartesians = list()     with open('longitude and latitude.txt') as f:         for line in f.readlines():             spherical_point = [float(v) for v in line.split()]             if len(spherical_point) == 2:                 spherical_point.append(fixed_radius)             cartesians.append(degree_spherical_to_cartesian(spherical_point))     return cartesians   def ok_reference_for_polygon(polygon):     point_count = len(polygon)     # fix the average of all vectors to minimize float skew     polyx, polyy, polyz = zip(*polygon)     # /10 is for visualization. Remove it to maximize accuracy     return (sum(polyx)/(point_count*10.0),             sum(polyy)/(point_count*10.0),             sum(polyz)/(point_count*10.0))   def points_are_equivalent(a, b, vague_tolerance):     # vague tolerance is something like a percentage tolerance (1% = 0.01)     (ax, ay, az), (bx, by, bz) = a, b     return all(((ax-bx)/ax < vague_tolerance,                 (ay-by)/ay < vague_tolerance,                 (az-bz)/az < vague_tolerance))   def degree_spherical_to_cartesian(point):     rad_lon, rad_lat, r = radians(point[0]), radians(point[1]), point[2]     x = r * cos(rad_lat) * cos(rad_lon)     y = r * cos(rad_lat) * sin(rad_lon)     z = r * sin(rad_lat)     return x, y, z   def area_of_spherical_triangle(r, a, b, c):     # points abc     # build an angle set: A(CAB), B(ABC), C(BCA)     # http://math.stackexchange.com/a/66731/25581     A, B, C = surface_points_to_surface_radians(a, b, c)     E = A + B + C - pi  # E is called the spherical excess     area = r**2 * E     # add or subtract area based on clockwise-ness of a-b-c     # http://stackoverflow.com/a/10032657/377366     if clockwise_or_counter(a, b, c) == 'counter':         area *= -1.0     return area   def surface_points_to_surface_radians(a, b, c):     """build an angle set: A(cab), B(abc), C(bca)"""     points = a, b, c     angles = list()     for i, mid in enumerate(points):         start, end = points[(i - 1) % 3], points[(i + 1) % 3]         x_startmid, x_endmid = xprod(start, mid), xprod(end, mid)         ratio = (dprod(x_startmid, x_endmid)                  / ((mag(x_startmid) * mag(x_endmid))))         angles.append(acos(ratio))     return angles   def clockwise_or_counter(a, b, c):     ab = diff_cartesians(b, a)     bc = diff_cartesians(c, b)     x = xprod(ab, bc)     if x < 0:         return 'clockwise'     elif x > 0:         return 'counter'     else:         raise RuntimeError('The reference point is in the polygon.')   def diff_cartesians(positive, negative):     return tuple(p - n for p, n in zip(positive, negative))   def xprod(v1, v2):     x = v1[1] * v2[2] - v1[2] * v2[1]     y = v1[2] * v2[0] - v1[0] * v2[2]     z = v1[0] * v2[1] - v1[1] * v2[0]     return [x, y, z]   def dprod(v1, v2):     dot = 0     for i in range(3):         dot += v1[i] * v2[i]     return dot   def mag(v1):     return sqrt(v1[0]**2 + v1[1]**2 + v1[2]**2)   def scale_v(scalar, v):     return tuple(scalar * vi for vi in v)   def plot(polygon, reference, _3d_centroid, surface_centroid, subcentroids):     fig = mpl.pyplot.figure()     ax = fig.add_subplot(111, projection='3d')     # plot the unit sphere     u = numpy.linspace(0, 2 * numpy.pi, 100)     v = numpy.linspace(-1 * numpy.pi / 2, numpy.pi / 2, 100)     x = numpy.outer(numpy.cos(u), numpy.sin(v))     y = numpy.outer(numpy.sin(u), numpy.sin(v))     z = numpy.outer(numpy.ones(numpy.size(u)), numpy.cos(v))     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', linewidth=0,                     alpha=0.3)     # plot 3d and flattened polygon     x, y, z = zip(*polygon)     ax.plot(x, y, z, c='b')     ax.plot(x, y, zs=0, c='g')     # plot the 3d centroid     x, y, z = _3d_centroid     ax.scatter(x, y, z, c='r', s=20)     # plot the spherical surface centroid and flattened centroid     x, y, z = surface_centroid     ax.scatter(x, y, z, c='b', s=20)     ax.scatter(x, y, 0, c='g', s=20)     # plot the full set of triangular centroids     x, y, z = zip(*subcentroids)     ax.scatter(x, y, z, c='r', s=4)     # plot the reference vector used to findsub centroids     x, y, z = reference     ax.plot((0, x), (0, y), (0, z), c='k')     ax.scatter(x, y, z, c='k', marker='^')     # display     ax.set_xlim3d(-1, 1)     ax.set_ylim3d(-1, 1)     ax.set_zlim3d(0, 1)     mpl.pyplot.show()  # run it in a function so the main code can appear at the top main() 

---

Here is the longitude and latitude data you can paste into longitude and latitude.txt

  -39.366295      -1.633460   -47.282630      -0.740433   -53.912136       0.741380   -59.004217       2.759183   -63.489005       5.426812   -68.566001       8.712068   -71.394853      11.659135   -66.629580      15.362600   -67.632276      16.827507   -66.459524      19.069327   -63.819523      21.446736   -61.672712      23.532143   -57.538431      25.947815   -52.519889      28.691766   -48.606227      30.646295   -45.000447      31.089437   -41.549866      32.139873   -36.605156      32.956277   -32.010080      34.156692   -29.730629      33.756566   -26.158767      33.714080   -25.821513      34.179648   -23.614658      36.173719   -20.896869      36.977645   -17.991994      35.600074   -13.375742      32.581447   -9.554027      28.675497   -7.825604      26.535234   -7.825604      26.535234   -9.094304      23.363132   -9.564002      22.527385   -9.713885      22.217165   -9.948596      20.367878   -10.496531      16.486580   -11.151919      12.666850   -12.350144       8.800367   -15.446347       4.993373   -20.366139       1.132118   -24.784805      -0.927448   -31.532135      -1.910227   -39.366295      -1.633460 

Answer 2

To clarify: the quantity of interest is the projection of the true 3d centroid (i.e. 3d center-of-mass, i.e. 3d center-of-area) onto the unit sphere.

Since all you care about is the direction from the origin to the 3d centroid, you don't need to bother with areas at all; it's easier to just compute the moment (i.e. 3d centroid times area). The moment of the region to the left of a closed path on the unit sphere is half the integral of the leftward unit vector as you walk around the path. This follows from a non-obvious application of Stokes' theorem; see http://www.owlnet.rice.edu/~fjones/chap13.pdf Problem 13-12.

In particular, for a spherical polygon, the moment is the half the sum of (a x b) / ||a x b|| * (angle between a and b) for each pair of consecutive vertices a,b. (That's for the region to the left of the path; negate it for the region to the right of the path.)

(And if you really did want the 3d centroid, just compute the area and divide the moment by it. Comparing areas might also be useful in choosing which of the two regions to call "the polygon".)

Here's some code; it's really simple:

#!/usr/bin/python  import math  def plus(a,b): return [x+y for x,y in zip(a,b)] def minus(a,b): return [x-y for x,y in zip(a,b)] def cross(a,b): return [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]] def dot(a,b): return sum([x*y for x,y in zip(a,b)]) def length(v): return math.sqrt(dot(v,v)) def normalized(v): l = length(v); return [1,0,0] if l==0 else [x/l for x in v] def addVectorTimesScalar(accumulator, vector, scalar):     for i in xrange(len(accumulator)): accumulator[i] += vector[i] * scalar def angleBetweenUnitVectors(a,b):     # http://www.plunk.org/~hatch/rightway.php     if dot(a,b) < 0:         return math.pi - 2*math.asin(length(plus(a,b))/2.)     else:         return 2*math.asin(length(minus(a,b))/2.)  def sphericalPolygonMoment(verts):     moment = [0.,0.,0.]     for i in xrange(len(verts)):         a = verts[i]         b = verts[(i+1)%len(verts)]         addVectorTimesScalar(moment, normalized(cross(a,b)),                                      angleBetweenUnitVectors(a,b) / 2.)     return moment  if __name__ == '__main__':     import sys     def lonlat_degrees_to_xyz(lon_degrees,lat_degrees):         lon = lon_degrees*(math.pi/180)         lat = lat_degrees*(math.pi/180)         coslat = math.cos(lat)         return [coslat*math.cos(lon), coslat*math.sin(lon), math.sin(lat)]      verts = [lonlat_degrees_to_xyz(*[float(v) for v in line.split()])              for line in sys.stdin.readlines()]     #print "verts = "+`verts`      moment = sphericalPolygonMoment(verts)     print "moment = "+`moment`     print "centroid unit direction = "+`normalized(moment)` 

For the example polygon, this gives the answer (unit vector):

[-0.7644875430808217, 0.579935445918147, -0.2814847687566214] 

This is roughly the same as, but more accurate than, the answer computed by @KobeJohn's code, which uses rough tolerances and planar approximations to the sub-centroids:

[0.7628095787179151, -0.5977153368303585, 0.24669398601094406] 

The directions of the two answers are roughly opposite (so I guess KobeJohn's code decided to take the region to the right of the path in this case).