¿Cómo calcular la intersección de 2 círculos?

Question

Estoy tratando de averiguar cómo derivar matemáticamente los puntos comunes de dos círculos que se cruzan en la superficie de la tierra dado un centro Lat/Lon y un radio para cada punto.

Ex.

Dada:

• Lat/Lon (37.673442, -90.234036) Radio 107.5 NM
• Lat/Lon (36.109997, -90.953669) Radio 145 NM

Debería encontrar dos puntos de intersección siendo uno de ellos (36,948, -088,158).

Sería trivial resolver esto en un plano, pero no tengo experiencia en resolver ecuaciones en una esfera imperfecta como la superficie de la Tierra. Se agradece cualquier ayuda.

Answer 1

No es mucho más difícil en la esfera que en el avión, una vez que reconoces que

1. Los puntos en cuestión son las intersecciones mutuas de tres esferas: una esfera centrada bajo el lugar x1 (en la superficie de la tierra) de un radio determinado, una esfera centrada bajo el lugar x2 (en la superficie de la tierra) de un radio determinado, y la propia tierra, que es una esfera centrada en O = (0,0,0) de un radio determinado.

2. La intersección de cada una de las dos primeras esferas con la superficie terrestre es un círculo, que define dos planos. Por tanto, las intersecciones mutuas de las tres esferas se encuentran en la intersección de esos dos planos: a línea .

En consecuencia, el problema se reduce a la intersección de una línea con una esfera, lo cual es fácil.

---

Aquí están los detalles. Las entradas son los puntos P1 = (lat1, lon1) y P2 = (lat2, lon2) en la superficie terrestre, considerada como una esfera, y dos radios correspondientes r1 y r2.

1. Convertir (lat, lon) en coordenadas geocéntricas (x,y,z). Como siempre, porque podemos elegir unidades de medida en las que la tierra tiene un radio unitario,

   x = cos(lon) cos(lat)
   y = sin(lon) cos(lat)
   z = sin(lat).

   En el ejemplo, P1 = (-90,234036 grados, 37,673442 grados) tiene coordenadas geocéntricas x1 = (-0,00323306, -0,7915, 0,61116) y P2 = (-90,953669 grados, 36,109997 grados) tiene coordenadas geocéntricas x2 = (-0,0134464, -0,807775, 0,589337).

2. Convierte los radios r1 y r2 (que se miden a lo largo de la esfera) en ángulos a lo largo de la esfera. Por definición, una milla náutica (NM) es 1/60 grados de arco (que es pi/180 * 1/60 = 0,0002908888 radianes). Por lo tanto, como ángulos,

   r1 = 107.5 / 60 Degree = 0.0312705 radian
   r2 = 145 / 60 Degree = 0.0421788 radian

3. El geodésico círculo de radio r1 alrededor de x1 es la intersección de la superficie terrestre con un Euclidiano esfera de radio sin(r1) centrada en cos(r1)*x1.

4. El plano determinado por la intersección de la esfera de radio sin(r1) alrededor de cos(r1)*x1 y la superficie terrestre es perpendicular a x1 y pasa por el punto cos(r1)*x1, por lo que su ecuación es x.x1 = cos(r1) (la "." representa la producto puntual habitual ); lo mismo para el otro plano. Habrá un único punto x0 en la intersección de esos dos planos que es una combinación lineal de x1 y x2. Escribiendo x0 = a*x1 + b*x2 las dos ecuaciones planas son

   cos(r1) = x.x1 = (a*x1 + b*x2).x1 = a + b*(x2.x1)
   cos(r2) = x.x2 = (a*x1 + b*x2).x2 = a*(x1.x2) + b

   Utilizando el hecho de que x2.x1 = x1.x2, que escribiré como q, la solución (si existe) viene dada por

   a = (cos(r1) - cos(r2)*q) / (1 - q^2),
   b = (cos(r2) - cos(r1)*q) / (1 - q^2).

   En el ejemplo actual, calculo a = 0,973503 y b = 0,0260194.

   Evidentemente necesitamos q^2 != 1. Esto significa que x1 y x2 no pueden ser ni el mismo punto ni puntos antípodas.

5. Ahora todos los demás puntos de la línea de intersección de los dos planos difieren de x0 en algún múltiplo de un vector n que es mutuamente perpendicular a ambos planos. El producto cruzado

   n = x1~Cross~x2

   hace el trabajo siempre que n sea distinto de cero: una vez más, esto significa que x1 y x2 no son ni coincidentes ni diametralmente opuestos. (Hay que tener cuidado de calcular el producto cruzado con mucha precisión, porque implica sustracciones con mucha cancelación cuando x1 y x2 están cerca el uno del otro). En el ejemplo, n = (0,0272194, -0,00631254, -0,00803124).

6. Por lo tanto, buscamos hasta dos puntos de la forma x0 + t*n que se encuentren en la superficie terrestre: es decir, su longitud es igual a 1. De forma equivalente, sus al cuadrado longitud es 1:

   1 = squared length = (x0 + t*n).(x0 + t*n) = x0.x0 + 2t*x0.n + t^2*n.n = x0.x0 + t^2*n.n

   El término con x0.n desaparece porque x0 (siendo una combinación lineal de x1 y x2) es perpendicular a n. Las dos soluciones fácilmente son

   t = sqrt((1 - x0.x0)/n.n)

   y su negativo. Una vez más se requiere alta precisión, porque cuando x1 y x2 están cerca, x0.x0 es muy cercano a 1, lo que conlleva una cierta pérdida de precisión en coma flotante. En el ejemplo, t = 1,07509 o t = -1,07509. Por lo tanto, los dos puntos de intersección son iguales a

   x0 + t*n = (0.0257661, -0.798332, 0.601666)
   x0 - t*n = (-0.0327606, -0.784759, 0.618935)

7. Por último, podemos volver a convertir estas soluciones en (lat, lon) convirtiendo las coordenadas geocéntricas (x,y,z) en coordenadas geográficas:

   lon = ArcTan(x,y)
   lat = ArcTan(Sqrt[x^2+y^2], z)

   Para la longitud, utilice la arctangente generalizada que devuelve valores en el rango de -180 a 180 grados (en aplicaciones informáticas, esta función toma ambos x e y como argumentos en lugar de sólo la relación y/x; a veces se llama "ATan2").

   Obtengo las dos soluciones (-88,151426, 36,989311) y (-92,390485, 38,238380), mostradas en la figura como puntos amarillos.

Los ejes muestran las coordenadas geocéntricas (x,y,z). La mancha gris es la parte de la superficie terrestre comprendida entre -95 y -87 grados de longitud y entre 33 y 40 grados de latitud (marcada con una retícula de un grado). La superficie terrestre se ha hecho parcialmente transparente para mostrar las tres esferas. La corrección de las soluciones calculadas es evidente por la forma en que los puntos amarillos se sitúan en las intersecciones de las esferas.

Answer 2

El elipsoidal caso:

Este problema es una generalización del de encontrar fronteras marítimas marítimos definidos como "líneas medianas" y hay una extensa literatura sobre este tema. Mi solución a este problema es aprovechar la proyección azimutal equidistante:

1. Adivina el punto de intersección
2. Proyectar los dos puntos base utilizando este punto de intersección adivinado como el centro de una proyección azimutal equidistante,
3. Resolver el problema de intersección en el espacio proyectado 2d.
4. Si el nuevo punto de intersección está demasiado lejos del anterior, vuelva al paso 2.

Este algoritmo converge cuadráticamente y da una solución precisa en un elipsoide. (La precisión es necesaria en el caso de las fronteras marítimas fronteras, ya que determina los derechos de pesca, petróleo y minerales).

Las fórmulas figuran en la sección 14 de Geodésicas en un elipsoide de revolución . La proyección azimutal equidistante elipsoidal es proporcionada por GeographicLib . Una versión de MATLAB está disponible en Proyecciones geodésicas de un elipsoide .

Answer 3

Aquí hay un código R para hacer esto:

p1 <- cbind(-90.234036, 37.673442) 
p2 <- cbind(-90.953669, 36.109997 )

library(geosphere)
steps <- seq(0, 360, 0.1)
c1 <- destPoint(p1, steps, 107.5 * 1852)
c2 <- destPoint(p2, steps, 145 * 1852)

library(raster)
s1 <- spLines(c1)
s2 <- spLines(c2)

i <- intersect(s1, s2)
coordinates(i)

#        x        y
# -92.38241 38.24267
# -88.15830 36.98740

s <- bind(s1, s2)
crs(s) <- "+proj=longlat +datum=WGS84"
plot(s)
points(i, col='red', pch=20, cex=2)

Answer 4

A raíz de Respuesta de @whuber Aquí hay un código Java que es útil por dos razones:

• pone de manifiesto un problema en relación con ArcTan (para Java, y tal vez otros lenguajes)
• maneja los posibles casos de borde, incluyendo uno no mencionado en la respuesta de @whuber.

No está optimizado ni completo (he dejado fuera clases obvias como Point ), pero debería servir.

public static List<Point> intersection(EarthSurfaceCircle c1, EarthSurfaceCircle c2) {

    List<Point> intersections = new ArrayList<Point>();

    // project to (x,y,z) with unit radius
    UnitVector x1 = UnitVector.toPlanar(c1.lat, c1.lon);
    UnitVector x2 = UnitVector.toPlanar(c2.lat, c2.lon);

    // convert radii to radians:
    double r1 = c1.radius / RadiusEarth;
    double r2 = c2.radius / RadiusEarth;

    // compute the unique point x0
    double q = UnitVector.dot(x1, x2);
    double q2 = q * q;
    if (q2 == 1) {
        // no solution: circle centers are either the same or antipodal
        return intersections;
    }
    double a = (Math.cos(r1) - q * Math.cos(r2)) / (1 - q2);
    double b = (Math.cos(r2) - q * Math.cos(r1)) / (1 - q2);
    UnitVector x0 = UnitVector.add(UnitVector.scale(x1, a), UnitVector.scale(x2, b));

    // we only have a solution if x0 is within the sphere - if not,
    // the circles are not touching.
    double x02 = UnitVector.dot(x0, x0);
    if (x02 > 1) {
        // no solution: circles not touching
        return intersections;
    }

    // get the normal vector:
    UnitVector n = UnitVector.cross(x1, x2);
    double n2 = UnitVector.dot(n, n);
    if (n2 == 0) {
        // no solution: circle centers are either the same or antipodal
        return intersections;
    }

    // find intersections:
    double t = Math.sqrt((1 - UnitVector.dot(x0, x0)) / n2);
    intersections.add(UnitVector.toPolar(UnitVector.add(x0, UnitVector.scale(n, t))));
    if (t > 0) {
        // there's only multiple solutions if t > 0
        intersections.add(UnitVector.toPolar(UnitVector.add(x0, UnitVector.scale(n, -t))));
    }
    return intersections;
}

Además, cabe destacar el uso de atan2 - es lo contrario de lo que se espera de la respuesta de @whuber (no sé por qué, pero funciona):

    public static Point toPolar(UnitVector a) {
        return new Point(
                Math.toDegrees(Math.atan2(a.z, Math.sqrt(a.x * a.x + a.y * a.y))),
                Math.toDegrees(Math.atan2(a.y, a.x)));          
    }

Answer 5

Código 'R' de trabajo para la respuesta de @wuhber.

P1 <- c(37.673442, -90.234036)
P2 <- c(36.109997, -90.953669) 

#1 NM nautical-mile is 1852 meters
R1 <- 107.5
R2 <- 145

x1 <- c(
  cos(deg2rad(P1[2])) * cos(deg2rad(P1[1])),  
  sin(deg2rad(P1[2])) * cos(deg2rad(P1[1])),
  sin(deg2rad(P1[1]))
);

x2 <- c(
  cos(deg2rad(P2[2])) * cos(deg2rad(P2[1])),
  sin(deg2rad(P2[2])) * cos(deg2rad(P2[1])),
  sin(deg2rad(P2[1]))
);

r1 = R1 * (pi/180) * (1/60)
r2 = R2 * (pi/180) * (1/60)

q = dot(x1,x2)
a = (cos(r1) - cos(r2) * q) / (1 - q^2)
b = (cos(r2) - cos(r1) * q)/ (1 - q^2)

n <- cross(x1,x2)

x0 = a*x1 + b*x2

t = sqrt((1 - dot(x0, x0))/dot(n,n))

point1 = x0 + (t * n)
point2 = x0 - (t * n)

lat1 = rad2deg(atan2(point1[2] ,point1[1]))
lon1= rad2deg(asin(point1[3]))
paste(lat1, lon1, sep=",")

lat2 = rad2deg(atan2(point2[2] ,point2[1]))
lon2 = rad2deg(asin(point2[3]))
paste(lat2, lon2, sep=",")