Make and aggregate your own dot density maps in R

Last post, I talked about dot density maps as a way of taking polygons (like zip codes) and distributing their population evenly across the area in the form of dots. If 1000 people live in a zip code, then place 1000 random dots inside it. That's useful for visualizing population density, but what if you could then add these dots back together across new polygons? You'd be able to then get a pretty good estimate of the population within a new boundary.

For example, we can take the locations of clinics in Pretoria and create "Thessian" polygons (otherwise known as a Voronoi diagram) that show the area for which that particular clinic is closest. In other words, if you are in the Thessian polygon for Clinic A, then you are closer to Clinic A than to any other clinic. It's a rudimentary way of calculating a catchment area for each clinic. Then, we can overlay our population dot density map. Finally, we can just count the dots inside each clinic catchment to estimate the total population served by that clinic.

I'll include my R code below so you can replicate these results if you'd like. A sidenote on mapping in R: I never understood why people would map in R, when you could map in QGIS, Tilemill, Mapbox, Leaflet or other programs that seemed more intuitive. What I now recognize is that R keeps everything reproducible, which is essential for publishing. Also, it makes in incredibly easy to repeat your entire script with new data or variables. For example, we can rerun this same analysis for the entire country of South Africa with only a few extra lines of code!

Larger resolution image: here.

Below is the R code. It is thoroughly commented but email me or write a comment if something doesn't make sense and I can update the Gist. The code is a bit long, but that's partly due to some data cleanup in the beginning, plus additional code for printing all four of the maps in the GIF above.