Skip to content

Instantly share code, notes, and snippets.

@jgbos

jgbos/README.md

Last active December 16, 2015 03:18
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jgbos/5368756 to your computer and use it in GitHub Desktop.
Save jgbos/5368756 to your computer and use it in GitHub Desktop.
Download ZIP
Distances from North Korea
Raw
README.md

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.

Raw
index.html
<!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="http://d3js.org/d3.v3.min.js"></script>
<script src="http://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>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment