D3 Examples - Distances from North Korea

Distances from North Korea

Jason Davies recently demostrated the power of map projections for visualizating the distance a missile from North Korea would have to travel to hit any point on earth. An assumption in his visualization is that the earth is spherical and not rotating during the trajectory. This visualization examines the distances given the following assumptions:

Spherical, non-rotating: Left plot copies Jason's creation , each contour is 1,000 km

Spherical, rotating: Middle plot assumes a rotating earth, but keeps the earth as a sphere

Elliptical, rotating: Right plot models the WGS-84 elliptical earth that rotates

As expected, the rotating earth makes a difference on the actual distance a missile would travel over the earth (keeping the range of the missile constant). Trajectories moving eastward along the direction the earth is rotating travel less distance and opposite for westward trajectories. There isn't much difference between spherical and elliptical.

Some notes:

Flight path angle is 45 degrees

In both visualations (Jason's and mine) the effects of atmospheric drag and the need to boost a target to the desired velocity are ignored

Crossing my fingers this is correct, there may be some integrated effects between flight time and distance traveled that I may not be taking into account. But the impact should be relatively small, at least it shouldn't change the projection to much.

I may try to look into producing a projection that preserves flight time as suggested by Jason. Since I can clearly model a spherical earth, the math can be simplified.

xxxxxxxxxx
 
<!DOCTYPE html>
<meta charset="utf-8">
<title>Distances from North Korea</title>
<style>
​
body {
  background: #fcfcfc;
}
​
.background {
  fill: #fff;
}
​
.outline {
  stroke: #000;
  stroke-width: 1px;
  fill: none;
}
​
.land {
  fill: #fff;
}
​
.land-glow {
  fill: #000;
  fill-opacity: .5;
  filter: url(#glow);
}
​
.border {
  stroke: #000;
  stroke-width: .5px;
  fill: none;
}
​
div {
    float: left;
}
​
.circle {
  stroke: #f00;
  stroke-width: .5px;
  fill: none;
}
​
.major {
  stroke-width: 1.5px;
  stroke-dasharray: none;
}
​
.graticule {
  stroke: #99f;
  stroke-width: .5px;
  stroke-opacity: .3;
  fill: none;
}
​
.caption {
  font-style: italic;
}
​
p {
    position: absolute;
    top: 50px;
    color: #000;
    font-family: Arial, Helvetica, sans-serif;
    font-weight: bold;
    font-size: 14px;
    text-align: center;
    width: 300px;
}
</style>
<body>
<div id="map1"><p> Spherical, Non-Rotating</div>
<div id="map2"><p> Spherical, Rotating</div>
<div id="map3"><p> Elliptical, Non-Rotating</div>
<script src="https://d3js.org/d3.v3.min.js"></script>
<script src="https://d3js.org/topojson.v1.min.js"></script>
<script>
​
var width = 300,
    height = 500;
​
var origin = [127.2565, 40.2012],
    degrees = 180 / Math.PI,
    radians = Math.PI / 180;
​
var projection = d3.geo.azimuthalEquidistant()
    .translate([150,250])
    .clipAngle(120)
    .scale(70)
    .rotate([-origin[0], -origin[1], 0])
    .precision(.1);
​
var angle = projection.clipAngle(),
    t = projection.translate();
​
var path = d3.geo.path().pointRadius(2.5).projection(projection),
    circle = d3.geo.circle().origin(origin),
    format = d3.format(",d");
​
var ellipse = function (data,φ,params){
    data.coordinates[0] = data.coordinates[0].map(function (d){
        var tf = flightTime(origin,φ*radians,params),
            xeci = toCartesian(d,params),
            xecef = rotate(xeci, tf,params); // rotates earth
​
        return toLLA(xecef,params);
    });
​
    return data;
}
​
var svg = d3.selectAll("#map1, #map2, #map3")
    .data(["spherical inertial", "spherical", "wgs84"])
    .append("svg");
​
svg.each(function (type){
    d3.select(this)
        .attr("width", width)
        .attr("height", height);
})
​
​
svg.append("filter")
    .attr("id", "glow")
    .append("feGaussianBlur")
    .attr("stdDeviation", 3);
​
svg.append("path")
    .datum({type: "Sphere"})
    .attr("class", "background");
​
var g = svg.append("g").attr("class","map");
​
svg.append("path")
    .attr("class", "graticule")
    .datum(d3.geo.graticule()());
​
svg.each(function (type){
    var g = d3.select(this),
        p = earth(type);
​
    var δ = 1e6 / p.Re * degrees;
    g.selectAll("path.circle")
        .data(d3.range(δ, angle, δ))
        .enter().append("path")
        .datum(function(d) { return ellipse(circle.angle(d)(),d,p); })
        .attr("class", function(_, i) { return (i + 1) % 5 ? "circle" : "major circle"; });
​
    g.selectAll("text")
        .data(projection.clipAngle() > 90 ? [5000, 10000] : [5000])
        .enter().append("text")
        .attr("dy", "-.35em")
        .append("textPath")
        .attr("xlink:href", function(_, i) { return "#text-" + angle + "-" + i; })
        .text(function(d) { return format(d) + "km"; });
});
​
svg.append("path")
    .datum({type: "Sphere"})
    .attr("class", "outline");
​
svg.append("path")
    .datum({type: "Point", coordinates: origin});
​
svg.each(redraw);
​
d3.json("/d/4090846/world-50m.json", function(error, world) {
    var land = topojson.feature(world, world.objects.land);
    g.append("path")
        .datum(land)
        .attr("class", "land-glow");
    g.append("path")
        .datum(land)
        .attr("class", "land");
    g.append("path")
        .datum(topojson.mesh(world, world.objects.countries))
        .attr("class", "border");
    g.each(redraw);
});
​
function redraw(type) {
    d3.select(this).selectAll("path").attr("d", path);
}
​
function earth(earth_type){
    var param = {
        // Semi-major axis of earth radus (equatorial radius) (m)
        Re: 6378137.0,
​
        // Earth gravitational constant (m^3/s^2)
        // Mass of Earth's Atmosphere not included
        mu: 3.986004418e14,
​
        // Earth rotation rate (rad/s)
        rotation_rate : 7292115.1467e-11,
​
        // Eccentricity (unitless)
        eccentricity: 0.08181919084262149
    }
​
    type = earth_type.split(" ");
    type.forEach(function (t){
        switch (t){
            case "wgs84":
                // Do nothing
                break;
            case "spherical":
                param.Re = 6371000; // Mean radius
                param.eccentricity = 0;
                break;
            case "inertial":
                param.rotation_rate = 0;
                break;
        }
    });
​
    return param;
}
​
function rotate(pos,time,p){
    var ω = p.rotation_rate,
        ct = Math.cos(time*ω),
        st = Math.sin(time*ω);
​
    return [
        ct * pos[0] + st*pos[1],
        -st * pos[0] + ct*pos[1],
        pos[2]
    ];
}
​
function flightTime(launch, φ, p){
    // Earth Parameters
    var μ = p.mu,
        Re = p.Re,
        γ = 45 * Math.PI/180;
​
    // Calculate earth center vector length at latitude, longitude
    var R0 = toCartesian(launch, p),
        r0 = Math.sqrt(R0.reduce(function (a,b){return a+b*b;}));
​
    // Velocity of Target
    var V = Math.sqrt(μ*(1-Math.cos(φ))/ (r0*Math.cos(γ)*(r0*Math.cos(γ)/Re - Math.cos(φ+γ))));
​
    // Constant to calculate flight time
    var λ = r0*V*V/μ;
​
    // Calculate Flight Time
    var cg = Math.cos(γ), 
        sg = Math.sin(γ), 
        tg = Math.tan(γ),
        cp = Math.cos(φ),
        sp = Math.sin(φ),
        cgp = Math.cos(γ + φ),
        cot = 1 / Math.tan(φ/2);
​
    var first_term = (tg*(1-cp) + (1-λ)*sp) / ((2-λ)*((1-cp)/(λ*cg*cg) + cgp/cg)),
        sec_term = 2*cg / (λ*Math.pow(2/λ-1,1.5)) * Math.atan2(Math.sqrt(2/λ-1),(cg*cot-sg));
​
    return r0 / (V*cg) * (first_term + sec_term);
}
​
function toCartesian(gc,p) {
    var ε = p.eccentricity,
        Re = p.Re;
​
    var lat = gc[1]*radians,
        lon = gc[0]*radians,
        alt = 0,
        slat = Math.sin(lat),
        clat = Math.cos(lat);
​
    var R = Re / Math.sqrt(1 - ε*ε*slat*slat);
    var x = (R + alt)*clat*Math.cos(lon);
    var y = (R + alt)*clat*Math.sin(lon);
    var z = (R + alt - ε*ε*R)*slat;
​
    return [x,y,z];
}
​
function toLLA(position,p){
    var ε = p.eccentricity,
        Re = p.Re;
​
    var x = position[0],
        y = position[1],
        z = position[2],
        xSq=x*x,
        ySq=y*y;
​
    // WGS84 ellipsoid constants:
    var eSq = ε*ε,
        r = Math.sqrt(xSq + ySq),
        b = Re*Math.sqrt(1-eSq);
​
    // Longitude
    var lon = Math.atan2(y, x);
​
    if (ε == 0){
        lat = Math.atan2(z,r);    
    }else {
        var l = eSq/2.0,
            lSq = l* l,
            m = Math.pow(r/Re, 2),
            n = Math.pow((1-eSq)*z/b, 2),
            i = -(2.0*lSq+m+n)/2.0,
            k = lSq*(lSq-m-n),
            q = Math.pow((m+n-4.0*lSq), 3)/216.0 + m*n*lSq,
            D = Math.sqrt((2.0*q-m*n*lSq)*m*n*lSq),
            beta = i/3.0 - Math.pow((q+D), (1.0/3.0)) - Math.pow((q-D), (1.0/3.0));
​
        var sign=0;
        if((m-n) > 0){
            sign = 1;
        }else if((m-n) < 0){
            sign = -1;
        }
​
        var t = Math.sqrt( Math.sqrt(beta*beta - k)
            - (beta+i)/2.0)
            - sign*Math.sqrt((beta-i)/2.0);
​
        var r0 = r/(t+l);
        var z0 = (1-eSq)*z/(t-l);
​
        // Latitude
        var lat = Math.atan(z0/((1-eSq)*r0));
    }
​
    return [lon*degrees, lat*degrees];
}
​
d3.select(self.frameElement).style("height", height + "px");
​
</script>
​
    