housing-price-maps-ggplot2.R

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")