¿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

¡Muchas gracias @whuber por tu brillante trabajo! Aquí comparto el código python que escribí basado en tu trabajo, por si puede ayudar a alguien.

'''
FINDING THE INTERSECTION COORDINATES (LAT/LON) OF TWO CIRCLES (GIVEN THE COORDINATES OF THE CENTER AND THE RADII)

The code below is based on whuber's brilliant work here:
https://gis.stackexchange.com/questions/48937/calculating-intersection-of-two-circles 

The idea is that;
  1. The points in question are the mutual intersections of three spheres: a sphere centered beneath location x1 (on the 
  earth's surface) of a given radius, a sphere centered beneath location x2 (on the earth's surface) of a given radius, and
  the earth itself, which is a sphere centered at O = (0,0,0) of a given radius.
  2. The intersection of each of the first two spheres with the earth's surface is a circle, which defines two planes.
  The mutual intersections of all three spheres therefore lies on the intersection of those two planes: a line.
  Consequently, the problem is reduced to intersecting a line with a sphere.

Note that "Decimal" is used to have higher precision which is important if the distance between two points are a few
meters.
'''
from decimal import Decimal
from math import cos, sin, sqrt
import math
import numpy as np

def intersection(p1, r1_meter, p2, r2_meter):
    # p1 = Coordinates of Point 1: latitude, longitude. This serves as the center of circle 1. Ex: (36.110174,  -90.953524)
    # r1_meter = Radius of circle 1 in meters
    # p2 = Coordinates of Point 2: latitude, longitude. This serves as the center of circle 1. Ex: (36.110174,  -90.953524)
    # r2_meter = Radius of circle 2 in meters
    '''
    1. Convert (lat, lon) to (x,y,z) geocentric coordinates.
    As usual, because we may choose units of measurement in which the earth has a unit radius
    '''
    x_p1 = Decimal(cos(math.radians(p1[1]))*cos(math.radians(p1[0])))  # x = cos(lon)*cos(lat)
    y_p1 = Decimal(sin(math.radians(p1[1]))*cos(math.radians(p1[0])))  # y = sin(lon)*cos(lat)
    z_p1 = Decimal(sin(math.radians(p1[0])))                           # z = sin(lat)
    x1 = (x_p1, y_p1, z_p1)

    x_p2 = Decimal(cos(math.radians(p2[1]))*cos(math.radians(p2[0])))  # x = cos(lon)*cos(lat)
    y_p2 = Decimal(sin(math.radians(p2[1]))*cos(math.radians(p2[0])))  # y = sin(lon)*cos(lat)
    z_p2 = Decimal(sin(math.radians(p2[0])))                           # z = sin(lat)
    x2 = (x_p2, y_p2, z_p2)
    '''
    2. Convert the radii r1 and r2 (which are measured along the sphere) to angles along the sphere.
    By definition, one nautical mile (NM) is 1/60 degree of arc (which is pi/180 * 1/60 = 0.0002908888 radians).
    '''
    r1 = Decimal(math.radians((r1_meter/1852) / 60)) # r1_meter/1852 converts meter to Nautical mile.
    r2 = Decimal(math.radians((r2_meter/1852) / 60))
    '''
    3. The geodesic circle of radius r1 around x1 is the intersection of the earth's surface with an Euclidean sphere
    of radius sin(r1) centered at cos(r1)*x1.

    4. The plane determined by the intersection of the sphere of radius sin(r1) around cos(r1)*x1 and the earth's surface
    is perpendicular to x1 and passes through the point cos(r1)x1, whence its equation is x.x1 = cos(r1)
    (the "." represents the usual dot product); likewise for the other plane. There will be a unique point x0 on the
    intersection of those two planes that is a linear combination of x1 and x2. Writing x0 = ax1 + b*x2 the two planar
    equations are;
       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
    Using the fact that x2.x1 = x1.x2, which I shall write as q, the solution (if it exists) is given by
       a = (cos(r1) - cos(r2)*q) / (1 - q^2),
       b = (cos(r2) - cos(r1)*q) / (1 - q^2).
    '''
    q = Decimal(np.dot(x1, x2))

    if q**2 != 1 :
        a = (Decimal(cos(r1)) - Decimal(cos(r2))*q) / (1 - q**2)
        b = (Decimal(cos(r2)) - Decimal(cos(r1))*q) / (1 - q**2)
        '''
        5. Now all other points on the line of intersection of the two planes differ from x0 by some multiple of a vector
        n which is mutually perpendicular to both planes. The cross product  n = x1~Cross~x2  does the job provided n is 
        nonzero: once again, this means that x1 and x2 are neither coincident nor diametrically opposite. (We need to 
        take care to compute the cross product with high precision, because it involves subtractions with a lot of
        cancellation when x1 and x2 are close to each other.)
        '''
        n = np.cross(x1, x2)
        '''
        6. Therefore, we seek up to two points of the form x0 + t*n which lie on the earth's surface: that is, their length
        equals 1. Equivalently, their squared length is 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
        '''
        x0_1 = [a*f for f in x1]
        x0_2 = [b*f for f in x2]
        x0 = [sum(f) for f in zip(x0_1, x0_2)]
        '''
          The term with x0.n disappears because x0 (being a linear combination of x1 and x2) is perpendicular to n.
          The two solutions easily are   t = sqrt((1 - x0.x0)/n.n)    and its negative. Once again high precision
          is called for, because when x1 and x2 are close, x0.x0 is very close to 1, leading to some loss of
          floating point precision.
        '''
        if (np.dot(x0, x0) <= 1) & (np.dot(n,n) != 0): # This is to secure that (1 - np.dot(x0, x0)) / np.dot(n,n) > 0
            t = Decimal(sqrt((1 - np.dot(x0, x0)) / np.dot(n,n)))
            t1 = t
            t2 = -t

            i1 = x0 + t1*n
            i2 = x0 + t2*n
            '''
            7. Finally, we may convert these solutions back to (lat, lon) by converting geocentric (x,y,z) to geographic
            coordinates. For the longitude, use the generalized arctangent returning values in the range -180 to 180
            degrees (in computing applications, this function takes both x and y as arguments rather than just the
            ratio y/x; it is sometimes called "ATan2").
            '''

            i1_lat = math.degrees( math.asin(i1[2]))
            i1_lon = math.degrees( math.atan2(i1[1], i1[0] ) )
            ip1 = (i1_lat, i1_lon)

            i2_lat = math.degrees( math.asin(i2[2]))
            i2_lon = math.degrees( math.atan2(i2[1], i2[0] ) )
            ip2 = (i2_lat, i2_lon)
            return [ip1, ip2]
        elif (np.dot(n,n) == 0):
            return("The centers of the circles can be neither the same point nor antipodal points.")
        else:
            return("The circles do not intersect")
    else:
        return("The centers of the circles can be neither the same point nor antipodal points.")

'''
Example: The output of below is  [(36.989311051533505, -88.15142628069133), (38.2383796094578, -92.39048549120287)]

         intersection_points = intersection((37.673442, -90.234036), 107.5*1852, (36.109997, -90.953669), 145*1852)
         print(intersection_points)
'''

Se agradece cualquier comentario.