Skip to content

Instantly share code, notes, and snippets.

@fitnr
Last active June 10, 2017 05:41
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save fitnr/9986268 to your computer and use it in GitHub Desktop.
Save fitnr/9986268 to your computer and use it in GitHub Desktop.
US Urban Areas Voronoi
(function() {
var π = Math.PI,
degrees = 180 / π,
radians = π / 180,
ε = 1e-15,
circle = d3.geo.circle().angle(90);
d3.geo.voronoi = function(points, triangles) {
if (arguments.length < 2) triangles = d3.geo.delaunay(points);
if (!triangles) triangles = [];
var n = points.length;
var edgeByStart = [];
triangles.forEach(function(t) {
var edges = t.edges;
edges.reverse();
for (var i = 0, prev = edges[2]; i < 3; ++i) {
var e = edges[i];
e.prev = prev;
e.triangle = t;
edgeByStart[e.b] = (prev = e);
}
});
return {
type: "GeometryCollection",
geometries: n === 1 ? [{type: "Sphere"}]
: n === 2 ? hemispheres(points[0], points[1])
: points.map(function(_, i) {
var cell = [],
neighbors = [],
o = {type: "Polygon", coordinates: [cell], neighbors: neighbors},
e00 = edgeByStart[i],
e0 = e00,
e = e0,
centre0 = e.triangle.centre;
do {
var centre = e.triangle.centre;
if (dot(centre, centre0) < ε - 1) {
var a = cartesian(points[e0.b]), b = cartesian(points[e.a]),
c = normalise([a[0] + b[0], a[1] + b[1], a[2] + b[2]]);
if (dot(centre, cross(a, b)) > 0) c[0] = -c[0], c[1] = -c[1], c[2] = -c[2];
cell.push(spherical(c));
}
cell.push(spherical(centre));
neighbors.push(e.a);
centre0 = centre;
if (e === e00 && e0 !== e00) break;
e = (e0 = e).prev.neighbour;
} while (1);
return o;
})
};
};
d3.geo.delaunay = function(points) {
var p = points.map(cartesian),
n = points.length,
triangles = convexhull3d(p);
if (triangles.length) return triangles.forEach(function(t) {
t.centre = circumcentre(t);
// TODO reuse original points.
t[0] = spherical(t[0]);
t[1] = spherical(t[1]);
t[2] = spherical(t[2]);
}), triangles;
if (n > 2) {
var edgeByEnds = {},
N = normalise(cross(subtract(p[1], p[0]), subtract(p[2], p[0]))),
M = [-N[0], -N[1], -N[2]],
p0 = p.shift(),
reference = subtract(p[0], p0);
p.sort(function(a, b) { return angle(reference, subtract(a, p0)) - angle(reference, subtract(b, p0)); });
p.unshift(p0);
var d = [0, 1, 2];
for (var i = 2; i < n; ++i) {
d[0] = 0, d[1] = i - 1, d[2] = i;
var t = orient(p, d, N);
t.edges = [
edgeByEnds[d[0] + "," + d[1]] = new Edge(d[0], d[1]),
edgeByEnds[d[1] + "," + d[2]] = new Edge(d[1], d[2]),
edgeByEnds[d[2] + "," + d[0]] = new Edge(d[2], d[0])
];
t.centre = circumcentre(t);
triangles.push(t);
t = orient(p, d, M);
t.edges = [
edgeByEnds[d[0] + "," + d[1]] = new Edge(d[0], d[1]),
edgeByEnds[d[1] + "," + d[2]] = new Edge(d[1], d[2]),
edgeByEnds[d[2] + "," + d[0]] = new Edge(d[2], d[0])
];
t.centre = circumcentre(t);
triangles.push(t);
}
for (var k in edgeByEnds) {
var edge = edgeByEnds[k],
ends = k.split(",");
edge.neighbour = edgeByEnds[ends.reverse().join(",")];
}
return triangles;
}
};
function hemispheres(a, b) {
var c = d3.geo.interpolate(a, b)(.5),
n = cross(cross(cartesian(a), cartesian(b)), cartesian(c)),
m = 1 / norm(n);
n[0] *= m, n[1] *= m, n[2] *= m;
var ring = circle.origin(spherical(n))().coordinates[0];
return [
{type: "Polygon", coordinates: [ring]},
{type: "Polygon", coordinates: [ring.slice().reverse()]}
];
}
function convexhull3d(points) {
var n = points.length;
if (n < 4) return []; // coplanar points
var t = points.slice(0, 3);
t.n = cross(subtract(t[1], t[0]), subtract(t[2], t[0]));
for (var i3 = 3; i3 < n && coplanar(t, points[i3]); ++i3);
if (i3 === n) return []; // coplanar points
var triangles = [], edgeByEnds = {};
// First find a tetrahedron.
[[1, 2, i3, 0], [0, 2, 1, i3], [0, i3, 2, 1], [0, 1, i3, 2]].forEach(function(d) {
var t = orient(points, d, points[d[3]]); // also sets normal
t.edges = [
edgeByEnds[d[0] + "," + d[1]] = new Edge(d[0], d[1]),
edgeByEnds[d[1] + "," + d[2]] = new Edge(d[1], d[2]),
edgeByEnds[d[2] + "," + d[0]] = new Edge(d[2], d[0])
];
triangles.push(t);
});
for (var k in edgeByEnds) {
var edge = edgeByEnds[k],
ends = k.split(",");
edge.neighbour = edgeByEnds[ends.reverse().join(",")];
}
for (var i = 3; i < n; ++i) {
if (i === i3) continue;
var p = points[i];
var H = [];
for (var j = 0, m = triangles.length; j < m; ++j) {
var t = triangles[j];
if (!t || !visible(t, p)) continue;
// remove edges of t
var e = t.edges;
e[0].remove(H), e[1].remove(H), e[2].remove(H);
triangles[j] = null;
}
edgeByEnds = {};
for (var j = 0, m = H.length, e; j < m; ++j) {
if ((e = H[j]).removed) continue;
var t = [points[e.b], points[e.a], p];
t.n = cross(subtract(t[1], t[0]), subtract(t[2], t[0]));
t.edges = [
new Edge(e.b, e.a),
edgeByEnds[e.a + "," + i] = new Edge(e.a, i),
edgeByEnds[i + "," + e.b] = new Edge(i, e.b)
];
(e.neighbour = t.edges[0]).neighbour = e;
triangles.push(t);
}
for (var k in edgeByEnds) {
var edge = edgeByEnds[k],
ends = k.split(",");
edge.neighbour = edgeByEnds[ends.reverse().join(",")];
}
}
return triangles.filter(Boolean);
}
function circumcentre(t) {
var p0 = t[0],
p1 = t[1],
p2 = t[2],
n = t.n,
m2 = 1 / norm2(n),
m = Math.sqrt(m2),
radius = asin(.5 * m * norm(subtract(p0, p1)) * norm(subtract(p1, p2)) * norm(subtract(p2, p0))),
α = .5 * m2 * norm2(subtract(p1, p2)) * dot(subtract(p0, p1), subtract(p0, p2)),
β = .5 * m2 * norm2(subtract(p0, p2)) * dot(subtract(p1, p0), subtract(p1, p2)),
γ = .5 * m2 * norm2(subtract(p0, p1)) * dot(subtract(p2, p0), subtract(p2, p1)),
centre = [
α * p0[0] + β * p1[0] + γ * p2[0],
α * p0[1] + β * p1[1] + γ * p2[1],
α * p0[2] + β * p1[2] + γ * p2[2]
],
k = norm2(centre);
if (k > ε) centre[0] *= (k = 1 / Math.sqrt(k)), centre[1] *= k, centre[2] *= k;
else n[0] *= m, n[1] *= m, n[2] *= m, centre = n;
if (!visible(t, centre)) centre[0] *= -1, centre[1] *= -1, centre[2] *= -1, radius = π - radius, centre.negative = true;
centre.radius = radius;
return centre;
}
function norm2(p) { return dot(p, p); }
function norm(p) { return Math.sqrt(norm2(p)); }
function dot(a, b) {
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
function cross(a, b) {
return [a[1] * b[2] - a[2] * b[1],
a[2] * b[0] - a[0] * b[2],
a[0] * b[1] - a[1] * b[0]];
}
function subtract(a, b) {
return [a[0] - b[0], a[1] - b[1], a[2] - b[2]];
}
function planeDistance(t, p) {
return dot(t.n, p) - dot(t.n, t[0]);
}
function visible(t, p) {
return dot(t.n, p) - dot(t.n, t[0]) > ε;
}
function coplanar(t, p) {
return Math.abs(dot(t.n, p) - dot(t.n, t[0])) < ε;
}
function asin(x) {
return Math.asin(Math.max(-1, Math.min(1, x)));
}
function spherical(cartesian) {
return [
Math.atan2(cartesian[1], cartesian[0]) * degrees,
asin(cartesian[2]) * degrees
];
}
function cartesian(spherical) {
var λ = spherical[0] * radians,
φ = spherical[1] * radians,
cosφ = Math.cos(φ);
return [
cosφ * Math.cos(λ),
cosφ * Math.sin(λ),
Math.sin(φ)
];
}
// Construct triangle from first three points, with last point behind.
function orient(points, d, p) {
var triangle = [points[d[0]], points[d[1]], points[d[2]]];
triangle.n = cross(subtract(triangle[1], triangle[0]), subtract(triangle[2], triangle[0]));
if (visible(triangle, p)) {
triangle.reverse();
var t = d[0];
d[0] = d[2];
d[2] = t;
triangle.n[0] *= -1, triangle.n[1] *= -1, triangle.n[2] *= -1;
}
return triangle;
}
function Edge(a, b) {
this.a = a;
this.b = b;
this.neighbour = this.prev = this.triangle = null;
this.removed = false;
}
Edge.prototype.remove = function(horizon) {
var neighbour = this.neighbour;
if (neighbour) {
horizon.push(neighbour);
neighbour.neighbour = null;
this.neighbour = null;
}
this.removed = true;
};
function normalise(d) {
var m = 1 / norm(d);
d[0] *= m, d[1] *= m, d[2] *= m;
return d;
}
function angle(a, b) {
return Math.acos(dot(a, b) / (norm(a) * norm(b)));
}
})();
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<title>US Urban Area Voronoi</title>
<link rel="stylesheet" href="./style.css">
<script src="http://d3js.org/d3.v3.min.js"></script>
<script src="http://d3js.org/queue.v1.min.js"></script>
<script src="http://d3js.org/topojson.v1.min.js"></script>
<script src="./d3.geo.voronoi.js?1"></script>
</head>
<body>
<div id="map"></div>
<script>
var width = 960,
height = 500;
var projection = d3.geo.albers()
.precision(0.1);
var path = d3.geo.path()
.projection(projection)
.pointRadius(1.5);
var svg = d3.selectAll("#map").append("svg")
.attr("width", width)
.attr("height", height);
var fill = d3.scale.category10();
queue()
.defer(d3.json, "./us-10m.json")
.defer(d3.csv, "./top-50-metros.csv")
.await(ready);
function ready(error, us, capitals) {
if (error) {
console.log(error);
}
var defs = svg.append("defs");
defs.append("path")
.datum(topojson.feature(us, us.objects.land))
.attr("id", "land")
.attr("d", path);
defs.append("clipPath")
.attr("id", "clip")
.append("use")
.attr("xlink:href", "#land");
svg.append("use")
.attr("xlink:href", "#land")
.attr("class", "land");
var coordinates = capitals.map(function(d) { return [+d.longitude, +d.latitude]; });
var voronoi = d3.geo.voronoi(coordinates, d3.geo.delaunay(coordinates)).geometries;
var g = svg.append("g").attr("clip-path", "url(#clip)");
g.selectAll(".voronoi")
.data(voronoi)
.enter().append("path")
.attr("class", "voronoi")
.style("fill", function(d, i) { return fill(d.color = d3.max(d.neighbors, function(n) { return voronoi[n].color; }) + 1 | 0); })
.attr("d", path)
.append("title")
.text(function(_, i) {
var d = capitals[i];
return d.name + ", " + d.description;
});
g.append("path")
.datum(topojson.mesh(us, us.objects.states, function(a, b) { return a !== b; }))
.attr("class", "states")
.attr("d", path);
g.selectAll(".voronoi-border")
.data(voronoi.map(function(cell) {
return {type: "LineString", coordinates: cell.coordinates[0]};
}))
.enter().append("path")
.attr("class", "voronoi-border")
.attr("d", path);
svg.append("path")
.datum({type: "MultiPoint", coordinates: coordinates})
.attr("class", "points")
.attr("d", path);
}
</script>
</body>
</html>
body {
font-family: "Helvetica Neue", Helvetica, Arial, sans-serif;
width: 960px;
height: 500px;
position: relative;
}
canvas {
cursor: move;
}
path {
stroke-linejoin: round;
}
.voronoi-border {
stroke: #fff;
fill: none;
}
.points {
pointer-events: none;
}
.states {
stroke: #000;
fill: none;
}
name description latitude longitude
New York--Newark NY--NJ--CT Urbanized Area 40.718357 -73.970221
Los Angeles--Long Beach--Anaheim CA Urbanized Area 33.982668 -118.104332
Chicago IL--IN Urbanized Area 41.827126 -87.895427
Miami FL Urbanized Area 26.175616 -80.231428
Philadelphia PA--NJ--DE--MD Urbanized Area 39.973331 -75.298163
Dallas--Fort Worth--Arlington TX Urbanized Area 32.812727 -96.972001
Houston TX Urbanized Area 29.784308 -95.393531
Washington DC--VA--MD Urbanized Area 38.897222 -77.189759
Atlanta GA Urbanized Area 33.824102 -84.331858
Boston MA--NH--RI Urbanized Area 42.373132 -71.140708
Detroit MI Urbanized Area 42.489752 -83.227211
Phoenix--Mesa AZ Urbanized Area 33.494067 -111.970041
San Francisco--Oakland CA Urbanized Area 37.690191 -122.128542
Seattle WA Urbanized Area 47.468576 -122.274888
San Diego CA Urbanized Area 32.928728 -117.128799
Minneapolis--St. Paul MN--WI Urbanized Area 44.9783 -93.28051
Tampa--St. Petersburg FL Urbanized Area 28.006247 -82.513188
Denver--Aurora CO Urbanized Area 39.710774 -104.955096
Baltimore MD Urbanized Area 39.234099 -76.64776
St. Louis MO--IL Urbanized Area 38.630681 -90.341011
Riverside--San Bernardino CA Urbanized Area 33.985049 -117.341972
Las Vegas--Henderson NV Urbanized Area 36.134893 -115.167218
Portland OR--WA Urbanized Area 45.520404 -122.651087
Cleveland OH Urbanized Area 41.44323 -81.608373
San Antonio TX Urbanized Area 29.514979 -98.475643
Pittsburgh PA Urbanized Area 40.456928 -79.950944
Sacramento CA Urbanized Area 38.638578 -121.324669
San Jose CA Urbanized Area 37.329845 -121.945145
Cincinnati OH--KY--IN Urbanized Area 39.185505 -84.462043
Orlando FL Urbanized Area 28.584434 -81.411184
Kansas City MO--KS Urbanized Area 39.039792 -94.594897
Indianapolis IN Urbanized Area 39.811262 -86.145521
Virginia Beach VA Urbanized Area 36.908627 -76.309596
Austin TX Urbanized Area 30.358509 -97.761515
Milwaukee WI Urbanized Area 43.055674 -88.100524
Columbus OH Urbanized Area 40.021161 -82.968532
Charlotte NC--SC Urbanized Area 35.249654 -80.815826
Providence RI--MA Urbanized Area 41.776323 -71.397401
Jacksonville FL Urbanized Area 30.240149 -81.642722
Memphis TN--MS--AR Urbanized Area 35.095071 -89.917857
Salt Lake City--West Valley City UT Urbanized Area 40.638204 -111.926483
Nashville-Davidson TN Urbanized Area 36.128138 -86.681006
Louisville/Jefferson County KY--IN Urbanized Area 38.231406 -85.674529
Richmond VA Urbanized Area 37.477942 -77.487983
Buffalo NY Urbanized Area 42.92493 -78.823675
Hartford CT Urbanized Area 41.740122 -72.69883
Bridgeport--Stamford CT--NY Urbanized Area 41.135297 -73.534005
New Orleans LA Urbanized Area 29.954653 -90.143683
Raleigh NC Urbanized Area 35.765702 -78.66551
Oklahoma City OK Urbanized Area 35.504826 -97.524305
Tucson AZ Urbanized Area 32.254891 -110.946141
Urban Honolulu HI Urbanized Area 21.364556 -157.939756
El Paso TX--NM Urbanized Area 31.780239 -106.389031
Birmingham AL Urbanized Area 33.469107 -86.791974
Albuquerque NM Urbanized Area 35.141344 -106.62889
Omaha NE--IA Urbanized Area 41.234859 -96.034641
Dayton OH Urbanized Area 39.754825 -84.16816
Display the source blob
Display the rendered blob
Raw
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment