Obliczanie przecięcia dwóch kół?

Próbuję dowiedzieć się, jak matematycznie wyprowadzić wspólne punkty dwóch przecinających się kół na powierzchni Ziemi, biorąc pod uwagę środek Lat / Lon i promień dla każdego punktu.

Na przykład biorąc pod uwagę:

• Lat / Lon (37,673442, -90,234036) Promień 107,5 NM
• Lat / Lon (36.109997, -90.953669) Promień 145 NM

Powinienem znaleźć dwa punkty przecięcia z jednym z nich (36,948, -088.158).

Rozwiązanie problemu na płaskiej płaszczyźnie byłoby banalnie proste, ale nie mam doświadczenia w rozwiązywaniu równań na niedoskonałej kuli, takiej jak powierzchnia ziemi.

Gdy tylko to rozpoznasz, nie będzie to trudniejsze na kuli niż w samolocie

1. Chodzi o wzajemne przecięcia trzech sfer: kuli wyśrodkowanej poniżej położenia x1 (na powierzchni ziemi) o danym promieniu, kuli wyśrodkowanej poniżej położenia x2 (na powierzchni ziemi) o danym promieniu i samej ziemi , która jest kulą wyśrodkowaną na O = (0,0,0) o danym promieniu.

2. Przecięcie każdej z dwóch pierwszych sfer z powierzchnią ziemi to okrąg, który definiuje dwie płaszczyzny. Wzajemne przecięcia wszystkich trzech sfer leżą zatem na przecięciu tych dwóch płaszczyzn: linii .

W rezultacie problem sprowadza się do przecięcia linii ze sferą, co jest łatwe.

Oto szczegóły. Wejściami są punkty P1 = (lat1, lon1) i P2 = (lat2, lon2) na powierzchni ziemi, uważane za kulę, i dwa odpowiadające im promienie r1 i r2.

1. Konwertuj (lat, lon) na (x, y, z) współrzędne geocentryczne. Jak zwykle, ponieważ możemy wybrać jednostki miary, w których Ziemia ma promień jednostki,

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

   W przykładzie P1 = (-90,234036 stopnia, 37,673442 stopień) ma współrzędne geocentryczne x1 = (-0,00323306, -0,7915, 0,61116), a P2 = (-90.953669 stopień, 36,109997 stopnia) ma współrzędne geocentryczne x2 = (-0,0134464, -0,8077 , 0,589337).

2. Przelicz promienie r1 i r2 (mierzone wzdłuż kuli) na kąty wzdłuż kuli. Z definicji jedna mila morska (NM) to 1/60 stopnia łuku (co oznacza pi / 180 * 1/60 = 0,0002908888 radianów). Dlatego jako kąty

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

3. Geodezyjnej okrąg o promieniu r1 około X1 przecięcie powierzchni Ziemi z euklidesowej kuli o promieniu sin (R1) wyśrodkowany cos (R1) * X1.

4. Płaszczyzna wyznaczona przez przecięcie kuli o promieniu sin (r1) wokół cos (r1) * x1, a powierzchnia ziemi jest prostopadła do x1 i przechodzi przez punkt cos (r1) x1, skąd jej równanie wynosi x.x1 = cos (r1) („.” oznacza zwykły iloczyn skalarny ); podobnie w przypadku drugiej płaszczyzny. Na przecięciu tych dwóch płaszczyzn będzie unikalny punkt x0, który jest liniową kombinacją x1 i x2. Zapisując x0 = a x1 + b * x2 są dwa równania płaskie

   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
   

   Biorąc pod uwagę fakt, że x2.x1 = x1.x2, które napiszę jako q, rozwiązanie (jeśli istnieje) jest podane przez

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

   W bieżącym przykładzie obliczam a = 0,973503 ib = 0,0260194.

   Najwyraźniej potrzebujemy q ^ 2! = 1. Oznacza to, że x1 i x2 nie mogą być ani tym samym punktem, ani antypodalnymi punktami.

5. Teraz wszystkie inne punkty na linii przecięcia dwóch płaszczyzn różnią się od x0 pewną wielokrotnością wektora n, który jest wzajemnie prostopadły do ​​obu płaszczyzn. Produkt krzyżowy

   n = x1~Cross~x2
   

   wykonuje zadanie pod warunkiem, że n jest niezerowe: ponownie oznacza to, że x1 i x2 nie są ani zbieżne ani diametralnie przeciwne. (Musimy zadbać o obliczenie iloczynu krzyżowego z wysoką precyzją, ponieważ obejmuje on odejmowanie z dużą ilością anulowania, gdy x1 i x2 są blisko siebie.) W przykładzie n = (0,0272194, -0,00631254, -0,00803124) .

6. Dlatego szukamy do dwóch punktów formy x0 + t * n, które leżą na powierzchni ziemi: to znaczy, że ich długość wynosi 1. Odpowiednio ich kwadratowa długość wynosi 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
   

   Termin z x0.n znika, ponieważ x0 (będący liniową kombinacją x1 i x2) jest prostopadły do ​​n. Oba rozwiązania łatwo są

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

   i jego negatywne. Ponownie wymagana jest wysoka precyzja, ponieważ gdy x1 i x2 są blisko, x0.x0 jest bardzo bliskie 1, co prowadzi do pewnej utraty precyzji zmiennoprzecinkowej. W przykładzie t = 1,07509 lub t = -1,07509. Dwa punkty przecięcia są zatem równe

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

7. Na koniec możemy przekonwertować te rozwiązania z powrotem na (lat, lon), przekształcając dane geocentryczne (x, y, z) na współrzędne geograficzne:

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

   Dla długości użyć uogólnionego tangens powrocie wartości w zakresie od -180 do 180 stopni (w zastosowaniach komputerowych, funkcja wykonuje zarówno X i Y, argumentów, a nie tylko w stosunku y / x, jest czasami nazywana „ATAN2”).

   Otrzymuję dwa rozwiązania (-88,151426, 36,989311) i (-92,390485, 38,238380), pokazane na rysunku jako żółte kropki.

Osie wyświetlają współrzędne geocentryczne (x, y, z). Szara łata to część powierzchni ziemi od -95 do -87 stopni długości geograficznej, od 33 do 40 stopni szerokości geograficznej (zaznaczona siatką o jeden stopień). Powierzchnia ziemi została częściowo przezroczysta, aby pokazać wszystkie trzy sfery. Poprawność obliczonych rozwiązań jest widoczna dzięki temu, jak żółte punkty znajdują się na przecięciach sfer.

Elipsoidalny przypadek:

Problem ten jest uogólnieniem problemu znajdowania granic morskich określonych jako „linie środkowe” i istnieje obszerna literatura na ten temat. Moim rozwiązaniem tego problemu jest wykorzystanie równoodległej projekcji azymutalnej:

1. Zgadnij w punkcie przecięcia
2. Rzutuj dwa punkty bazowe, używając tego odgadniętego punktu przecięcia jako środka równoodległego rzutu azymutalnego,
3. Rozwiąż problem przecięcia w rzutowanej przestrzeni 2d.
4. Jeśli nowy punkt przecięcia znajduje się za daleko od starego, wróć do kroku 2.

Algorytm ten zbiega się kwadratowo i daje dokładne rozwiązanie dla elipsoidy. (Dokładność jest wymagana w przypadku granic morskich, ponieważ określa prawa do połowów, ropy naftowej i minerałów).

Wzory podano w rozdziale 14 Geodezji na elipsoidzie obrotu . Elipsoidalna równoodległa projekcja azymutalna jest dostarczana przez GeographicLib . Wersja MATLAB jest dostępna w rzutach geodezyjnych dla elipsoidy .

Oto trochę kodu R, aby to zrobić:

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)

Po odpowiedzi na @ whuber , oto kod Java, który jest użyteczny z dwóch powodów:

• podkreśla gotcha dotyczącą ArcTan (dla Javy i może innych języków?)
• obsługuje możliwe przypadki krawędzi, w tym jeden nie wymieniony w odpowiedzi @ whuber.

Nie jest zoptymalizowany ani ukończony (pominąłem oczywiste klasy takie jak Point), ale powinien załatwić sprawę.

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

Co ważne, zwróć uwagę na użycie atan2- jest to odwrotność tego, czego można oczekiwać od odpowiedzi @ whubera (nie wiem dlaczego, ale działa):

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

Działający kod „R” dla odpowiedzi @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=",")

Jeśli jednym z kół jest Nortstar, to najłatwiejszym sposobem jest użycie kuli jednostkowej.

Możesz zmierzyć swoją szerokość geograficzną za pomocą Nortstar. Następnie masz względną pozycję na tej sferze. v1 (0, sin (la), cos (la)) Znasz pozycję (kąt) innej gwiazdy (star2), z almanach. v2 (sin (lo2) * cos (la2), sin (la2), cos (lo2) * cos (la2)) Jego wektory. Z równania sfery.

lo2 jest względną długością geograficzną. Nie wiadomo .

Kąt między tobą a gwiazdą 2 również możesz zmierzyć, (m) I wiesz, iloczyn wewnętrzny wektora dwóch jednostek to cos (kąt) pomiędzy. cos (m) = kropka (v1, v2) możesz teraz obliczyć względną długość geograficzną (lo2). lo2 = acos ((cos (m) -sin (la) * sin (la2)) / (cos (la) * cos (la2)))

W końcu dodajesz prawdziwą długość geograficzną star2 do lo2. (lub okręt podwodny, zależą od jego zachodniej strony od ciebie lub wschodniej.) lo2 teraz jest twoją długością geograficzną.

Przepraszam za mój angielski, nigdy nie uczę się tego języka.

2 rzeczy: Northstar oznacza gwiazdę polarną.

Inne. Ponieważ kąt mierzony względem poziomu jest konieczny, zawsze potrzebna jest korekcja 90-kątowa. Obowiązuje również kąt m.

ps: średni kąt rzeczywisty: pozycja gwiazdy - korekta czasowa.