Wrapping great circles with Mercator maps with D3.js, Leaflet, or Mapbox

Question

The question, in short: How can I accurately project a wrapping great circle radius in Mercator using something other than the Google Maps API?

The question, in long:

So I've got a conundrum. I run a mapping application that uses Google Maps API to project huge circles onto the Mercator map — it is an attempt to show very large, accurate radii, say, on the order of 13,000 km. But I don't want to use Google Maps API anymore, because Google's new pricing scheme is insane. So I'm trying to convert the code to Leaflet, or Mapbox, or anything non-Google, and nothing can handle these circles correctly.

Here's how Google Maps API handles a geodesic circle with a 13,000 km radius centered just north of Africa:

This looks intuitively weird but is correct. The wavy pattern is caused by the circle wrapping around the Earth.

D3.js can render this correctly in an orthographic projection. So here's the same circle rendered in D3.js with d3.geo.circle() on a globe, in two rotations:

Which makes the 2D-"wavy" pattern make more sense, right? Right. I love it. Totally works for my purposes of science communication and all that.

But when I convert my code to Leaflet, it doesn't work at all. Why? Because Leaflet's circle class is not a great circle at all. Instead it seems like it is just an ellipse that is distorted a bit with latitude, but not in a true geodesic way. Same circle, same radius, same origin point, and we get this:

So wrong, so wrong! Aside from looking totally unrealistic, it's just incorrect — Australia would not be inside such a circle radius. This matters for my application! This can't do.

OK, I thought, maybe the trick is to just try and implement my own great circle class. The approach I took was to iterate over circle points as distances from an origin point, but to calculate the distances using the "Destination point given distance and bearing from start point" calculations from this very helpful website, and then project those as a polygon in Leaflet. This is what I get as a result:

This looks bad but is actually much closer to being accurate! We're getting the wave effect, which is correct. Like me, you might ask, "what's really going on here?" So I did a version that allowed me to highlight every iterated point:

And you can see quite clearly that it's correctly rendered the circle, but that the polygon is incorrectly joining it. What it ought to be doing (one might naively think) is wrapping that wave figure around the multiple instances of the Mercator map projection, not naively joining them at the top, but joining them spherically. Like this crude Photoshop rendering:

And then the polygon would "close" in a way that indicated that everything above the polygon was enclosed in it, as well.

I have no idea how to implement something like this in Leaflet, though. Or anything else for that matter. Maybe I have to somehow process the raw SVG myself, taking into account the zoom state? Or something? Before I go off into those treacherous weeds, I thought I'd ask for any suggestions/ideas/etc. Maybe there's some more obvious approach.

Oh, and I tried one other thing: using the same d3.geo.circle constructor that worked so well in an orthographic projection for a Mercator/Leaflet projection. It produces more or less the same results as my "naive" Leaflet great circle implementation:

Which is promising, I guess. If you move the longitude of the origin point, though, the D3.js version wraps in a much weirder way (D3.js in red, my Leaflet class in turquoise):

I wouldn't be surprised if there was some way in D3.js to change how that worked, but I haven't gone fully down the D3.js rabbit hole. I was hoping that D3.js would make this "easier" (since it is the more full-formed cartographic tool than Leaflet), so I'm going to keep looking into this.

I have not tried to do this in Mapbox-gl yet (I guess that's next on the "attempt" list).

Anyway. Thanks for reading. To reiterate the question: How can I accurately project a wrapping great circle radius in Mercator using something other than the Google Maps API?

Answer 1

It is antimeridian cutting that's GeoJSON needs to be correctly drawn in leaflet or mapbox.

For d3 its simple d3.geoCircle(), for other mapping services, that does not handle antimeridian cutting you can use d3 to properly calculate input json.

The main idea is to unproject coordinates calculated by d3 back to lat-lng using same projection on its calculated, unprotected features will be splitted by antimeridian by d3.

See projection.invert()

I've developed example, run code snippet and drag the circles on d3 plot.

Here is result screenshot:

<!DOCTYPE html>
<html>
<head>
    <script src="https://d3js.org/d3.v5.min.js"></script>
    <script src="https://unpkg.com/[email protected]/dist/leaflet.js"></script>
    <link rel="stylesheet" href="https://unpkg.com/[email protected]/dist/leaflet.css"/>
</head>
<body style="margin: 0">
<svg style="display: inline-block"></svg>
<div id="map" style="display: inline-block; width: 350px; height: 260px"></div>
<script>
    var tileLayer = L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png');
    var map = L.map('map').setView([0,0], 0).addLayer(tileLayer);
    var bounds = [260, 260];
    var colorGenerator = d3.scaleOrdinal(d3.schemeCategory10);
    var projection = d3.geoMercator().translate([bounds[0] / 2, bounds[1] / 2]).scale(40);
    var geoPath = d3.geoPath().projection(projection);
    var geoCircle = d3.geoCircle();

    var svg = d3.select('svg')
        .attr("width", bounds[0])
        .attr("height", bounds[1])
        .attr("viewbox", "0 0 " + bounds[0] + " " + bounds[1])
        .append('g');

    svg.append("g")
        .append("path")
        .datum(d3.geoGraticule())
        .attr("stroke", "gray")
        .attr('d', geoPath);

    function addCircle(center, radius, color) {

        var g = svg.append("g");
        var drag = d3.drag().on("drag", dragged);
        var xy = projection(center);

        var path = g.append("path")
            .datum({
                type: "Polygon",
                coordinates: [[]],
                x: xy[0],
                y: xy[1]
            })
            .classed("zone", "true")
            .attr("fill", color)
            .attr("stroke", color)
            .attr("fill-opacity", 0.3)
            .call(drag);

        update(path.datum());

        function dragged(d) {
            g.raise();
            d.x = d3.event.x;
            d.y = d3.event.y;
            update(d)
        }

        function update(d) {
            center = projection.invert([d.x, d.y]);
            var poly = geoCircle.center(center).radius(radius)();
            d.coordinates[0] = poly.coordinates[0];
            path.attr('d', geoPath);
            d.geojson && d.geojson.remove();
            d.geojson = L.geoJSON(unproject(path.attr('d')), {
                color: color,
            }).addTo(map);
        }

        function unproject(d) {
            var features = d.toLowerCase().split('z').join('').split('m');
            features.shift();
            var coords = features.map(function (feature) {
                return feature.split('l').map(function (pt) {
                    var xy = pt.split(',');
                    return projection.invert([+xy[0], +xy[1]]);
                });
            });
            return {
                type: 'MultiPolygon',
                coordinates: [coords]
            }
        }
    }

    d3.range(0, 4).forEach(function (i) {
        addCircle([-120 + i * 60, 0], i * 10 + 10, colorGenerator(i));
    });
</script>
</body>
</html>

following function outputs geojson with features splitted by +-180 meridian, argument is svg path's 'd' attribute, calculated by d3:

function unproject(d, projection) {
    var features = d.toLowerCase().split('z').join('').split('m');
    features.shift();
    var coords = features.map(function (feature) {
        return feature.split('l').map(function (pt) {
            var xy = pt.split(',');
            return projection.invert([+xy[0], +xy[1]]);
        });
    });
    return {
        type: 'MultiPolygon',
        coordinates: [coords]
    }
}

Also this effect can be achieved with d3-geo-projection extesion with following code:

function unproject(geojson) {
    var projected = d3.geoProject(geojson, projection);
    if (projected.type === "MultiPolygon") {
        projected.coordinates = projected.coordinates.map(function(arr) {
            return [invert(arr[0])];
        });
    } else {
        projected.coordinates[0] = invert(projected.coordinates[0]);
    }
    return projected;
}

function invert(coords) {
    return coords.map(function(c) {
        return projection.invert(c);
    });
}

Both approaches is not handle polygons with holes, but point transformations will be same in other cases

Thank you for reading!