计算两个圆的交点

#1 楼

一旦您意识到


，所讨论的点就是三个球体的相互交点，那么在球体上的难度就比在平面上难得多。一个球体位于x1位置下方（在地球的给定半径的球），以给定半径的位置x2（在地球表面上）为中心的球体和以球半径为O =（0,0,0）为中心的球体本身。 br />前两个球体与地球表面的交集是一个圆，它定义了两个平面。因此，所有三个球体的相互交点位于这两个平面的交点上：一条线。

，因此，问题简化为将一条线与一个球体相交，这很容易。


以下是详细信息。输入是地球表面上的点P1 =（lat1，lon1）和P2 =（lat2，lon2），以及两个对应的半径r1和r2。



将（lat，lon）转换为（x，y，z）地心坐标。像往常一样，因为我们可以选择地球具有单位半径的度量单位，所以在示例中，P1 =（-90.234036度，37.673442度）具有地心坐标x1 =（-0.00323306，-0.7915、0.61116）和P2 =（-90.953669度，36.109997度）具有地心坐标x2 =（-0.0134464，-0.807775、0.589337）。


将半径r1和r2（沿球体测量）转换为沿球体的角度。根据定义，一海里（NM）是1/60弧度（即pi / 180 * 1/60 = 0.0002908888弧度）。因此，作为角度，围绕x1的半径r1的测地线圆是地球表面与以cos为中心的半径sin（r1）的欧几里得球面的交点。 （r1）* x1。

由围绕cos（r1）* x1的半径sin（r1）的球面与地球表面的交点确定的平面垂直于x1并经过点cos（r1）x1，因此其等式为x.x1 = cos （r1）（“。”代表通常的点积）；对于另一架飞机也是如此。在这两个平面的交点上将有一个唯一点x0，它是x1和x2的线性组合。写x0 = ax1 + b * x2这两个平面方程是

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

利用x2.x1 = x1.x2的事实，我将其写为q，解决方案（如果存在）由

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

在运行的示例中，我计算a = 0.973503和b = 0.0260194。

显然我们需要q ^ 2！=1。这意味着x1和x2既不能是同一点，也不能是对映点。


现在将所有其他点都放在两个平面的相交线上与x0的差为向量n的倍数，向量n垂直于两个平面。叉积

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

做工作，只要n不为零：再次，这意味着x1和x2既不是重合也不是完全相反。 （我们需要小心地计算叉积，因为当x1和x2彼此接近时，它涉及很多相减的减法。）在示例中，n =（0.0272194，-0.00631254，-0.00803124） 。


因此，我们寻找位于地球表面的x0 + t * n形式的最多两个点：即它们的长度等于1。等效地，它们的平方长度为1：

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

因为x0（是x1和x2的线性组合）垂直于n，所以x0.n的项消失了。这两个解决方案很容易

n = x1~Cross~x2

及其负面影响。再次需要高精度，因为当x1和x2接近时，x0.x0非常接近1，从而导致浮点精度下降。在示例中，t = 1.07509或t = -1.07509。因此，两个相交点相等

1 = squared length = (x0 + t*n).(x0 + t*n) = x0.x0 + 2t*x0.n + t^2*n.n = x0.x0 + t^2*n.n

最后，我们可以通过转换地心线（x，y， z）到地理坐标：

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

对于经度，请使用-180至180度范围内的广义反正切返回值（在计算应用程序中，此函数同时使用x和y作为参数，而不仅仅是y / x；有时称为“ ATan2”。

我得到了两个解（-88.151426，36.989311）和（-92.390485，38.238380），在图中以黄色圆点显示。




轴显示地心（x，y，z）坐标。灰色斑块是地球表面经度-95至-87度，纬度33至40度（以1度标线标记）的部分。地球表面已经部分透明，可以显示所有三个球体。计算出的解的正确性通过黄点在球体的交点处可见。

#2 楼

椭圆形情况：

这个问题是寻找被定义为“中线”的海事边界的一个概括，关于这一主题有大量文献
。我对这个问题的解决方案是利用等距的方位角投影：


在相交点进行猜测
使用这个猜测的相交点将两个基点投影为
等距方位角投影的中心，
解决2d投影空间中的相交问题。
新相交点与旧相交点太远，返回
第2步。

该算法二次收敛，并在椭圆体上产生精确的解
。 （对于海洋
边界，要求精确度，因为它决定了捕鱼，石油和矿产权。）

公式在《大地测量学》第14节中给出了椭圆形
革命。椭球等距方位角投影
由GeographicLib提供。 MATLAB的版本可在
大地投影中获得椭球。

#3 楼

以下是一些R代码可以执行此操作：

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)

#4 楼

根据@whuber的回答，以下是一些有用的Java代码，有两个原因：


它突出显示了有关ArcTan的陷阱（适用于Java，也许还有其他语言？）
它处理可能的极端情况，包括@whuber的答案中未提及的情况。

尚未优化或完成（我遗漏了诸如Point之类的明显类），但应该可以解决问题。

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

此外，重要的是，请注意atan2的使用-这与您从@whuber答案中得到的期望相反（我不知道为什么，但是有效）：

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

#5 楼

@wuhber答案的有效“ R”代码。

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=",")

#6 楼

非常感谢@whuber的出色工作！在此，我分享了根据您的工作输入的python代码，以防它对某人有所帮助。

 '''
FINDING THE INTERSECTION COORDINATES (LAT/LON) OF TWO CIRCLES (GIVEN THE COORDINATES OF THE CENTER AND THE RADII)

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)
'''
 

任何反馈是感谢。

#7 楼

如果圆弧之一是Nortstar，则存在单位球的最简单方法。

您可以使用Nortstar测量纬度。然后，您在该球体上具有相对位置。
v1（0，sin（la），cos（la））
您可以从年历中知道另一颗恒星（star2）的位置（角度）。
v2（sin（lo2）* cos（la2），sin（la2），cos（lo2）* cos（la2））
其向量。从球面方程开始。

lo2是相对经度。它是未知数。

您和star2之间的角度也可以测量，（m）
您知道，两个单位向量的内积是cos（角度）之间的。
cos（m）= dot（v1，v2）
u现在可以计算相对经度（lo2）。
lo2 = acos（（cos（m）-sin（la）* sin（la2） ）/（cos（la）* cos（la2）））

毕竟将star2的真实经度添加到lo2中。 （或潜艇，取决于它在您的西边还是向东。）
lo2现在是您的经度。

对不起，我的英语，我永远不会学习这种语言。


2件事：
北极星指极星。

另一个。因为相对于水平角测量的角度总是需要90角校正。
对m角也有效。
ps：实际角度均值：恒星位置-时间校正。