Plot lines in folium with correct Mercator distortion

Question

I’ve been looking around and not seen many answers / discussion about this, though maybe I am missing the right vocabulary.

I am using the Python Folium library to plot some lines and polygons over a map. For my use-case, I need the lines to be drawn along the earths curve, meaning they would show up as curved due to the Mercator projection.

Folium seems to only support straight lines and I haven’t been able to find any discussion about this issue, if anyone has any advice I would greatly appreciate it.

My backup option was to treat each line independently (break the polygons up into lines), then draw each line as a bezier curve to achieve the needed distortion, however I have also not had luck identifying how to calculate Bezier control points based on the distortion caused by the Mercator projection.

Edit: It just occurred to me that my above suggestion of using a Bezier curve, while possible, is probably not the most efficient, instead it would probably be better to convert each point in the line / polygon into a x/y/z representation, calculate N points equally spaced between two points, then project those onto the globe and convert back to lon/lat.

Any advice / help would be greatly appreciated!

Answer 1

Ended up chasing down the approach I mentioned in the Edit, seems to be fairly effective, I was able to get it working in the following script:

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()

Its definitely not perfect, but does a decent job:

Hopefully someone else benefits from this!