Jon Zelner Social epidemiology, infectious diseases, and statistics

Voronoi intensity mapping with R and d3: Part 2

In the last post, we focused on how to generate a voronoi diagram from a set of points in R, clip the voronoi diagram to a polygon - in this case the city boundaries of Philadelphia - and write all of this goodness out to geojson.

What we’re going to do now is to focus on how to get this information into D3 and to plot it like so:

Preprocessing

The first thing I did was to convert my voronoi polygons into topojson to make a smaller, faster-loading file using Mike Bostock’s topojson package:

topojson -o phl_voronoi.json phl_voronoi.json --id-property=id

Loading

So, the first housekeeping thing I did was to make a div for the map to hang off of:

<div id="map"></div>

Then, using a d3 selector, I go ahead and create an SVG element attached to the map div:

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

I used Mike Bostock’s queue function to pre-fetch the data before calling the function that will go ahead and draw the map. We’re going to load both a topojson full of the voronoi polygons we just created and the geojson from OpenDataPhilly with the configuration of Philadelphia neighborhoods:

    queue()
        .defer(d3.json, "/assets/phl_voronoi.json")
        .defer(d3.json, "/assets/Neighborhoods_Philadelphia.geojson")
        .await(ready);

Drawing

Now that we’ve gotten everything loaded up, it’s time to do what we came here to do: make the map. But before we do this, we have to fix the projection so that our map is scaled and centered properly. I took the approach suggested in this stackoverflow post, to start with some default values and then center and shift so that the map is centered and takes up the maximum space within the area provided:

    var json = topojson.feature(phila, phila.objects.phl_voronoi_2, function(a,b) { return a === b;});
    var center = d3.geo.centroid(json);


    var scale = 1050;
    var offset = [width/2, height/2];
    var projection = d3.geo.mercator()
        .scale(1)
        .translate([0,0]);


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

    var b = path.bounds(json);

    s = .95 / Math.max((b[1][0] - b[0][0]) / width, (b[1][1] - b[0][1]) / height),
    t = [(width - s * (b[1][0] + b[0][0])) / 2, (height - s * (b[1][1] + b[0][1])) / 2];

    // Update the projection to use computed scale & translate.
    projection
        .scale(s)
        .translate(t);

Ok, now we’ve gotten to the point where we’re going to actually do something that impacts our visualization. I’m a big fan of lodash for functional-style evaluation in javascript. So wherever you see something like _.forEach, we’re using a lodash function. What I did below was to iterate through each of the polygons in the FeatureCollection in the voronoi topojson and 1) project them so they’ll display properly on the map and more accurately represent physical space and 2) calculate their density by converting them into d3 Polygons and accessing the area property of each. Recall that we’re going to estimate the density of homicides within each voronoi cell as 1/area, so we end up with the following:

    density_vals = []
    features = topojson.feature(phila, phila.objects.phl_voronoi_2).features
    projected_coords = _.forEach(features, function(d) {

        p = []
        c = d.geometry.coordinates[0]
        for (i=0; i<c.length-1; i++) {
            p.push(projection(c[i]));
        }

        d.density = 1.0/Math.abs(d3.geom.polygon(p).area());
        density_vals.push(d.density);
    });

Note the call to Math.abs. This is because the polygon.area() function will return a negative value for area if vertices are listed in a counterclockwise fashion. So this ensures we get a positive value for area.

Now that we’ve gone ahead and attached the density values to each polygon, we can calculate the color scale to vary from white to red. There are decidedly more elegant ways to do this, but this is what I did:

    var color_domain = _(density_vals)
        .thru(function(x) {
            low_val = _.min(x);
            high_val = _.max(x);
            return [low_val, high_val];
        })
        .value()

    var color = d3.scale.log()
        .domain(color_domain)
        .range(["white","red"])

Finally, we’ll attach our voronoi polygons to the DOM:

    svg.selectAll("voronoi")
        .data(features)
        .enter()
        .append("path")
        .attr("d", path)
        .attr("fill", function(d) { return color(d.density);})
        .attr("stroke", function(d) { return color(d.density);})

Tooltips

One thing that I wanted to be able to do was to show the neighborhoods corresponding to a given outline to make it easy to understand where the most concentrated risk was. To do this, I used the d3 tip package, which makes it easy to make simple tooltips. You can add css, etc, but I went with a simple html one. So, first you make the tooltip object and tell it what to do with the data. In my case, I just want it to pop-up the name property. After we make this, we have to call it to attach it to the DOM and define some functions that will be called on mouseover and mouseout:

    tip = d3.tip().html(function(d) {
        return d.properties.listname;});

    svg.call(tip)

    function mouseover(d) {
        n = svg.select(" #"+d.properties.name);
        n.attr("class", "neighborhoods--show")
        tip.show(d)
    }

    function mouseout(d) {
        n = svg.select(" #"+d.properties.name);
        n.attr("class", "neighborhoods")
        tip.hide(d)
    }

Note that the mouseover and mouseout functions manipulate the class of the polygon. Changing the class to ‘neighborhood–show’ makes the outline red like so:

    .neighborhoods {
    fill: none;
    stroke: black;
    stroke-width: 0.2;
    pointer-events: all;
    }

    .neighborhoods--show {
    stroke: red;
    fill: none;
    pointer-events:all;
    stroke-width:2.0;

    }

An important thing to note here that I didn’t understand at first, is that ‘pointer-events:all’ lets us track mouseovers even when the polygon is unfilled. Otherwise, the mouseover function will only be called when the mouse is over the very thin borders of the neighborhood polygons.

And now that we’re all ready for them, we’ll go ahead and attach our neighborhood outline to the DOM:

    svg.selectAll("neighborhoods")
        .data(neigh.features)
        .enter()
        .append("path")
        .attr("d", path)
        .attr("id", function(d) { return d.properties.name;})
        .attr("class", "neighborhoods")
        .on("mouseover", mouseover)
        .on("mouseout", mouseout);

And that’s pretty much it for the Voronoi intensity map!

Conclusion

Again, this one is definitely a work-in-progress and by no means a perfect implementation. But it does the job pretty well and was not all that difficult to cobble together. In the future, a better implementation would probably pre-compute the density values and store them in the input topojson, or cache them server-side. For an SVG this complex, it might even make sense to pre-render the SVG element using a headless browser like PhantomJS and then just send that to the browser.

Another valid knock against this kind of map is that it doesn’t give a good sense of proportionality. We see areas with a high raw number of homicides on the map, but it’s not clear how these relate to some kind of baseline, whether it’s population density or some historical stats, like murder rates from the 1990s. This kind of spatial relative risk mapping is an interest of mine and something I hope to explore soon using these kinds of techniques.