python - irregular - Localización del centroide(centro de masa) de polígonos esféricos

calcular centroide poligono qgis

Question

Estoy tratando de encontrar la mejor manera de ubicar el centroide de una forma arbitraria desplegada sobre una esfera unitaria, con la entrada ordenada (en sentido horario o anti-cw) vértices para el contorno de la forma. La densidad de los vértices es irregular a lo largo del límite, por lo que las longitudes de arco entre ellos no son generalmente iguales. Debido a que las formas pueden ser muy grandes (la mitad de un hemisferio) generalmente no es posible proyectar los vértices a un plano y usar métodos planos, como se detalla en Wikipedia (lo siento, no tengo más de 2 hipervínculos como recién llegado). Un enfoque un poco mejor implica el uso de geometría plana manipulada en coordenadas esféricas, pero de nuevo, con grandes polígonos este método falla, como se ilustra muy bien here . En esa misma página, ''Cffk'' destacó este artículo que describe un método para calcular el centroide de triángulos esféricos. Intenté implementar este método, pero sin éxito, y espero que alguien pueda detectar el problema.

He mantenido las definiciones de variables similares a las del documento para que sea más fácil compararlas. La entrada (datos) es una lista de coordenadas de longitud / latitud, convertidas en coordenadas [x, y, z] por el código. Para cada uno de los triángulos, he fijado arbitrariamente un punto para que sea el polo + z, y los otros dos vértices se componen de un par de puntos vecinos a lo largo del límite del polígono. El código avanza a lo largo del límite (comenzando en un punto arbitrario), usando cada segmento de límite del polígono como un lado del triángulo a su vez. Se determina un sub-centroide para cada uno de estos triángulos esféricos individuales y se ponderan de acuerdo con el área del triángulo y se agregan para calcular el centroide del polígono total. No obtengo ningún error al ejecutar el código, pero los centroides totales devueltos son claramente incorrectos (he ejecutado algunas formas muy básicas donde la ubicación del centroide no es ambigua). No he encontrado ningún patrón sensible en la ubicación de los centroides devueltos ... por lo que en este momento no estoy seguro de qué es lo que está mal, ya sea en el código matemático (aunque, la sospecha es la matemática).

El código a continuación debería funcionar copiar y pegar como está si le gustaría probarlo. Si tiene matplotlib y numpy instalados, trazará los resultados (ignorará el trazado si no lo hace). Solo tiene que poner los datos de longitud / latitud debajo del código en un archivo de texto llamado 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))

Gracias de antemano por cualquier sugerencia o idea.

EDITAR: Aquí hay una figura que muestra una proyección de la esfera de la unidad con un polígono y el centroide resultante que calculo a partir del código. Claramente, el centroide está equivocado ya que el polígono es bastante pequeño y convexo, pero el centroide cae fuera de su perímetro.

EDITAR: Aquí hay un conjunto de coordenadas muy similar a las anteriores, pero en el formato original [lon, lat] que normalmente uso (que ahora se convierte en [x, y, z] por el código actualizado).

  -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

EDITAR: Un par de ejemplos más ... con 4 vértices que definen un cuadrado perfecto centrado en [1,0,0] obtengo el resultado esperado: Sin embargo, a partir de un triángulo no simétrico obtengo un centroide que no está en ninguna parte cerca ... el centroide en realidad cae en el lado más alejado de la esfera (aquí proyectado en el lado frontal como la antípoda): Curiosamente, la estimación del centroide parece "estable" en el sentido de que si invierto la lista (vaya del orden de las agujas del reloj al sentido antihorario o viceversa), el centroide se invierte correspondientemente.

Answer 1

Creo que esto lo hará. Debería poder reproducir este resultado copiando y pegando el siguiente código.

• Tendrá que tener los datos de latitud y longitud en un archivo llamado longitude and latitude.txt . Puede copiar y pegar los datos de muestra originales que se incluyen debajo del código.
• Si tiene mplotlib, también producirá la gráfica a continuación
• Para cálculos no obvios, incluí un enlace que explica lo que está sucediendo
• En el siguiente gráfico, el vector de referencia es muy corto (r = 1/10), por lo que los centroides 3d son más fáciles de ver. Puede quitar fácilmente la escala para maximizar la precisión.
• Nota para op: reescribí casi todo, así que no estoy seguro de dónde estaba funcionando el código original. Sin embargo, al menos creo que no estaba teniendo en cuenta la necesidad de manejar vértices triangulares en sentido horario / antihorario.

Leyenda:

• (línea negra) vector de referencia
• (pequeños puntos rojos) triángulo esférico 3d-centroides
• (gran punto rojo / azul / verde) 3d-centroide / proyectado a la superficie / proyectado al plano xy
• (líneas azules / verdes) el polígono esférico y la proyección en el plano xy

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()

Aquí están los datos de longitud y latitud que puede pegar en 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

Creo que una buena aproximación sería calcular el centro de masa usando coordenadas cartesianas ponderadas y proyectando el resultado en la esfera (suponiendo que el origen de las coordenadas es (0, 0, 0)^T ).

Let be (p[0], p[1], ... p[n-1]) los n puntos del polígono. El centroide aproximativo (cartesiano) se puede calcular de la siguiente manera:

c = 1 / w * (sum of w[i] * p[i])

mientras que w es la suma de todos los pesos y considerando que p[i] es un punto de polígono y w[i] es un peso para ese punto, por ejemplo

w[i] = |p[i] - p[(i - 1 + n) % n]| / 2 + |p[i] - p[(i + 1) % n]| / 2

mientras que |x| es la longitud de un vector x . Es decir, un punto se pondera con la mitad de la longitud que el anterior y la mitad de la longitud hasta el siguiente punto del polígono.

Este centroide c ahora puede proyectarse en la esfera de la siguiente manera:

c'' = r * c / |c| 

mientras que r es el radio de la esfera.

Para considerar la orientación del polígono (ccw, cw) el resultado puede ser

c'' = - r * c / |c|. 

Answer 3

Lo siento (como usuario recién registrado) tuve que escribir una nueva publicación en lugar de solo votar / comentar la respuesta anterior de Don Hatch. La respuesta de Don, creo, es la mejor y más elegante. Es matemáticamente riguroso al calcular el centro de masa (primer momento de masa) de una manera sencilla cuando se aplica al polígono esférico.

La respuesta de Kobe John es una buena aproximación pero solo satisfactoria para áreas más pequeñas. También noté algunas fallas en el código. En primer lugar, el punto de referencia debe proyectarse a la superficie esférica para calcular el área esférica real. En segundo lugar, la función points_are_equivalent () podría necesitar ser refinada para evitar dividir por cero.

El error de aproximación en el método de Kobe radica en el cálculo del centroide de triángulos esféricos. El sub-centroide NO es el centro de masa del triángulo esférico sino el plano. Esto no es un problema si se trata de determinar ese único triángulo (el signo puede voltearse, ver más abajo). Tampoco es un problema si los triángulos son pequeños (por ejemplo, una triangulación densa del polígono).

Algunas pruebas simples podrían ilustrar el error de aproximación. Por ejemplo, si usamos solo cuatro puntos:

10 -20

10 20

-10 20

-10 -20

La respuesta exacta es (1,0,0) y ambos métodos son buenos. Pero si agrega algunos puntos más a lo largo de un borde (por ejemplo, agregue {10, -15}, {10, -10} ... al primer borde), verá que los resultados del método de Kobe comienzan a cambiar. Además, si aumenta la longitud de [10, -10] a [100, -100], verá que el resultado de Kobe cambia la dirección. Una posible mejora podría ser agregar otro (s) nivel (es) para el cálculo de los sub-centroides (básicamente refinar / reducir tamaños de triángulos).

Para nuestra aplicación, el límite del área esférica está compuesto por arcos múltiples y, por lo tanto, no es un polígono (es decir, el arco no forma parte del gran círculo). Pero esto solo requerirá un poco más de trabajo para encontrar el n-vector en la integración de la curva.

EDITAR: Reemplazar el cálculo del subcentroide con el dado en el documento de Brock debería corregir el método de Kobe. Pero no lo intenté.

Answer 4

Para aclarar: la cantidad de interés es la proyección del centroide 3d verdadero (es decir, el centro de masa 3d, es decir, el centro-de-área 3d) en la esfera de la unidad.

Como lo único que le importa es la dirección desde el origen hasta el centroide 3d, no necesita preocuparse por las áreas en absoluto; es más fácil calcular el momento (es decir, el área de los tiempos del centroide 3D). El momento de la región a la izquierda de un camino cerrado en la esfera de la unidad es la mitad de la integral del vector unitario hacia la izquierda mientras camina alrededor del camino. Esto se desprende de una aplicación no obvia del teorema de Stokes; ver http://www.owlnet.rice.edu/~fjones/chap13.pdf Problema 13-12.

En particular, para un polígono esférico, el momento es la mitad de la suma de (axb) / || axb || * (ángulo entre a y b) para cada par de vértices consecutivos a, b. (Eso es para la región a la izquierda de la ruta; anularla para la región a la derecha de la ruta).

(Y si realmente deseaba el centroide 3d, simplemente calcule el área y divida el momento. Comparando áreas también podría ser útil elegir a cuál de las dos regiones llamar "el polígono").

Aquí hay un código; es realmente 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)`

Para el polígono de ejemplo, esto da la respuesta (vector unitario):

[-0.7644875430808217, 0.579935445918147, -0.2814847687566214]

Esto es más o menos el mismo, pero más preciso que, la respuesta calculada por el código de @ KobeJohn, que usa tolerancias aproximadas y aproximaciones planas a los sub-centroides:

[0.7628095787179151, -0.5977153368303585, 0.24669398601094406]

Las direcciones de las dos respuestas son más bien opuestas (así que supongo que el código de KobeJohn decidió llevar la región a la derecha del camino en este caso).