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?

Myeloid

For me d3.geoCircle() works just fine. Have a look at this demo forked from the one I set up for my answer to "D3 Topojson Circle with Radius Scaled in Miles". – Cralg

That's very promising! I'll check it out. I've not been using d3 v4, much less v5, and maybe that's the issue. Thank you. – Myeloid

So I'm struggling to adapt this to a Leaflet environment. To use d3 with Leaflet, the projection must be something like: var transform = d3.geoTransform({point: projectPoint}); (per examples like this one: bl.ocks.org/shimizu/749df041c1945aef78fd992c7dfbe0e1). However when you use that with svg.append("path").datum({type: "Sphere"}).attr("id", "sphere").attr("d", path); you get an error: "this.stream.sphere is not a function." If you omit the "Sphere" path, it gives exactly the same (wrong) results as before. Oh d3, and your mysterious ways... – Myeloid

Here's what I'm talking about: jsfiddle.net/9r7x2ond/1 Without the Sphere datum part, it clearly doesn't work. But with it, it throws an error. I'm either missing something or there's some fundamental incompatibility? – Myeloid

My current thinking is that this isn't a sphere projection issue, it's more fundamental again: when you render a shape like this in d3.js, it understands that the map is never going to wrap, and seems to know what the upper and lower bounds are an act accordingly. It doesn't do that in a tiled map environment: it doesn't wrap (or know how to — it really needs 3X more curves). Which is kind of my original problem again. I haven't seen anything (after much Googling) that suggests anyone has done this with Leaflet or D3 before. :-/ – Myeloid

Answer 1

So it ended up being not something with an easy solution. To achieve the desired Google Maps-like behavior, I ended up having to code a Leaflet plugin from scratch that extended the L.Polygon object. This is because the desired behavior including a "wrapping" polygon, and there's no "magic" way to do that in Leaflet.

What I ended up doing was creating a plugin that could detect whether it ought (based on the zoom level) to create numerous wrapped "copies," and then use a bit of logic to determine whether it should be splicing together polygons or not. It's not especially elegant (it is more logic than math) but that's my programming in a nutshell.

Anyway, here is the final plugin. It can be dropped in like a regular L.Circle object (just change it to L.greatCircle) without too many other changes. You can see it in action on my MISSILEMAP (which also features a geodesic polyline class I had to write, which was a lot easier).

Thanks for those who gave advice and suggestions.

Myeloid

Answer 2

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!

Bakemeier

This is interesting. I haven't fully understood the code, but it's an interesting approach to use D3 and then generate GeoJSON from it and import it to the Lefalet tileset. – Myeloid

of course you can dont render d3, and only use it for calculations – Bakemeier

@Myeloid you are welcome, ask me if you need code clarification – Bakemeier