Menghitung persimpangan dua Lingkaran?

Saya mencoba mencari cara untuk mendapatkan secara matematis poin-poin umum dari dua lingkaran berpotongan pada permukaan bumi yang diberi pusat Lat / Lon dan jari-jari untuk setiap titik.

Misalnya, diberikan:

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

Saya harus menemukan dua titik persimpangan dengan salah satunya adalah (36.948, -088.158).

Akan mudah untuk menyelesaikan ini pada bidang datar tapi saya tidak punya pengalaman memecahkan persamaan pada bidang yang tidak sempurna seperti permukaan bumi.

Ini tidak jauh lebih sulit di bidang daripada di pesawat, setelah Anda menyadari itu

1. Poin yang dimaksud adalah persimpangan timbal balik dari tiga bola: bola yang berpusat di bawah lokasi x1 (di permukaan bumi) dari jari-jari yang diberikan, bola yang berpusat di bawah lokasi x2 (di permukaan bumi) dari jari-jari yang diberikan, dan bumi itu sendiri , yang merupakan bola yang berpusat pada O = (0,0,0) dari radius yang diberikan.

2. Perpotongan masing-masing dari dua bola pertama dengan permukaan bumi adalah sebuah lingkaran, yang mendefinisikan dua bidang. Oleh karena itu, persimpangan timbal balik dari ketiga bola terletak di persimpangan dua bidang: garis .

Akibatnya, masalah dikurangi menjadi memotong garis dengan bola, yang mudah.

Berikut detailnya. Inputnya adalah titik P1 = (lat1, lon1) dan P2 = (lat2, lon2) di permukaan bumi, dianggap sebagai bola, dan dua jari-jari yang bersesuaian r1 dan r2.

1. Konversi (lat, lon) ke (x, y, z) koordinat geosentris. Seperti biasa, karena kita dapat memilih satuan pengukuran di mana bumi memiliki satuan jari,

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

   Pada contoh, P1 = (-90.234036 Derajat, 37.673442 Derajat) memiliki koordinat geosentris x1 = (-0.00323306, -0.7915, 0.61116) dan P2 = (-90.953669 Derajat, 36.109997 Derajat) memiliki koordinat geosentris x2 = (-0.0134464, -0.80777575 , 0,589337).

2. Konversikan jari-jari r1 dan r2 (yang diukur di sepanjang bola) menjadi sudut di sepanjang bola. Menurut definisi, satu mil laut (NM) adalah 1/60 derajat busur (yaitu pi / 180 * 1/60 = 0,0002908888 radian). Oleh karena itu, sebagai sudut,

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

3. The geodesik lingkaran radius r1 sekitar x1 adalah persimpangan permukaan bumi dengan Euclidean lingkup radius dosa (r1) berpusat di cos (r1) * x1.

4. Pesawat ditentukan oleh persimpangan bola radius sin (r1) di sekitar cos (r1) * x1 dan permukaan bumi tegak lurus dengan x1 dan melewati titik cos (r1) x1, di mana persamaannya adalah x.x1 = cos (r1) ("." mewakili produk titik biasa ); demikian juga untuk pesawat lainnya. Akan ada titik unik x0 pada persimpangan dua bidang yang merupakan kombinasi linear dari x1 dan x2. Menulis x0 = a x1 + b * x2 dua persamaan planar

   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
   

   Menggunakan fakta bahwa x2.x1 = x1.x2, yang akan saya tulis sebagai q, solusinya (jika ada) diberikan oleh

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

   Dalam contoh yang sedang berjalan, saya menghitung a = 0,973503 dan b = 0,0260194.

   Jelas kita membutuhkan q ^ 2! = 1. Ini berarti bahwa x1 dan x2 tidak bisa menjadi titik yang sama atau titik antipodal.

5. Sekarang semua titik lain pada garis perpotongan kedua bidang berbeda dari x0 dengan kelipatan beberapa vektor n yang saling tegak lurus terhadap kedua bidang. Produk silang

   n = x1~Cross~x2
   

   apakah pekerjaan yang diberikan n adalah bukan nol: sekali lagi, ini berarti bahwa x1 dan x2 tidak bertepatan atau bertolak belakang. (Kita perlu berhati-hati untuk menghitung produk silang dengan presisi tinggi, karena itu melibatkan pengurangan dengan banyak pembatalan ketika x1 dan x2 berdekatan satu sama lain.) Dalam contoh, n = (0,0272194, -0.00631254, -0.00803124) .

6. Oleh karena itu, kami mencari hingga dua titik dari bentuk x0 + t * n yang terletak di permukaan bumi: yaitu, panjangnya sama dengan 1. Secara ekivalen, panjang kuadratnya adalah 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
   

   Istilah dengan x0.n menghilang karena x0 (menjadi kombinasi linear dari x1 dan x2) adalah tegak lurus terhadap n. Dua solusi mudah adalah

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

   dan negatifnya. Sekali lagi presisi tinggi diperlukan, karena ketika x1 dan x2 dekat, x0.x0 sangat dekat dengan 1, yang menyebabkan hilangnya presisi floating point. Dalam contoh, t = 1.07509 atau t = -1.07509. Oleh karena itu, dua titik persimpangan sama

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

7. Akhirnya, kami dapat mengonversi solusi ini kembali ke (lat, lon) dengan mengubah geosentris (x, y, z) ke koordinat geografis:

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

   Untuk bujur, gunakan arctangent umum kembali nilai-nilai dalam kisaran -180 sampai 180 derajat (dalam aplikasi komputasi, fungsi ini mengambil kedua x dan y sebagai argumen bukan hanya rasio y / x; kadang-kadang disebut "atan2").

   Saya mendapatkan dua solusi (-88.151426, 36.989311) dan (-92.390485, 38.238380), ditunjukkan dalam gambar sebagai titik kuning.

Sumbu menampilkan koordinat geosentris (x, y, z). Patch abu-abu adalah bagian dari permukaan bumi dari -95 ke -87 derajat bujur, 33 hingga 40 derajat lintang (ditandai dengan graticule satu derajat). Permukaan bumi telah dibuat sebagian transparan untuk menunjukkan ketiga bola. Kebenaran dari solusi yang dihitung terbukti dengan bagaimana titik kuning duduk di persimpangan bola.

Kasus ellipsoidal :

Masalah ini adalah generalisasi dari salah satu dari menemukan batas-batas laut yang didefinisikan sebagai "garis tengah" dan ada literatur yang luas tentang topik ini. Solusi saya untuk masalah ini adalah dengan memanfaatkan proyeksi azimut yang sama:

1. Tebak di titik persimpangan
2. Proyeksikan dua titik dasar menggunakan titik persimpangan yang dapat ditebak ini sebagai pusat dari proyeksi azimut yang sama,
3. Memecahkan masalah persimpangan di ruang proyeksi 2d.
4. Jika titik persimpangan baru terlalu jauh dari yang lama, kembali ke langkah 2.

Algoritma ini menyatu secara kuadrat dan menghasilkan solusi yang akurat pada ellipsoid. (Akurasi diperlukan dalam hal batas laut, karena menentukan hak penangkapan, minyak, dan mineral.)

Rumus diberikan dalam Bagian 14 Geodesics pada ellipsoid revolusi . Proyeksi azimut equidistant ellipsoidal disediakan oleh GeographicLib . Versi MATLAB tersedia di proyeksi Geodesic untuk ellipsoid .

Berikut adalah beberapa kode R untuk melakukan ini:

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)

Mengikuti dari jawaban @ whuber , inilah beberapa kode Java yang berguna karena dua alasan:

• itu menyoroti gotcha tentang ArcTan (untuk Java, dan mungkin bahasa lain?)
• itu menangani kemungkinan kasus tepi, termasuk yang tidak disebutkan dalam jawaban @ whuber.

Itu tidak dioptimalkan atau lengkap (saya sudah meninggalkan kelas yang jelas suka Point), tetapi harus melakukan trik.

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

Juga, yang penting, perhatikan penggunaan atan2- ini kebalikan dari apa yang Anda harapkan dari jawaban @ whuber (saya tidak tahu mengapa, tetapi berfungsi):

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

Bekerja kode 'R' untuk jawaban @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=",")

Jika salah satu lingkaran adalah Nortstar, maka ada cara termudah dengan unit sphere.

Anda dapat mengukur garis lintang Anda dengan Nortstar. Maka Anda memiliki posisi relatif di bidang ini. v1 (0, sin (la), cos (la)) Anda tahu posisi (sudut) bintang lain (star2), dari almanach. v2 (sin (lo2) * cos (la2), sin (la2), cos (lo2) * cos (la2)) vektor-vektornya. Dari persamaan bola.

lo2 adalah garis bujur relatif. Itu tidak diketahui .

Sudut antara Anda dan star2, Anda dapat mengukur juga, (m) Dan Anda tahu, produk dalam dari dua satuan vektor adalah cos (sudut) di antaranya. cos (m) = titik (v1, v2) u dapat menghitung sekarang garis bujur relatif (lo2). lo2 = acos ((cos (m) -sin (la) * sin (la2)) / (cos (la) * cos (la2)))

Setelah semua, Anda menambahkan bujur nyata dari star2 ke lo2. (atau sub, tergantung pada sisi barat dari Anda, atau timur.) lo2 sekarang adalah bujur Anda.

Maaf untuk bahasa Inggris saya, saya tidak pernah belajar bahasa ini.

2 hal: Bintang kutub utarastar.

Lain. Karena sudut yang diukur untuk mendatar secara relatif, selalu perlu koreksi 90 sudut. Ini berlaku untuk sudut m juga.

ps: sudut pandang sebenarnya: posisi bintang - koreksi waktu.