Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to fix map boundaries on d3 cartographic raster reprojection?

I try to use the raster reprojection of a map following this example. If I change the example kavrayskiy7 projection by the Azimuthal Equidistant projection,

var projection = d3.geo.azimuthalEquidistant()
    .scale(90)
    .translate([width / 2, height / 2])
    .clipAngle(180 - 1e-3)
    .precision(.1);

it should project the Earth onto a disc (the image of the projection map). However, the raster reprojection goes beyond that disc and fills the entire canvas with an extended picture (the inverse projection function is not injective, several x/y points on the map correspond to a single lon/lat coordinates). In the original example, this should be avoided with the line

if (λ > 180 || λ < -180 || φ > 90 || φ < -90) { i += 4; continue; }

but for this example that does not work. I found other glitches for instance when using the Mollweide projection (two lines appear at the poles) due to the same efect.

To solve this, one way would be to fix the inverse projections so they return error or None when the x/y input is out of range. My attempt was to check if a point is on range using the forward projection of the whole sphere to obtain a SVG path with the boundary of the map, as given by this code:

var path = d3.geo.path()
    .projection(projection);

var bdry = svg.append("defs").append("path")
    .datum({type: "Sphere"})
    .attr("id", "sphere")
    .attr("d", path);

(see for instance this example). However, I found no easy method to check whether a point [x,y] is inside a SVG closed path.

So my questions are:

  • Is there a bug on the inverse projections, or am I not using them correctly?
  • How could I find if a [x,y] point is inside the svg path, assuming that this is the best approach?
  • By curiosity, where is the algorithm code of the d3 path function to obtain the boundary profile of the map? I could not find it on the github repo.

Thanks.

Edit: I went through all the 44 projections in this example and I found glitches on the following 25:

Albers, Bromley, Collignon, Eckert II, Eckert IV, Eckert VI, Hammer, Hill, Goode Homolosine, Lambert cylindrical equal-area, Larrivée, Laskowski, McBryde–Thomas Flat-Polar Parabolic, McBryde–Thomas Flat-Polar Quartic, McBryde–Thomas Flat-Polar Sinusoidal, Mollweide, Natural Earth, Nell–Hammer, Polyconic, Sinu-Mollweide, van der Grinten, van der Grinten IV, Wagner IV, Wagner VII, Winkel Tripel.

like image 446
Daniel Ramos Avatar asked Jan 24 '17 15:01

Daniel Ramos


1 Answers

While I'm pretty sure you are using the projection.inverse function correctly, relying on :

if (λ > 180 || λ < -180 || φ > 90 || φ < -90) { i += 4; continue; }

to clip the projection will always fail as projection.inverse appears to always return angles within 180 degrees east/west. While there may be a way to modify the projection itself to return values greater than 180 degrees, it is likely more difficult than other approaches (And honestly, that goes well beyond any answer I can give).

Likewise, using an SVG path to represent the world outline and then using that as the basis to determine if a point should be drawn is probably complicating the matter.

Instead, assuming a circular disc, you can easily compute the radius of the disc and from there determine if a pixel should be drawn:

var edge = {};
var center = {};

edge.x = projection([180 - 1e-6, 0])[0];
edge.y = projection([180 - 1e-6, 0])[1];

center.x = width/2;
center.y = height/2;

var radius = Math.pow( Math.pow(center.x - edge.x,2) + Math.pow(center.y - edge.y,2) , 0.5 )

Using the radius of the disc, we can then compute if a pixel falls on the disc or beyond it in the for loop:

for (var y = 0, i = -1; y < height; ++y) {
    for (var x = 0; x < width; ++x) {

    var p = projection.invert([x, y]), λ = p[0], φ = p[1];

        if (Math.pow( Math.pow(center.x-x,2) + Math.pow(center.y-y,2), 0.5) < radius) {

            if ( λ > 180 || λ < -180 || φ > 90 || φ < -90 ) { i += 4; continue; }
                var q = ((90 - φ) / 180 * dy | 0) * dx + ((180 + λ) / 360 * dx | 0) << 2;
                targetData[++i] = sourceData[q];
                targetData[++i] = sourceData[++q];
                targetData[++i] = sourceData[++q];
                targetData[++i] = 255;
        }
        else {
            targetData[++i] = 0;
            targetData[++i] = 0;
            targetData[++i] = 0;
            targetData[++i] = 0;
        }
    }
}

Together, these gave me:

enter image description here

For aesthetic effect, it might be worthwhile trimming the radius down by a certain percent. Of course for different projections this approach may be either difficult or impossible.

I've put the code into a bl.ock here (I moved it to d3 v4 in the process, partly to see if the behavior of projection.inverse is the same).

For the third part of your question, you could try d3's graticule (graticule.outline) function for some info on how d3 gets the boundary profile of a projection.

like image 121
Andrew Reid Avatar answered Oct 28 '22 13:10

Andrew Reid