Calcul de l'intersection de deux cercles?

J'essaie de comprendre comment dériver mathématiquement les points communs de deux cercles qui se croisent à la surface de la terre, étant donné un centre Lat / Lon et un rayon pour chaque point.

Par exemple, étant donné:

• Lat / Lon (37.673442, -90.234036) Rayon 107,5 NM
• Lat / Lon (36.109997, -90.953669) Rayon 145 NM

Je devrais trouver deux points d'intersection avec l'un d'eux étant (36,948, -088,158).

Il serait trivialement facile de résoudre ce problème sur un plan plat, mais je n'ai aucune expérience en résolution d'équations sur une sphère imparfaite telle que la surface de la Terre.

Ce n'est pas beaucoup plus difficile sur la sphère que dans l'avion, une fois que vous reconnaissez que

1. Les points en question sont les intersections mutuelles de trois sphères: une sphère centrée sous l'emplacement x1 (sur la surface de la terre) d'un rayon donné, une sphère centrée sous l'emplacement x2 (sur la surface de la terre) d'un rayon donné et la terre elle-même , qui est une sphère centrée sur O = (0,0,0) d'un rayon donné.

2. L'intersection de chacune des deux premières sphères avec la surface de la terre est un cercle qui définit deux plans. Les intersections mutuelles des trois sphères se situent donc à l'intersection de ces deux plans: une ligne .

Par conséquent, le problème se réduit à couper une ligne avec une sphère, ce qui est facile.

Voici les détails. Les entrées sont des points P1 = (lat1, lon1) et P2 = (lat2, lon2) à la surface de la terre, considérés comme une sphère, et deux rayons correspondants r1 et r2.

1. Convertissez (lat, lon) en (x, y, z) coordonnées géocentriques. Comme d'habitude, parce que nous pouvons choisir des unités de mesure dans lesquelles la terre a un rayon unitaire,

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

   Dans l'exemple, P1 = (-90.234036 degrés, 37.673442 degrés) a des coordonnées géocentriques x1 = (-0.00323306, -0.7915, 0.61116) et P2 = (-90.953669 degrés, 36.109997 degrés) a des coordonnées géocentriques x2 = (-0.0134464, -0.807775 , 0,589337).

2. Convertissez les rayons r1 et r2 (mesurés le long de la sphère) en angles le long de la sphère. Par définition, un mille marin (NM) correspond à 1/60 degré d'arc (soit pi / 180 * 1/60 = 0,0002908888 radians). Par conséquent, en tant qu'angles,

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

3. Le cercle géodésique de rayon r1 autour de x1 est l'intersection de la surface terrestre avec une sphère euclidienne de rayon sin (r1) centrée en cos (r1) * x1.

4. Le plan déterminé par l'intersection de la sphère de rayon sin (r1) autour de cos (r1) * x1 et la surface de la terre est perpendiculaire à x1 et passe par le point cos (r1) x1, d'où son équation est x.x1 = cos (r1) (le "." représente le produit scalaire habituel ); de même pour l'autre avion. Il y aura un point unique x0 à l'intersection de ces deux plans qui est une combinaison linéaire de x1 et x2. En écrivant x0 = a x1 + b * x2 les deux équations planes sont

   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
   

   En utilisant le fait que x2.x1 = x1.x2, que j'écrirai comme q, la solution (si elle existe) est donnée par

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

   Dans l'exemple en cours d'exécution, je calcule a = 0,973503 et b = 0,0260194.

   Évidemment, nous avons besoin de q ^ 2! = 1. Cela signifie que x1 et x2 ne peuvent être ni le même point ni les points antipodaux.

5. Maintenant, tous les autres points sur la ligne d'intersection des deux plans diffèrent de x0 par un multiple d'un vecteur n qui est mutuellement perpendiculaire aux deux plans. Le produit croisé

   n = x1~Cross~x2
   

   le travail fourni n est-il non nul: encore une fois, cela signifie que x1 et x2 ne sont ni coïncidents ni diamétralement opposés. (Nous devons prendre soin de calculer le produit croisé avec une grande précision, car cela implique des soustractions avec beaucoup d'annulation lorsque x1 et x2 sont proches l'un de l'autre.) Dans l'exemple, n = (0,0272194, -0,00631254, -0,00803124) .

6. Par conséquent, nous recherchons jusqu'à deux points de la forme x0 + t * n qui se trouvent à la surface de la terre: c'est-à-dire que leur longueur est égale à 1. De manière équivalente, leur longueur au carré est 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
   

   Le terme avec x0.n disparaît car x0 (étant une combinaison linéaire de x1 et x2) est perpendiculaire à n. Les deux solutions sont facilement

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

   et son négatif. Encore une fois, une haute précision est requise, car lorsque x1 et x2 sont proches, x0.x0 est très proche de 1, ce qui entraîne une perte de précision en virgule flottante. Dans l'exemple, t = 1,07509 ou t = -1,07509. Les deux points d'intersection sont donc égaux

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

7. Enfin, nous pouvons reconvertir ces solutions en (lat, lon) en convertissant les coordonnées géocentriques (x, y, z) en coordonnées géographiques:

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

   Pour la longitude, utilisez l'arctangente généralisée qui renvoie des valeurs comprises entre -180 et 180 degrés (dans les applications informatiques, cette fonction prend à la fois x et y comme arguments plutôt que simplement le rapport y / x; elle est parfois appelée "ATan2").

   J'obtiens les deux solutions (-88.151426, 36.989311) et (-92.390485, 38.238380), représentées sur la figure par des points jaunes.

Les axes affichent les coordonnées géocentriques (x, y, z). La tache grise est la partie de la surface de la terre de -95 à -87 degrés de longitude, 33 à 40 degrés de latitude (délimitée par un réticule d'un degré). La surface de la terre a été rendue partiellement transparente pour montrer les trois sphères. La justesse des solutions calculées est évidente par la position des points jaunes aux intersections des sphères.

Le cas ellipsoïdal :

Ce problème est une généralisation de celui de la recherche de frontières maritimes définies comme des «lignes médianes» et il existe une littérature abondante sur ce sujet. Ma solution à ce problème consiste à tirer parti de la projection azimutale équidistante:

1. Devinez au point d'intersection
2. Projetez les deux points de base en utilisant ce point d'intersection deviné comme le centre d'une projection azimutale équidistante,
3. Résoudre le problème d'intersection dans l'espace projeté 2D.
4. Si le nouveau point d'intersection est trop éloigné de l'ancien, revenez à l'étape 2.

Cet algorithme converge quadratique et fournit une solution précise sur un ellipsoïde. (La précision est requise dans le cas des frontières maritimes, car elle détermine les droits de pêche, de pétrole et de minéraux.)

Les formules sont données dans la section 14 de Géodésiques sur un ellipsoïde de révolution . La projection azimutale équidistante ellipsoïdale est fournie par GeographicLib . Une version MATLAB est disponible dans les projections géodésiques pour un ellipsoïde .

Voici un code R pour ce faire:

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)

À la suite de la réponse de @ whuber , voici du code Java qui est utile pour deux raisons:

• il met en évidence un piège concernant ArcTan (pour Java, et peut-être d'autres langues?)
• il gère les éventuels cas marginaux, dont un non mentionné dans la réponse de @ whuber.

Ce n'est pas optimisé ou complet (j'ai laissé de côté des classes évidentes comme Point), mais devrait faire l'affaire.

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;
}

Notez également l'utilisation de atan2- c'est l'inverse de ce que vous attendez de la réponse de @ whuber (je ne sais pas pourquoi, mais cela fonctionne):

    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)));          
    }

Code 'R' de travail pour la réponse @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=",")

Si l'un des cercles est le Nortstar, il existe un moyen le plus simple avec la sphère unitaire.

Vous pouvez mesurer votre latitude avec Nortstar. Vous avez alors une position relative sur cette sphère. v1 (0, sin (la), cos (la)) Vous connaissez la position (angle) d'une autre étoile (star2), d'Almanach. v2 (sin (lo2) * cos (la2), sin (la2), cos (lo2) * cos (la2)) Ses vecteurs. De l'équation de la sphère.

lo2 est la longitude relative. C'est inconnu .

L'angle entre vous et star2, vous pouvez aussi le mesurer, (m) Et vous savez, le produit interne de deux unités vectorielles est cos (angle) entre les deux. cos (m) = point (v1, v2) u peut maintenant calculer la longitude relative (lo2). lo2 = acos ((cos (m) -sin (la) * sin (la2)) / (cos (la) * cos (la2)))

Après tout, vous ajoutez la longitude réelle de star2 à lo2. (ou sous, dépend de son côté ouest de vous, ou à l'est.) lo2 est maintenant votre longitude.

Désolé pour mon anglais, je n'apprends jamais cette langue.

2 choses: Northstar signifie étoile polaire.

Un autre. Parce que l'angle mesuré à l'horizontale relativement, il faut toujours une correction de 90 angles. C'est également valable pour l'angle m.

ps: moyenne de l'angle réel: position de l'étoile - correction du temps.