Найти координаты пересечения двух кругов на земле?
Я пытаюсь найти вторую точку пересечения двух кругов. Одна из точек, которые я уже знаю, использовалась для вычисления расстояния, а затем использовалась в качестве радиуса круга ( пример). Проблема в том, что я не знаю точку, я получаю две новые координаты, даже если они похожи. Проблема, вероятно, связана с искривлением Земли, но я искал какое-то решение и ничего не нашел.
Радиус окружностей рассчитывается с учетом кривизны земли. И вот код, который у меня есть:
function GET_coordinates_of_circles(position1,r1, position2,r2) {
var deg2rad = function (deg) { return deg * (Math.PI / 180); };
x1=position1.lng;
y1=position1.lat;
x2=position2.lng;
y2=position2.lat;
var centerdx = deg2rad(x1 - x2);
var centerdy = deg2rad(y1 - y2);
var R = Math.sqrt(centerdx * centerdx + centerdy * centerdy);
if (!(Math.abs(r1 - r2) <= R && R <= r1 + r2)) { // no intersection
console.log("nope");
return []; // empty list of results
}
// intersection(s) should exist
var R2 = R*R;
var R4 = R2*R2;
var a = (r1*r1 - r2*r2) / (2 * R2);
var r2r2 = (r1*r1 - r2*r2);
var c = Math.sqrt(2 * (r1*r1 + r2*r2) / R2 - (r2r2 * r2r2) / R4 - 1);
var fx = (x1+x2) / 2 + a * (x2 - x1);
var gx = c * (y2 - y1) / 2;
var ix1 = fx + gx;
var ix2 = fx - gx;
var fy = (y1+y2) / 2 + a * (y2 - y1);
var gy = c * (x1 - x2) / 2;
var iy1 = fy + gy;
var iy2 = fy - gy;
// note if gy == 0 and gx == 0 then the circles are tangent and there is only one solution
// but that one solution will just be duplicated as the code is currently written
return [[iy1, ix1], [iy2, ix2]];
}
Переменная deg2rad должна корректировать другие расчеты с кривизной земли.
Спасибо за любую помощь.
1 ответ
Ваши вычисления для R и т. Д. Неверны, потому что плоская формула Пифагора не работает для сферической тригонометрии (например, у нас может быть треугольник со всеми тремя прямыми углами на сфере!). Вместо этого мы должны использовать специальные формулы. Некоторые из них взяты с этой страницы.
Сначала найдите большие круговые дуги в радианах для обоих радиусов, используя R = Earth radius = 6,371km
a1 = r1 / R
a2 = r2 / R
И расстояние (снова дуга в радианах) между центром круга по формуле haversine
var R = 6371e3; // metres
var φ1 = lat1.toRadians();
var φ2 = lat2.toRadians();
var Δφ = (lat2-lat1).toRadians();
var Δλ = (lon2-lon1).toRadians();
var a = Math.sin(Δφ/2) * Math.sin(Δφ/2) +
Math.cos(φ1) * Math.cos(φ2) *
Math.sin(Δλ/2) * Math.sin(Δλ/2);
var ad = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
А подшипник из положения 1 в положение 2:
//where φ1,λ1 is the start point, φ2,λ2 the end point
//(Δλ is the difference in longitude)
var y = Math.sin(λ2-λ1) * Math.cos(φ2);
var x = Math.cos(φ1)*Math.sin(φ2) -
Math.sin(φ1)*Math.cos(φ2)*Math.cos(λ2-λ1);
var brng = Math.atan2(y, x);
Теперь посмотрите на картинку из моего ответа, рассматривая случай равных радиусов. (Здесь радиусы окружности могут быть различны, и мы должны использовать другой подход, чтобы найти необходимые дуги)
У нас есть сферические прямоугольные треугольники ACB и FCB (аналогично плоскому случаю BD перпендикулярно AF в точке C, а угол BCA - прямой).
Сферическая теорема Пифагора (из книги по sph. Trig) гласит, что
cos(AB) = cos(BC) * cos(AC)
cos(FB) = cos(BC) * cos(FC)
или (используя x для AC, y для BC и (ad-x) для FC)
cos(a1) = cos(y) * cos(x)
cos(a2) = cos(y) * cos(ad-x)
разделить уравнения, чтобы исключить cos(y)
cos(a1)*cos(ad-x) = cos(a2) * cos(x)
cos(a1)*(cos(ad)*cos(x) + sin(ad)*sin(x)) = cos(a2) * cos(x)
cos(ad)*cos(x) + sin(ad)*sin(x) = cos(a2) * cos(x) / cos(a1)
sin(ad)*sin(x) = cos(a2) * cos(x) / cos(a1) - cos(ad)*cos(x)
sin(ad)*sin(x) = cos(x) * (cos(a2) / cos(a1) - cos(ad))
TAC = tg(x) = (cos(a2) / cos(a1) - cos(ad)) / sin(ad)
Имея гипотенузу и катетус треугольника ACB, мы можем найти угол между направлениями AC и AB (правила Нейпира для правильных сферических треугольников) - заметьте, мы уже знаем TAC = tg(AC) а также a1 = AB
cos(CAB)= tg(AC) * ctg(AB)
CAB = Math.acos(TAC * ctg(a1))
Теперь мы можем рассчитать точки пересечения - они лежат на расстоянии дуги а1 от положения 1 вдоль опор brng-CAB а также brng+CAB
B_bearing = brng - CAB
D_bearing = brng + CAB
Координаты точек пересечения:
var latB = Math.asin( Math.sin(lat1)*Math.cos(a1) +
Math.cos(lat1)*Math.sin(a1)*Math.cos(B_bearing) );
var lonB = lon1.toRad() + Math.atan2(Math.sin(B_bearing)*Math.sin(a1)*Math.cos(lat1),
Math.cos(a1)-Math.sin(lat1)*Math.sin(lat2));
и то же самое для D_bearing
latB, lonB в радианах
У меня была аналогичная потребность ( координаты пересечения (широта / долгота) двух кругов (с учетом координат центра и радиуса) на Земле), и поэтому я делюсь решением на python на случай, если оно может кому-то помочь:
'''
FINDING THE INTERSECTION COORDINATES (LAT/LON) OF TWO CIRCLES (GIVEN THE COORDINATES OF THE CENTER AND THE RADII)
Many thanks to Ture Pålsson who directed me to the right source, 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)
'''
Любая обратная связь приветствуется.