Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
library(readxl)
library(tigris)
library(sp)
library(maptools)
library(ggplot2)
library(rgdal)
library(ggthemes)
library(viridis)
library(extrafont)
library(rgeos)
library(raster)
library(dplyr)
library(ggsn)
# Process tabular data as in WaPo article
# Get the data from http://www.fhfa.gov/DataTools/Downloads/Documents/HPI/HPI_AT_ZIP5.xlsx
orig <- read_excel('HPI_AT_ZIP5.xlsx', skip = 6)
averages <- orig %>%
select(zip = `Five-Digit ZIP Code`, year = Year, change = `Annual Change (%)`) %>%
filter(year > 1989, change != '.') %>%
mutate(change = as.numeric(change)) %>%
group_by(zip) %>%
summarize(avg = mean(change, na.rm = TRUE)) %>%
mutate(avgf = cut(avg, breaks = c(-100, 0:5, 100), right = FALSE,
labels = c('Decrease', '0', '1', '2', '3', '4', '+5%')))
# Get the spatial data with tigris and sp
zips <- zctas(cb = TRUE)
ctys <- counties('TX', cb = TRUE)
dfw_metro <- ctys[ctys$NAME %in% c('Dallas', 'Tarrant', 'Collin', 'Denton'), ]
over_zips <- bind_rows(over(dfw_metro, zips, returnList = TRUE))
dfw_zips <- spTransform(zips[zips$ZCTA5CE10 %in% over_zips$ZCTA5CE10, ],
CRS("+init=epsg:26914"))
pri <- spTransform(primary_roads(), CRS("+init=epsg:26914"))
cities <- places('TX', cb = TRUE)
dfw_cities <- spTransform(cities[cities$NAME %in% c('Fort Worth', 'Dallas', 'Denton', 'Plano'), ],
CRS("+init=epsg:26914"))
dfw_cities$long <- coordinates(dfw_cities)[,1]
dfw_cities$lat <- coordinates(dfw_cities)[,2]
# Subset tabular data
dfw_averages <- averages[averages$zip %in% dfw_zips$ZCTA5CE10, ]
# Create a circle
radius <- ( (bbox(dfw_zips)[3] - bbox(dfw_zips)[1]) / 2 )
circle <- gBuffer(gCentroid(dfw_zips), width = radius, quadsegs = 500)
# Restrict the zips to the extent of the circle
dfw_clipped <- gIntersection(dfw_zips, circle, byid = TRUE, id = dfw_zips$ZCTA5CE10)
dfw_clipped$id <- row.names(dfw_clipped)
plot(circle)
plot(dfw_clipped, add = TRUE)
# Clip function for the roads thanks to Robin Lovelace (faster)
gClip <- function(shp, bb){
if(class(bb) == "matrix") b_poly <- as(extent(as.vector(t(bb))), "SpatialPolygons")
else b_poly <- as(extent(bb), "SpatialPolygons")
proj4string(b_poly) <- proj4string(shp)
gIntersection(shp, b_poly, byid = T)
}
pri_clipped <- gClip(pri, circle)
# Convert to SpatialLinesDataFrame for ggplot2
pri_clipped$id <- 1:length(pri_clipped)
# Build ggplot2 visualization
dfwf <- fortify(dfw_clipped, region = 'id')
circlef <- fortify(circle)
roadsf <- fortify(pri_clipped)
ggplot() +
geom_polygon(data = circlef, aes(x = long, y = lat, group = group),
fill = '#eaeaea') +
geom_map(data = dfw_averages, map = dfwf,
aes(fill = avgf, map_id = zip)) +
geom_path(data = roadsf, aes(x = long, y = lat, group = group), color = 'white',
size = 0.5) +
geom_text(data = dfw_cities@data, aes(x = long, y = lat, label = NAME),
color = 'black', fontface = 'bold', family = 'Tahoma') +
theme_map(base_family = 'Tahoma', base_size = 14) +
coord_equal() +
theme(legend.position = 'top',
legend.key = element_blank()) +
scale_fill_manual(values = viridis(7)[2:7],
guide = guide_legend(nrow = 1, direction = 'horizontal',
label.hjust = 0, label.position = 'bottom',
keywidth = 5.51, keyheight = 0.75, title = "")) +
labs(title = "Home prices in Dallas-Fort Worth, Tex.",
subtitle = "Annual average home-price change from 1990 to 2015",
caption = "Map by @kyle_e_walker | Source: Federal Housing Finance Agency") +
scalebar(data = circlef, dist = 20, location = "bottomleft")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment