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

Calculate arc coordinates based on original point, end point, center, distance and bearing using python

有以下信息:

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

我希望能够用尽可能少的点生成一个尽可能接近的圆弧,得到这个近似的坐标。

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

我们必须记住,我们正在处理坐标,不能忽略表面曲率。

预期输出将是一个函数,该函数获取原点、终点、中心点、距离、方位角和可选的误差容限,returns作为输出一系列点坐标原点到终点大致形成所需的圆弧。

相关链接:

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

如有任何帮助,我们将不胜感激。

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)