Estoy tratando de encontrar la mejor manera de localizar el centroide de una forma arbitraria que se extiende sobre una esfera unitaria, con la entrada ordenada (en el sentido de las agujas del reloj o en sentido contrario a las agujas del reloj) vértices para el límite 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 (medio hemisferio), generalmente no es posible proyectar simplemente los vértices a un plano y usar métodos planares, como se detalla en Wikipedia (lo siento, no se me permiten más de 2 hipervínculos como recién llegado). Un enfoque ligeramente 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 aquí . En esa misma página, 'Cffk' resaltó este documento que describe un método para calcular el centroide de los triángulos esféricos. He tratado de implementar este método, pero sin éxito, y espero que alguien pueda detectar el problema
He mantenido las definiciones de las variables similares a las del documento para facilitar la comparación. La entrada (datos) es una lista de coordenadas de longitud/latitud, convertidas a coordenadas [x,y,z] por el código. Para cada uno de los triángulos he fijado arbitrariamente un punto como el polo +z, los otros dos vértices están compuestos por un par de puntos vecinos a lo largo del límite del polígono. El código camina a lo largo del límite (comenzando en un punto arbitrario), usando cada segmento del límite del polígono como un lado del triángulo a su vez. Se determina un subcentroide para cada uno de estos triángulos esféricos individuales y se ponderan de acuerdo con el área del triángulo y se suman para calcular el centroide total del polígono. No obtengo ningún error al ejecutar el código, pero los centros totales devueltos están claramente equivocados (he ejecutado algunas formas muy básicas en las que la ubicación del centroide es inequívoca). No he encontrado ningún patrón sensato en la ubicación de los centroides devueltos... así que por el momento no estoy seguro de lo que va mal, ni en las matemáticas ni en el código (aunque, la sospecha es la matemática).
El código de abajo debería funcionar con copiar y pegar como si quisiera probarlo. Si tienes instalado matplotlib y numpy, éste trazará los resultados (ignorará el trazado si no lo haces). Sólo tienes 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 perspicacia.
EDITORIAL: Aquí hay una figura que muestra una proyección de la unidad esfera con un polígono y el centroide resultante que calculo del código. Claramente, el centroide está equivocado ya que el polígono es bastante pequeño y convexo, pero aún así el centroide cae fuera de su perímetro. 
EDITAR: Aquí hay un conjunto de coordenadas muy similar a las de arriba, pero en el formato original [lon,lat] que utilizo normalmente (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.633460EDITORIAL: Un par de ejemplos más... con 4 vértices definiendo un cuadrado perfecto centrado en [1,0,0] obtengo el resultado esperado: 
Sin embargo, 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 lejano 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 (pasar del sentido de las agujas del reloj al sentido contrario o viceversa) el centroide se invierte exactamente.
