Skip to content

Instantly share code, notes, and snippets.

@blehman

blehman/README.md

Created March 12, 2017 01:00
Download ZIP

Revisions

  1. blehman created this gist Mar 12, 2017.
    3 changes: 3 additions & 0 deletions README.md
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,3 @@
    ## Voronoi of the 50 largest urban areas in the US (d3.v4)

    Updating [Neil Freeman’s Block](http://bl.ocks.org/fitnr/9986268) to d3.v4 that was riffing on the [State Capitals Voronoi](http://jasondavies.com/maps/voronoi/us-capitals/) from Jason Davies, and using his [spherical voronoi](https://www.jasondavies.com/maps/voronoi/) script.
    297 changes: 297 additions & 0 deletions d3.geo.voronoi.js
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,297 @@
    (function() {

    var π = Math.PI,
    degrees = 180 / π,
    radians = π / 180,
    ε = 1e-15,
    circle = d3.geoCircle();

    geoVoronoi = function(points, triangles) {
    if (arguments.length < 2) triangles = geoDelaunay(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;
    })
    };
    };

    geoDelaunay = 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.geoInterpolate(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)));
    }

    })();
    93 changes: 93 additions & 0 deletions index.html
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,93 @@
    <!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.v4.min.js"></script>
    <script src="//d3js.org/d3-queue.v3.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.geoAlbers()
    .precision(0.1);

    var path = d3.geoPath()
    .projection(projection)
    .pointRadius(1.5);

    var svg = d3.selectAll("#map").append("svg")
    .attr("width", width)
    .attr("height", height);

    var fill = d3.scaleOrdinal(d3.schemeCategory10);

    d3.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 = geoVoronoi(coordinates, geoDelaunay(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>
    28 changes: 28 additions & 0 deletions style.css
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,28 @@
    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;
    }
    58 changes: 58 additions & 0 deletions top-50-metros.csv
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,58 @@
    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
    13,128 changes: 13,128 additions & 0 deletions us-10m.json
    View file
    13,128 additions, 0 deletions not shown because the diff is too large. Please use a local Git client to view these changes.