Trilateración usando 3 puntos de latitud y longitud, y 3 distancias

Question

Quiero averiguar la ubicación de un objetivo desconocido (coordenadas de latitud y longitud). Hay 3 puntos conocidos (pares de coordenadas de latitud y longitud) y para cada punto una distancia en kilómetros a la ubicación del objetivo. ¿Cómo puedo calcular las coordenadas de la ubicación del objetivo?

Por ejemplo, digamos que tengo los siguientes puntos de datos

37.418436,-121.963477   0.265710701754km
37.417243,-121.961889   0.234592423446km
37.418692,-121.960194   0.0548954278262km

Lo que me gustaría es saber cuál es la matemática de una función que toma eso como entrada y devuelve 37.417959,-121.961954 como salida.

Entiendo cómo calcular la distancia entre dos puntos, desde http://www.movable-type.co.uk/scripts/latlong.html Entiendo el principio general de que con tres círculos como estos se obtiene exactamente un punto de superposición. Lo que no entiendo es la matemática necesaria para calcular ese punto con esta entrada.

Answer 1

Después de mirar un poco en Wikipedia y la misma pregunta/respuesta en StackOverflow Me imaginé que iba a intentarlo, y tratar de llenar los vacíos.

En primer lugar, no estoy seguro de dónde has sacado la salida, pero parece que está mal. Tracé los puntos en ArcMap, los almacené en un buffer a las distancias especificadas, ejecuté intersect to en los buffers, y luego capturé el vértice de intersección para obtener las soluciones. Su salida propuesta es el punto en verde. He calculado el valor en el cuadro de llamada, que es de unos 3 metros de lo que ArcMap dio para la solución derivada de la intersección.

La matemática de la página de la wikipedia no está tan mal, sólo hay que convertir las coordenadas geodésicas a las cartesianas ECEF, que se pueden encontrar aquí . los términos a/x +h se pueden sustituir por el radio de la esfera auténtica, si no se utiliza un elipsoide.

Probablemente lo más fácil sea darle un código bien documentado, así que aquí está en python

import math
import numpy

#assuming elevation = 0
earthR = 6371
LatA = 37.418436
LonA = -121.963477
DistA = 0.265710701754
LatB = 37.417243
LonB = -121.961889
DistB = 0.234592423446
LatC = 37.418692
LonC = -121.960194
DistC = 0.0548954278262

#using authalic sphere
#if using an ellipsoid this step is slightly different
#Convert geodetic Lat/Long to ECEF xyz
#   1. Convert Lat/Long to radians
#   2. Convert Lat/Long(radians) to ECEF
xA = earthR *(math.cos(math.radians(LatA)) * math.cos(math.radians(LonA)))
yA = earthR *(math.cos(math.radians(LatA)) * math.sin(math.radians(LonA)))
zA = earthR *(math.sin(math.radians(LatA)))

xB = earthR *(math.cos(math.radians(LatB)) * math.cos(math.radians(LonB)))
yB = earthR *(math.cos(math.radians(LatB)) * math.sin(math.radians(LonB)))
zB = earthR *(math.sin(math.radians(LatB)))

xC = earthR *(math.cos(math.radians(LatC)) * math.cos(math.radians(LonC)))
yC = earthR *(math.cos(math.radians(LatC)) * math.sin(math.radians(LonC)))
zC = earthR *(math.sin(math.radians(LatC)))

P1 = numpy.array([xA, yA, zA])
P2 = numpy.array([xB, yB, zB])
P3 = numpy.array([xC, yC, zC])

#from wikipedia
#transform to get circle 1 at origin
#transform to get circle 2 on x axis
ex = (P2 - P1)/(numpy.linalg.norm(P2 - P1))
i = numpy.dot(ex, P3 - P1)
ey = (P3 - P1 - i*ex)/(numpy.linalg.norm(P3 - P1 - i*ex))
ez = numpy.cross(ex,ey)
d = numpy.linalg.norm(P2 - P1)
j = numpy.dot(ey, P3 - P1)

#from wikipedia
#plug and chug using above values
x = (pow(DistA,2) - pow(DistB,2) + pow(d,2))/(2*d)
y = ((pow(DistA,2) - pow(DistC,2) + pow(i,2) + pow(j,2))/(2*j)) - ((i/j)*x)

# only one case shown here
z = numpy.sqrt(pow(DistA,2) - pow(x,2) - pow(y,2))

#triPt is an array with ECEF x,y,z of trilateration point
triPt = P1 + x*ex + y*ey + z*ez

#convert back to lat/long from ECEF
#convert to degrees
lat = math.degrees(math.asin(triPt[2] / earthR))
lon = math.degrees(math.atan2(triPt[1],triPt[0]))

print lat, lon

Answer 2

No estoy seguro de si estoy siendo ingenuo, pero, ¿si amortiguas cada punto por tamaño, y luego intersectas los tres círculos eso te daría la ubicación correcta?

Puede calcular la intersección utilizando las APIs espaciales. Ejemplos:

• GeoScript
• Suite de Topología Java
• NET Topology Suite
• GEOS

Answer 3

Los siguientes apuntes utilizan la geometría planarítmica (es decir, tendrías que proyectar tus coordenadas en un sistema de coordenadas local adecuado).

Mi razonamiento, con un ejemplo trabajado en Python, es el siguiente:

Tome 2 de los puntos de datos (llámelos a y b ). Llamamos a nuestro punto de destino x . Ya conocemos las distancias ax y bx . Podemos calcular la distancia ab utilizando el teorema de Pitágoras.

>>> import math
>>> a = (1, 4)
>>> b = (3, 6)
>>> dist_ax = 3
>>> dist_bx = 5.385
# Pythagoras's theorem
>>> dist_ab = math.sqrt(abs(a[0]-b[0])**2 + abs(a[1]-b[1])**2)
>>> dist_ab
2.8284271247461903

Ahora, puedes calcular los ángulos de estas líneas:

>>> angle_abx = math.acos((dist_bx * dist_bx + dist_ab * dist_ab - dist_ax * dist_ax)/(2 * dist_bx * dist_ab))
>>> math.degrees(angle_abx)
23.202973815040256
>>> angle_bax = math.acos((dist_ax * dist_ax + dist_ab * dist_ab - dist_bx * dist_bx)/(2 * dist_ax * dist_ab))
>>> math.degrees(angle_bax)
134.9915256259537
>>> angle_axb = math.acos((dist_ax * dist_ax + dist_bx * dist_bx - dist_ab * dist_ab)/(2 * dist_ax * dist_bx))
>>> math.degrees(angle_axb)
21.805500559006095

Sin embargo, ahora que conoces los ángulos, puedes calcular dos posibles ubicaciones para x . Entonces, utilizando el tercer punto c se puede calcular cuál es la ubicación correcta.

Answer 4

Esto podría funcionar. Rápidamente de nuevo en python, podrías poner esto en el cuerpo de una función xN,yN = coordenadas de los puntos, r1 & r2 = valores del radio

dX = x2 - x1
dY = y2 - y1

centroidDistance = math.sqrt(math.pow(e,2) + math.pow(dY,2)) #distance from centroids
distancePL = (math.pow(centroidDistance,2) + (math.pow(r1,2) - math.pow(r2,2))) / (2 * centroidDistance) #distance from point to a line splitting the two centroids

rx1 = x1 + (dX *k)/centroidDistance + (dY/centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))
ry1 = y1 + (dY*k)/centroidDistance - (dX /centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))

rx2 = x1 + (dX *k)/centroidDistance - (dY/centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))
ry2 = y1 + (dY*k)/centroidDistance + (dX /centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))

Los valores rx y ry son los valores de retorno (deberían estar en un array) de los dos puntos de intersección en un círculo, si eso ayuda a aclarar las cosas.

Haga esto para los 2 primeros círculos, y luego otra vez para el primero y el último. Si cualquiera de los resultados de la primera iteración se compara con los resultados de la segunda (dentro de cierta tolerancia, probablemente), entonces tienes el punto de intersección. No es una gran solución, especialmente cuando se empieza a añadir más de los puntos en el proceso, pero es el más simple que puedo ver sin ir a resolver un sistema de ecuaciones.

Answer 5

Puede utilizar la API espacial de Postgis (funciones St_Intersection, St_buffer). Como fmark notó, también debes recordar que Postgis usa algoritmos planares, pero para áreas pequeñas, usar la pryección equidistante no introduce mucho error.