求地球上两个圆的交点坐标?

Find intersection coordinates of two circles on earth?

我正在寻找两个圆的第二个交点。我已经知道的其中一个点用于计算距离,然后用作圆半径(exemple)。问题是我没有得到知识点,我得到了两个新坐标,即使它们很相似。这个问题可能与地球曲率有关,但我已经搜索了一些解决方案,但一无所获。

圆的半径是用地球曲率计算的。这是我的代码:

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 变量应该用地球曲率调整其他计算。

感谢您的帮助。

你对 R 等的计算是错误的,因为平面毕达哥拉斯公式不适用于球面三角学(例如 - 我们可以在球体上得到所有三个直角的三角形!)。相反,我们应该使用特殊的公式。其中一些取自this page.

首先使用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);

现在看看我的 answer considering equal radii case 图片。 (这里的圆半径可能不同,我们应该使用另一种方法来找到所需的圆弧)

我们有球面直角三角形ACB和FCB(类似于平面情况BD在C点垂直于AF,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 沿方位 brng-CABbrng+CAB

的弧距 a1 处
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 以弧度为单位

我有类似的需求(Intersection coordinates (lat/lon) of two circles (given the coordinates of the center and the radius) on earth),特此我在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)
'''

欢迎任何反馈。