我本善良8584314
问题描述
我在 python 方面没有那么丰富的经验，但感谢这个社区，我改进了它！我迫切需要一个接受输入并给出以下输出的函数:
I am not that experienced in python but improving it thanks to this community! I desperately need a function which takes the input and gives the ouput below:
输入:
1- 圆心 1 的经纬度坐标(例如 (50.851295, 5.667969) )
1- Latitude/longitude coordinates of the center of circle 1 (e.g. (50.851295, 5.667969) )
2- 以米为单位的圆 1 的半径(例如 200)
2- The radius of circle 1 in meters (e.g. 200)
3- 圆心 2 的经纬度坐标(例如 (50.844101, 5.725889) )
3- Latitude/longitude coordinates of the center of circle 2 (e.g. (50.844101, 5.725889) )
4- 以米为单位的圆 2 的半径(例如 300)
4- The radius of circle 2 in meters (e.g. 300)
输出:可能的输出示例可以是；
Output: Possible output examples can be;
- 交点为 (50.848295, 5.707969) 和 (50.849295, 5.717969)
- 圆圈重叠
- 圆相切，交点为 (50.847295, 5.705969)
- 圆圈不相交
我检查了这个平台、其他平台、库中的类似主题，尝试组合不同的解决方案但未能成功.非常感谢任何帮助！
I have examined the similar topics in this platform, other platforms, libraries, tried to combine different solutions but couldn't succeed. Any help is much appreciated!
非常感谢 Ture P?lsson 解决了这个问题，他在下面发表评论并指导我了解 whuber 在此链接中的出色工作 https://gis.stackexchange.com/questions/48937/calculating-intersection-of-two-circles基于这项工作，我编写了下面的代码，据我测试它可以工作.我想在这里分享它，以防有人发现它有帮助.感谢您提供任何反馈.
The problem is solved many thanks to Ture P?lsson who commented below and directed me to whuber's brilliant work in this link https://gis.stackexchange.com/questions/48937/calculating-intersection-of-two-circles Based on that work, I wrote the code below and as far as I tested it works. I want to share it here in case someone might find it helpful. Any feedback is appreciated.
'''
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 logic and
explanation 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)
'''
推荐答案
根据您需要的精度，您可能会也可能不会将地球视为一个球体.在第二种情况下，计算变得更加复杂.
Depending on the precision you need, you may or may not consider the Earth as a sphere. In the second case, calculations become more complex.
精确测量的最佳选择当半径很小(如您的示例中)是使用投影(例如 UTM)，然后应用常见的平面欧几里得计算.
The best option for precise measurements when the radius is small (as in your example) is to use a projection (UTM for example) and then apply the common flat euclidean calculations.
我们先从https://stackoverflow.com/a/55817881/2148416复制平圆相交函数:
Let's first copy the flat circle intersection function from https://stackoverflow.com/a/55817881/2148416:
def circle_intersection(x0, y0, r0, x1, y1, r1):
d = math.sqrt((x1 - x0) ** 2 + (y1 - y0) ** 2)
if d > r0 + r1: # non intersecting
return None
if d < abs(r0 - r1): # one circle within other
return None
if d == 0 and r0 == r1: # coincident circles
return None
a = (r0 ** 2 - r1 ** 2 + d ** 2) / (2 * d)
h = math.sqrt(r0 ** 2 - a ** 2)
x2 = x0 + a * (x1 - x0) / d
y2 = y0 + a * (y1 - y0) / d
x3 = x2 + h * (y1 - y0) / d
y3 = y2 - h * (x1 - x0) / d
x4 = x2 - h * (y1 - y0) / d
y4 = y2 + h * (x1 - x0) / d
return (x3, y3), (x4, y4)
可以借助 utm 库.它处理了关于地球更像是椭圆体而不是球体这一事实的所有复杂情况:
The precise calculation for a small radius (up to a few kilometers) can be done in UTM coordinates with the help of the utm library. It handles all the complications regarding the fact the the Earth is more an ellipsoid than a sphere:
import utm
def geo_circle_intersection(latlon0, radius0, latlon1, radius1):
# Convert lat/lon to UTM
x0, y0, zone, letter = utm.from_latlon(latlon0[0], latlon0[1])
x1, y1, _, _ = utm.from_latlon(latlon1[0], latlon1 [1], force_zone_number=zone)
# Calculate intersections in UTM coordinates
a_utm, b_utm = circle_intersection(x0, y0, r0, x1, y1, r1)
# Convert intersections from UTM back to lat/lon
a = utm.to_latlon(a_utm[0], a_utm[1], zone, letter)
b = utm.to_latlon(b_utm[0], b_utm[1], zone, letter)
return a, b
使用您的示例(半径稍大):
Using your example (with slightly larger radii):
>>> p0 = 50.851295, 5.667969 >>> r0 = 2000 >>> p1 = 50.844101, 5.725889 >>> r1 = 3000 >>> a, b = geo_circle_intersection(p0, r0, p1, r1) >>> print(a) (50.836848562566004, 5.684869539768468) >>> print(b) (50.860635308778285, 5.692236858407678)
