使用python根据原点、终点、圆心、距离、方位计算圆弧坐标

Question

拥有以下信息：

• 原点：点(lat_origin, long_origin)
• 终点：点(lat_end, long_end)
• 中心点：点(lat_center, long_center)
• 距离：100
• 轴承：90°

from shapely.geometry import Point
origin_point = Point(...,...)
end_point = Point(...,...)
center_point = Point(...,...)
distance = 100
bearing = 90

我希望能够用尽可能少的点生成尽可能近的圆弧，并获得该近似值的坐标。

一个好的功能是能够控制误差容限并能够动态分级点数以近似弧线。

我们必须记住，我们正在使用坐标，我们不能忽视表面曲率。

预期的输出将是一个函数，该函数获取原点、终点、中心点、距离、方位以及可选的误差容限作为输入，并返回从原始点到终点的一系列点坐标作为输出大致形成所需的弧度。

相关链接：

https://gis.stackexchange.com/questions/326871/generate-arc-from-projection-cooperatives

任何帮助将不胜感激。

Answer 1

https://www.igismap.com/formula-to-find-bearing-or-heading-angle- Between-two-points-latitude-longitude/

import math
import numpy as np
from shapely.geometry import Point, LineString

def get_bearing(center_point, end_point):
    
    lat3 = math.radians(end_point[0])
    long3 = math.radians(end_point[1])
    lat1 = math.radians(center_point[0])
    long1 = math.radians(center_point[1])
    
    dLon = long3 - long1
    
    X = math.cos(lat3) * math.sin(dLon)
    Y = math.cos(lat1) * math.sin(lat3) - math.sin(lat1) * math.cos(lat3) * math.cos(dLon)
    
    end_brng = math.atan2(X, Y)
    
    return end_brng

def get_arc_coordinates(center_point, origin_point, end_point, brng_init, distance):
    '''
    center_point: (center_latitude, center_long) 
    origin_point: (origin_latitude, origin_long) 
    end_point: (end_latitude, end_long)
    brng_init: degrees
    distance: aeronautical miles
    '''
    
    brng_init = math.radians(brng_init) #Bearing in degrees converted to radians.
    d = distance * 1.852 #Distance in km
    
    R = 6378.1 #Radius of the Earth
    brng = get_bearing(center_point,end_point) #1.57 #Bearing is 90 degrees converted to radians.
    
    list_bearings = np.arange(brng, brng_init, 0.1) # 0.1 value to be modify with tolerance
    
    coordinates = []
    
    for i in list_bearings:
        lat1 = math.radians(center_point[0]) #Center lat point converted to radians
        lon1 = math.radians(center_point[1]) #Center long point converted to radians
        brng = i
        lat2 = math.asin( math.sin(lat1)*math.cos(d/R) +
             math.cos(lat1)*math.sin(d/R)*math.cos(brng))
        
        lon2 = lon1 + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat1),
                     math.cos(d/R)-math.sin(lat1)*math.sin(lat2))
        
        lat2 = math.degrees(lat2)
        lon2 = math.degrees(lon2)
        
        coordinates.append(Point(lat2, lon2))

    return LineString(coordinates)

Answer 2

def get_arc_coordinates(center_point, origin_point, end_point, start_angle, distance, step_size=0.1): ''' 中心点：（中心纬度，中心长） 原点：（原点纬度，原点长） 终点：（终点纬度，终点长） 起始角度：度 距离： 海里 step_size：度数，调整该参数可以控制沿圆弧的点数 ''' start_angle = math.radians(start_angle) # 将起始角度转换为弧度 d = 距离 * 1.852 # 距离（公里）

R = 6378.1  # Radius of the Earth

# Calculate the bearing between origin and end points
brng = get_bearing(origin_point, end_point)

# Calculate the arc coordinates starting from the given start angle
coordinates = []
while True:
    lat1 = math.radians(center_point[0])  # Center lat point converted to radians
    lon1 = math.radians(center_point[1])  # Center long point converted to radians
    
    # Calculate the point on the arc using the start angle
    lat2 = math.asin(math.sin(lat1) * math.cos(d / R) +
                      math.cos(lat1) * math.sin(d / R) * math.cos(start_angle))
    
    lon2 = lon1 + math.atan2(math.sin(start_angle) * math.sin(d / R) * math.cos(lat1),
                              math.cos(d / R) - math.sin(lat1) * math.sin(lat2))
    
    lat2 = math.degrees(lat2)
    lon2 = math.degrees(lon2)
    
    coordinates.append((lon2, lat2))
    
    # Stop when the distance from the end point is less than the tolerance
    if math.sqrt((lon2 - end_point[1])**2 + (lat2 - end_point[0])**2) < 0.01:
        break
    
    # Update start angle for the next iteration
    start_angle += math.radians(step_size)  # Adjust the step size
    
# Append the end point to the coordinates
coordinates.append((end_point[1], end_point[0]))

return coordinates

我这样做是顺时针和逆时针 起始角度 -= math.radians(step_size)。但输入的起点和计算的起点略有不同。我该如何纠正？我也遇到过重叠问题，geom 已创建