This post produces an interactive map that features analysis output mapped onto the geographical region to which it applies. This specific map answers the question: which areas in and around New Orleans exhibited the most growth (or loss) in active addresses over the past two years?
The interactive map itself, which this code produces, is accessible below in an iFrame.
1. Supply a shapefile that represents the geographical regions onto which the output should be mapped.
A zip code shapefile can be downloaded from https://www.census.gov/geo/maps-data/data/cbf/cbf_zcta.html. Technically, zip codes represent mail routes, not geographic areas, thus this file contains zip code tabulation areas (ZCTA).
library(maptools) # set directory to the folder where you saved the downloaded shapefile setwd("") zipcodes.shp <- readShapeSpatial("cb_2014_us_zcta510_500k.shp") # Re-set working directory to wherever you’d like to save the outputted GeoJSON and map setwd("") class(zipcodes.shp) names(zipcodes.shp) head(zipcodes.shp$ZCTA5CE10)
# Subset shapefile to select only those zip codes in the output file. # The output file I'm using was produced as part of the previous post. zips <- unique(output.2$zipcode) nola <- subset(zipcodes.shp, ZCTA5CE10 %in% zips) remove(zips)
2. Merge analysis output with shapefile and perform minor editing.
You might remember from the sequence of posts on transforming this data set for analysis that zip code 70068 appears twice – part of the zip code is in St. John the Baptist Parish, and part is in St. Charles Parish. These two parts of the zip code have different numbers of addresses, and different estimates for change over the past two years.
Unfortunately, it’s not possible to supply two different sets of information for the same zip code. There are many more 70068 addresses in the St. John the Baptist parish, so for convenience and simplicity here, I’m going to drop the 70068 observation in St. Charles Parish. (A more accurate, but also more lengthy, approach would be to return to the data wrangling stage, aggregate by zip code and date (thereby summing all 70068 addresses for each time point), then re-run the analysis to obtain only one estimate of change over time for 70068.)
# For convenience, drop the 70068 instance in St. Charles Parish. delete <- subset(output.2, zipcode == 70068) remove(delete) temp.data <- subset(output.2, unique.ID != 57) # Merge output.2 data into shapefile # The trend variable represents change in active addresses over time temp.data <- subset(temp.data, select= c("zipcode","parish","trend")) temp.data$ZCTA5CE10 <- temp.data$zipcode nola <- merge(nola, temp.data, by="ZCTA5CE10") remove(temp.data) names(nola) # Rename and clean up variables library(plyr) # renaming variables, ddply function nola <- rename(nola, c("ZCTA5CE10"="Zip_Code", "parish"="Parish")) nola$trend <- round(nola$trend, 2) # Create a character variable for trend that can be edited with text nola$Addresses_Added_Per_Day <- as.character(nola$trend) names(nola) print(nola$Addresses_Added_Per_Day) nola$Addresses_Added_Per_Day[nola$Addresses_Added_Per_Day=="0"] <- "0 (no change)"
3. Create the GeoJSON file.
Note that if you already have a GeoJSON file by the same name in the same location (for example, if you're editing a map and checking the edits), this command will throw an error. Deleting the old GeoJSON file will resolve the issue.
library(rgdal) writeOGR(nola, "nola_geojson", layer='nola', driver="GeoJSON")
4. Use Leaflet to render the map with choropleth colors and pop-ups.
# Set the demarcation points to associate trend values # with different colors on the choropleth graph cuts <- c(-0.7, -0.001, 0.0001, 0.25, 0.5, 0.75, 1.0, 1.25) # Fields to include in the popup popup <- c("Parish","Zip_Code","Addresses_Added_Per_Day") library(leafletR) # Create color gradient for choropleth based upon levels “trend” # I’ve set zip codes with negative change (which are losing addresses over time) # to blue, and zip codes with no change to white. sty <- styleGrad(prop="trend", breaks = cuts, # from the "cuts" creation, above right=FALSE, style.par="col", style.val=c("blue", "white",rev(heat.colors(5))), leg="Addresses Added / Day", lwd=1) # Create the map and load into browser map <- leaflet(data = 'nola_geojson', dest = "”,# specify your working directory, from above style=sty, title = "interactive_choropleth", base.map = "positron", incl.data = TRUE, popup = popup) remove(map, cuts, nola, popup, sty) # remove(zipcodes.shp)
This produces an .html file, saved to the folder you specified, which can be viewed as a stand-alone webpage or embedded in another page as an iFrame.