无法安装Cartopy，求手动处理Shapely LineStrings大地线的替代方案

当然可以手动实现，核心思路是把两点间的笛卡尔直线替换为地球表面的大圆航线（测地线），通过生成足够多的中间点来近似这条曲线，再用这些点构建Shapely LineString。以下是两种可行方案：

方案一：用pyproj计算测地线（推荐） #

如果环境能安装pyproj（比Cartopy轻量很多），可以直接用它的Geod类计算测地线坐标点，步骤简单可靠：

1. 初始化WGS84椭球的Geod对象（对应Cartopy默认的地球模型）
2. 输入起点和终点的经纬度，生成测地线上的密集点
3. 用这些点构建Shapely LineString

示例代码：

from shapely.geometry import LineString
import pyproj

# 初始化WGS84椭球的Geod实例
geod = pyproj.Geod(ellps='WGS84')

# 起点、终点经纬度（格式：经度, 纬度）
start = (116.40, 39.90)  # 北京
end = (121.47, 31.23)    # 上海

# 生成测地线上的中间点，点越多曲线越平滑
n_points = 100
lon_list, lat_list = geod.npts(start[0], start[1], end[0], end[1], n_points)

# 拼接起点、中间点、终点，构建LineString
all_lon = [start[0]] + lon_list + [end[0]]
all_lat = [start[1]] + lat_list + [end[1]]
geodesic_line = LineString(zip(all_lon, all_lat))

方案二：纯手动近似计算（无额外依赖） #

如果无法安装任何第三方库，可通过球面三角学公式手动计算大圆航线的中间点。以下是简化实现（基于球面近似，精度略低于椭球计算）：

from shapely.geometry import LineString
import math

def calculate_geodesic_points(start_lon, start_lat, end_lon, end_lat, n_points=100):
    # 经纬度转弧度
    lon1, lat1 = math.radians(start_lon), math.radians(start_lat)
    lon2, lat2 = math.radians(end_lon), math.radians(end_lat)

    # 计算大圆距离与初始方位角
    d = math.acos(math.sin(lat1)*math.sin(lat2) + math.cos(lat1)*math.cos(lat2)*math.cos(lon2-lon1))
    bearing1 = math.atan2(math.sin(lon2-lon1)*math.cos(lat2), math.cos(lat1)*math.sin(lat2)-math.sin(lat1)*math.cos(lat2)*math.cos(lon2-lon1))

    points = []
    for i in range(n_points+1):
        fraction = i / n_points
        a = math.sin((1-fraction)*d) / math.sin(d)
        b = math.sin(fraction*d) / math.sin(d)
        # 计算三维坐标
        x = a*math.cos(lat1)*math.cos(lon1) + b*math.cos(lat2)*math.cos(lon2)
        y = a*math.cos(lat1)*math.sin(lon1) + b*math.cos(lat2)*math.sin(lon2)
        z = a*math.sin(lat1) + b*math.sin(lat2)
        # 转回经纬度
        lat = math.atan2(z, math.sqrt(x**2 + y**2))
        lon = math.atan2(y, x)
        points.append((math.degrees(lon), math.degrees(lat)))
    return points

# 示例使用
start = (116.40, 39.90)
end = (121.47, 31.23)
geodesic_points = calculate_geodesic_points(*start, *end)
geodesic_line = LineString(geodesic_points)

注意事项 #

• 生成的点数量越多，曲线越接近真实大圆航线，但计算量也会增加，可根据可视化需求调整n_points
• 手动计算基于球面模型，若对精度要求极高，优先选择pyproj的椭球计算方案

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