如何高效计算地图上线条与圆形区域的首个交点？

我目前用geopy计算地图上直线与圆形区域的首个交点，现有代码能正常运行，但因为要逐点计算直线上的经纬度直到进入圆的半径范围，效率很低。

请问能不能通过geopy或其他同类库实现更高效的计算方式？

示例场景：直线起点为51.137933, -0.267017，沿257度方向直行，直到距离圆心51.053953, -0.625003不超过10.5海里。

现有代码：

import math
from geopy.distance import distance as dist

# Starting position in degrees
lat1 = 51.137933
lon1 = -0.267017

target_lat = 51.053953
target_lon = -0.625003

# Direction in radians
direction = math.radians(257)

# Earth's radius in miles
R = 3963.1676

# Distance to move in miles
distance = 0.1

d = dist((lat1, lon1), (target_lat, target_lon)).nautical
last_d = d

while d > 10.5:
    # Convert latitude and longitude to radians
    lat1 = math.radians(lat1)
    lon1 = math.radians(lon1)

    # Calculate new latitude and longitude
    lat2 = math.asin(math.sin(lat1) * math.cos(distance / R) + math.cos(lat1) * math.sin(distance / R) * math.cos(direction))
    lon2 = lon1 + math.atan2(math.sin(direction) * math.sin(distance / R) * math.cos(lat1), math.cos(distance / R) - math.sin(lat1) * math.sin(lat2))

    # Convert latitude and longitude back to degrees
    lat1 = math.degrees(lat2)
    lon1 = math.degrees(lon2)

    d = dist((lat1, lon1), (target_lat, target_lon)).nautical
    if d > last_d:
        print("NOT FOUND")
        break

print(lat1, lon1)

高效解决方案：直接通过球面几何计算交点 #

逐点迭代的方式效率低，核心是做了大量冗余计算。我们可以通过球面几何公式直接推导交点位置，无需循环遍历，计算复杂度从O(n)降到O(1)。

核心思路 #

1. 计算起点到圆心的初始球面距离，以及两者连线与直线前进方向的夹角。
2. 判断直线是否会与圆相交：如果圆心到直线的最短距离大于圆半径，直接判定无交点。
3. 计算从起点到交点的精确距离，再用geopy的destination方法直接生成交点坐标。

实现代码 #

import math
from geopy.distance import distance as dist
from geopy.point import Point
from geopy.distance import Distance

# 初始参数
start_point = Point(51.137933, -0.267017)
center_point = Point(51.053953, -0.625003)
direction_deg = 257  # 前进方位角（度）
radius_nautical = 10.5  # 圆半径（海里）

# 1. 计算起点到圆心的球面距离（海里）
D = dist(start_point, center_point).nautical

# 2. 计算起点到圆心的方位角，再求与前进方向的夹角θ
start_to_center_bearing = start_point.bearing_to(center_point)
theta_rad = math.radians(abs(direction_deg - start_to_center_bearing))
# 处理夹角超过180度的情况
theta_rad = theta_rad if theta_rad <= math.pi else 2 * math.pi - theta_rad

# 3. 判断直线是否与圆相交
if D * math.sin(theta_rad) > radius_nautical:
    print("NOT FOUND")
else:
    # 4. 计算从起点到交点的距离（海里）
    d_intersect = D * math.cos(theta_rad) - math.sqrt(radius_nautical**2 - (D * math.sin(theta_rad))**2)
    
    # 若起点已在圆内，直接返回起点坐标
    if d_intersect < 0:
        print(start_point.latitude, start_point.longitude)
    else:
        # 5. 直接生成交点坐标
        intersect_point = Distance(nautical=d_intersect).destination(start_point, bearing=direction_deg)
        print(intersect_point.latitude, intersect_point.longitude)

优势说明 #

• 完全避免循环迭代，效率提升显著，尤其适合长距离或需要高频计算的场景。
• 复用geopy内置的bearing_to和destination方法，无需手动编写复杂的球面坐标转换逻辑，减少出错概率。

内容的提问来源于stack exchange，提问作者KillerKode