在叶子中绘制具有正确墨卡托畸变的线条

Question

我环顾四周，没有看到很多关于此问题的答案/讨论，尽管也许我缺少正确的词汇。

我正在使用 Python Folium 库在地图上绘制一些线条和多边形。对于我的用例，我需要沿着地球曲线绘制线条，这意味着由于墨卡托投影，它们会显示为弯曲的。

Folium 似乎只支持直线，我还没有找到任何关于这个问题的讨论，如果有人有任何建议我将不胜感激。

我的备份选项是独立处理每条线（将多边形分成线），然后将每条线绘制为贝塞尔曲线以实现所需的变形，但是我也没有运气确定如何根据墨卡托投影造成的畸变。

编辑：我刚刚想到，我上面使用贝塞尔曲线的建议虽然可能，但可能不是最有效的，相反，将线/多边形中的每个点转换为 x/y/z 表示可能会更好，计算两点之间等距的 N 个点，然后将它们投影到地球上并转换回经度/纬度。

任何建议/帮助将不胜感激！

Answer 1

最终追寻了我在编辑中提到的方法，似乎相当有效，我能够在以下脚本中使其工作：

import folium
import pyproj
import itertools
import math

SAMPLE_COORDINATES = [[65.614249, -50.759657], [61.336557, 72.357013], [-29.119194, 80.783022], [-10.117030, -74.530119]]

def pairwise(iterable):
    "s -> (s0, s1), (s1, s2), (s2, s3), ..."
    a, b = itertools.tee(iterable)
    next(b, None)
    return zip(a, b)

def render_coordinates(new_coords: list, old_coords: list):
    m = folium.Map(tiles="cartodb positron")

    bounds = compute_bounds(new_coords + old_coords)

    m.fit_bounds(bounds)

    folium.Polygon(
        locations=new_coords,
        color="red"
    ).add_to(m)

    folium.Polygon(
        locations=old_coords,
        color="blue"
    ).add_to(m)

    m.save("test.html")

def compute_bounds(coords: list):
    max_lat, max_lon = None, None
    min_lat, min_lon = None, None
    for (lat, lon) in coords:
        if min_lat is None or lat < min_lat:
            min_lat = lat
        if max_lat is None or lat > max_lat:
            max_lat = lat
        if min_lon is None or lat < min_lon:
            min_lon = lon
        if max_lon is None or lat > max_lon:
            max_lon = lon
    return ((min_lat, min_lon), (max_lat, max_lon))

def convert_LL_to_ECEF(ll_coords: list):
    transformer = pyproj.Transformer.from_crs(
        {"proj":'latlong', "ellps":'WGS84', "datum":'WGS84'},
        {"proj":'geocent', "ellps":'WGS84', "datum":'WGS84'},
    )
    ecef_coords = []
    for coords in ll_coords:
        x, y, z = transformer.transform(*(coords[::-1]), 0, radians=False)
        ecef_coords.append((x, y, z))
        
    return ecef_coords

def convert_ECEF_to_LL(ecef_coords: list):
    transformer = pyproj.Transformer.from_crs(
        {"proj":'geocent', "ellps":'WGS84', "datum":'WGS84'},
        {"proj":'latlong', "ellps":'WGS84', "datum":'WGS84'},
    )
    ll_coords = []
    for coords in ecef_coords:
        lon, lat, _ = transformer.transform(*coords, radians=False)
        ll_coords.append((lat, lon))

    return ll_coords

def compute_intermediary_points(ecef_coords: list, N: int):
    all_coords = []
    ecef_coords = ecef_coords + [ecef_coords[0]]
    for i, (coord_set_a, coord_set_b) in enumerate(pairwise(ecef_coords)):
        coord_set = []

        displacement = math.sqrt(
            (coord_set_a[0] - coord_set_b[0])**2 +
            (coord_set_a[1] - coord_set_b[1])**2 +
            (coord_set_a[2] - coord_set_b[2])**2)
            
        distanceBetweenPoints = displacement / (N + 1)

        for i in range(N):
            t = (distanceBetweenPoints * i) / displacement

            intermediate_coord = (
                    (1 - t) * coord_set_a[0] + t * coord_set_b[0],
                    (1 - t) * coord_set_a[1] + t * coord_set_b[1],
                    (1 - t) * coord_set_a[2] + t * coord_set_b[2]
            )

            coord_set.append(intermediate_coord)

        coord_set.append(coord_set_b)

        all_coords += coord_set

    return all_coords

def main():
    print("Processing coordinates.")

    ecef_coords = convert_LL_to_ECEF(SAMPLE_COORDINATES)

    new_coords = compute_intermediary_points(ecef_coords, 30)

    ll_coords = convert_ECEF_to_LL(new_coords)

    render_coordinates(ll_coords, SAMPLE_COORDINATES)
    

if __name__ == "__main__":
    main()

它绝对不完美，但做得不错：

希望其他人能从中受益！