Overlaying population density in Google Maps

This example takes data from the US Census in the form of shape codes and population, massages it into a merged GeoJSON file, and then overlays it on to a Google Map. Specifically, I wanted to look at the population density of Portland, Oregon with granularity down to the block group.

This work is based heavily on Mike Bostock’s examples for California population density, particularly this one. I also recommend reading the “Let’s Make a Map” series as it covers most of the concepts used with much greater depth.

In order to get this in to Google Maps as an overlay I only made a few changes to Mike’s template. First, since I was only interested in the Portland metro area I reduced my county list to only Multnomah, Clackamas, and Washington counties. Second, since Google Maps would be handling projection for me I didn’t worry about applying a projection to my GeoJSON/TopoJSON data. Third, I played around with the density thresholds primarily through trial and error to improve contrast within the bulk of the city (with the drawback of saturating most rural areas in to the lowest bin). Fourth, I wanted to retain the block group definitions visually so I omitted the command to merge based on equal densities (at the expense of final file size). Fifth, I used different factors for simplifying the TopoJSON – again this was done through trial and error (note: I realized after the fact that the “-s” option is more appropriate for simplifying non-projected data rather than “-p”). As a word of warning, watch your GeoJSON file closely for size – if you are covering larger areas or more complex features loading times may become detrimental to the user experience unless you simplify the features further.

So, let’s see how these modifications play out in code. Below are the modified areas of the bash script from the California population density example linked above.

The first modification is in line 4 below, where instead of taking all items except for the header row we’ll grep for the lines of interest.

  1. if [ ! -f cb_${YEAR}_${STATE}_bg_B01003.ndjson ]; then
  2.   for COUNTY in $(ndjson-cat cb_${YEAR}_${STATE}_counties.json \
  3.       | ndjson-split \
  4.       | grep -E 'Multnomah|Clackamas|Washington' \
  5.       | ndjson-map 'd[2]' \
  6.       | cut -c 2-4); do
  7.     echo ${COUNTY}
  8.     if [ ! -f cb_${YEAR}_${STATE}_${COUNTY}_bg_B01003.json ]; then
  9.       curl -o cb_${YEAR}_${STATE}_${COUNTY}_bg_B01003.json \
  10.         "http://api.census.gov/data/${YEAR}/acs5?get=B01003_001E&for=block+group:*&in=state:${STATE}+county:${COUNTY}&key=${CENSUS_KEY}"
  11.     fi
  12.     ndjson-cat cb_${YEAR}_${STATE}_${COUNTY}_bg_B01003.json \
  13.       | ndjson-split \
  14.       | tail -n +2 \
  15.       >> cb_${YEAR}_${STATE}_bg_B01003.ndjson
  16.   done
  17. fi

Form the TopoJSON file, but I’ve taken out the steps to project and to merge on equal densities. I also modified the simplification factor (as mentioned above, using -s probably would have been better for toposimplify).

  1. geo2topo -n \
  2.   blockgroups=<(ndjson-join 'd.id' \
  3.     <(shp2json cb_${YEAR}_${STATE}_bg_500k.shp \
  4.       | ndjson-split 'd.features' \
  5.       | ndjson-map 'd.id = d.properties.GEOID.slice(2), d') \
  6.     <(ndjson-map < cb_${YEAR}_${STATE}_bg_B01003.ndjson '{id: d[2] + d[3] + d[4], B01003: +d[0]}') \
  7.     | ndjson-map -r d3=d3-array 'd[0].properties = {density: d3.bisect([500, 1000, 2000, 4000, 6000, 8000, 10000, 12000], (d[1].B01003 / d[0].properties.ALAND || 0) * 2589975.2356)}, d[0]') \
  8.   | topomerge -k 'd.id.slice(0, 3)' counties=blockgroups \
  9.   | topomerge --mesh -f 'a !== b' counties=counties \
  10.   | toposimplify -p 0.0005 -f \
  11.   > topo.json

And then extract out the feature collections.

  1. topo2geo \
  2.   < topo.json \
  3.   blockgroups=blockgroups.json \
  4.   counties=counties.json

The county features weren’t of much importance for me, so I only continued on with the blockgroups from here. Getting the blockgroups in to Google Maps and modifying the styles is a pretty easy exercise. First, load the file in to a data layer and second, style to your preference. We quantized our population densities previously based on 8 threshold values so we need 9 colors to code these values. To keep this light weight you can define an array of colors directly (ColorBrewer.com can generate these arrays for you automatically), but since I knew I would be using d3 to create a nice legend anyway I went with one of the built-in scales.

Initialize the map, load and style the json data

  1. function initMap() {
  2. 	map = new google.maps.Map(document.getElementById('map'), {
  3. 		zoom: 11,
  4. 		center: {
  5. 			lat: 45.490,
  6. 			lng: -122.715
  7. 		}
  8. 	});
  9. 	map.data.loadGeoJson(
  10. 		'blockgroups.json');
  11.  
  12. 	map.data.setStyle(function(feature) {
  13. 		var density = feature.getProperty('density');
  14. 		var colors = d3.schemeOrRd[9]
  15. 		return {
  16. 			fillColor: colors[density],
  17. 			fillOpacity: 0.7,
  18. 			strokeWeight: 1
  19. 		};
  20. 	});
  21. 	...
  22. };

Now we have a map loaded with the color-coded blockgroups. The last thing that is missing is a legend to convey what the different colors actually mean. Again, I borrowed from Mike Bostock’s California population density example. Main adjustments were to modify the values to match the thresholds I previously set, and adjust positional parameters such that the SVG created would be small enough to embed as a legend in to the Google Map (only partially successful – I’ve found that I clip the legend if viewing on mobile in portrait mode).

  1. <div id="legend">
  2. 	<svg width="400" height="38" id="testsvg"></svg>
  3. </div>
  4. <script>
  5. 	var svg = d3.select("svg"),
  6. 		width = +svg.attr("width"),
  7. 		height = +svg.attr("height");
  8.  
  9. 	var color = d3.scaleThreshold()
  10. 		.domain([500, 1000, 2000, 4000, 6000, 8000, 10000, 12000])
  11. 		.range(d3.schemeOrRd[9]);
  12.  
  13. 	var x = d3.scaleSqrt()
  14. 		.domain([300, 14000])
  15. 		.rangeRound([0, 400]);
  16.  
  17. 	var g = svg.append("g")
  18. 		.attr("class", "key")
  19. 		.attr("transform", "translate(0,13)");
  20.  
  21. 	g.selectAll("rect")
  22. 		.data(color.range().map(function(d) {
  23. 			d = color.invertExtent(d);
  24. 			if (d[0] == null) d[0] = x.domain()[0];
  25. 			if (d[1] == null) d[1] = x.domain()[1];
  26. 			return d;
  27. 		}))
  28. 		.enter().append("rect")
  29. 		.attr("height", 8)
  30. 		.attr("x", function(d) {
  31. 			return x(d[0]);
  32. 		})
  33. 		.attr("width", function(d) {
  34. 			return x(d[1]) - x(d[0]);
  35. 		})
  36. 		.attr("fill", function(d) {
  37. 			return color(d[0]);
  38. 		});
  39.  
  40. 	g.append("text")
  41. 		.attr("class", "caption")
  42. 		.attr("x", x.range()[0])
  43. 		.attr("y", -6)
  44. 		.attr("fill", "#000")
  45. 		.attr("text-anchor", "start")
  46. 		.attr("font-weight", "bold")
  47. 		.text("Population per square mile");
  48.  
  49. 	g.call(d3.axisBottom(x)
  50. 			.tickSize(13)
  51. 			.tickValues(color.domain()))
  52. 		.select(".domain")
  53. 		.remove();
  54. </script>

The SVG gets wrapped in to a div element, which is then provided to the Google Maps API as a new control in the initMap() function.

  1. 	var legend = document.getElementById('legend');
  2. 	map.controls[google.maps.ControlPosition.BOTTOM_RIGHT].push(legend);

See the end result here.

Leave a Reply

Your email address will not be published. Required fields are marked *