I’ve been working on a project on state-level farmland values and appreciation rates, and colleagues have expressed interest in knowing how I achieved some of the elementary spatial data exploration for it, particularly the eye-catching animated choropleth maps that I built. You can see the code and results below!

Setup

library(expm)
library(extrafont)
library(httr)
library(jsonlite)
library(lubridate)
library(magick)
library(rgdal)
library(spdep)
library(spdplyr)
library(tidyverse)

loadfonts(quiet=TRUE)

Data

Data from St. Louis FRED

As I’m analyzing real farmland values and appreciation rates, I’ll need data on the general price level in the US over time. I choose the GDP deflator as it is a good broad indicator of the price level independent of particular consumer goods baskets.

fred.url <- paste(
  'https://api.stlouisfed.org/fred/series/observations',
  '?series_id=GDPDEF',
  '&api_key=', api.token.fred,
  '&file_type=json',
  sep=''
)

fred.return <- GET(fred.url)

fred.data <- 
  fromJSON(content(fred.return, 'text'))$observations %>%
  type_convert() %>%
  mutate(year = year(date)) %>%
  group_by(year) %>%
  summarize(gdpdef = mean(value))

write_csv(fred.data, '../data/fred-data.csv')

Data from USDA NASS

Data on farmland value is available directly from the USDA National Agricultural Statistics Service via their API:

# NASS Data
nass.url <- paste(
  'http://quickstats.nass.usda.gov/api/api_GET/?key=', api.token.nass,
  '&source_desc=SURVEY',
  '&reference_period_desc=YEAR',
  '&agg_level_desc=STATE',
  '&commodity_desc=AG LAND',
  '&statisticcat=ASSET VALUE',
  '&short_desc=AG LAND, INCL BUILDINGS - ASSET VALUE, MEASURED IN $ / ACRE',
  '&domain_desc=TOTAL',
  sep=''
)

nass.return <- GET(nass.url)

nass.data <-
  fromJSON(content(nass.return, 'text'))$data %>%
  as_tibble() %>%
  type_convert() %>%
  mutate(year = parse_integer(year)) %>%
  rename(landvalac = Value)

write_csv(nass.data, '../data/nass-data.csv')

Shapefile from US Census Bureau

A shapefile of the United States is publicly available from the US Census Bureau:

download.file(
  'http://www2.census.gov/geo/tiger/GENZ2016/shp/cb_2016_us_state_20m.zip',
  destfile = '../data/cb_2016_us_state_20m.zip')
unzip('../data/cb_2016_us_state_20m.zip', exdir = '../data/')

Manipulation

With the data downloaded from their respective sources, I can prepare it for my map generating procedures and spatial regressions. First, I can access the GDP deflator for 2010 (our base year for real values). Then generating real values and appreciation rates is fairly straightfoward.

For the shapefile, there are some districts that should be removed if they are not included in the NASS data. Then for the spatial regressions later on, I make sure that the states are alphabetized in order to ensure the spatial weights matrices I generate are in accordance with the order of the NASS data.

# Set a range of years
year.range <- 1960:2016

# Read in FRED data from local copy
fred.data  <- read_csv('../data/fred-data.csv')
gdpdef2010 <-
  fred.data %>%
  filter(year == 2010) %>%
  select(gdpdef) %>%
  as_vector()

# Function to calculate appreciation rates
pct <- function(x) 100 * (x / lag(x) - 1)
  
# Read in NASS data from local copy, adjust land values to real (2010) values,
# and add appreciation rates
nass.data <- read_csv('../data/nass-data.csv') %>%
  select(state_alpha, year, landvalac) %>%
  arrange(state_alpha, year) %>%
  inner_join(fred.data, by = 'year') %>%
  mutate(realvalac = landvalac * gdpdef / gdpdef2010) %>%
  group_by(state_alpha) %>%
  mutate(realvalac.PD1 = pct(realvalac))

# Read in US state shapefile
states <-
  readOGR(dsn='../data', layer='cb_2016_us_state_20m',
          stringsAsFactors = FALSE) %>%
  filter(!(STUSPS %in% c('DC', 'AK', 'HI', 'PR'))) %>%
  arrange(STUSPS)
## OGR data source with driver: ESRI Shapefile 
## Source: "../data", layer: "cb_2016_us_state_20m"
## with 52 features
## It has 9 fields
## Integer64 fields read as strings:  ALAND AWATER
states.tbl <- fortify(states, region='STUSPS')

Choropleth Maps

Farmland Values

Legend Preparation

A problem with an animated map is that I’m making 57 separate maps (one for each year) and stitching them together. Because of that, I need to explicitly tell ggplot2 that I want the legend to use the data available for all years in the animation, so that the colors have constant meaning year by year. There’s no built-in color scale in ggplot2 that accounts for the legend breaks that I want to hold constant across year, so I will use the scale_fill_manual function. This option, however, requires us to explicitly state my desired breakpoints and the colors I will use between them, so I define both below:

# Create bins and color palette for legend
# Use K-State Purple and an accent tan from branding guide
quantile(nass.data$realvalac, seq(0, 1, 0.1))
##           0%          10%          20%          30%          40% 
##     1.527728    11.524702    22.362028    45.977944   118.194399 
##          50%          60%          70%          80%          90% 
##   300.622549   552.699499   932.764775  1698.832913  3326.439259 
##         100% 
## 16473.475800
realvalac.levels <- c(0,   25,   50,   125,  275,   450,
                      750, 1150, 2000, 3500, 17000)
realvalac.labels <- c('0 - 25',    '25 - 50',   '50 - 125',   '125 - 275',
                      '275 - 450', '450 - 750', '750 - 1150', '1150 - 2000',
                      '2000 - 3500', '3500 - 17000')
colfunc <- colorRampPalette(c('#ffffff', '#cfa35c','#512888'))
realvalac.colors <- colfunc(10)
names(realvalac.colors) <- realvalac.labels
nass.data$realvalac.fc <- factor(
  cut(nass.data$realvalac, realvalac.levels),
  labels = realvalac.labels
)

Map Generation

Now that the legend breaks and colors have been created, I can generate my year-by-year choropleth maps using ggplot2:

# Generate a choropleth map for each year
maplist <- lapply(year.range, function(yr) {
  ggmap <- ggplot(nass.data %>% filter(year == yr)) +
    geom_map(aes(map_id=state_alpha, fill=realvalac.fc),
             map=states.tbl, color='#cccccc', size=0.20) +
    coord_map() +
    expand_limits(x = states.tbl$long, y = states.tbl$lat) +
    labs(title=paste('Real Farmland Value (', yr, ')', sep=''),
         subtitle='Including Buildings, Constant 2010 US Dollars/Acre',
         caption='Data: USDA National Agricultural Statistics Service') +
    scale_fill_manual(element_blank(),
                      values = realvalac.colors, drop=FALSE) +
    theme_void() +
    theme(text=element_text(family='Myriad Pro'),
          plot.title=element_text(face='bold'),
          plot.subtitle=element_text(size=8),
          plot.caption=element_text(size=6),
          legend.margin=margin(t=-0.6, r=0, b=-0, l=0, unit='cm'))
  
  ggsave(paste('../png/landvalues-levels-01-choropleth-', yr, '.png', sep=''),
         width = 6, height = 3.1, units = 'in', dpi=150, pointsize=14)
  
  return(ggmap)
  })

Animated Map Compilation

With a choropleth map generated and stored for each year in my sample, I can use ImageMagick to stitch them together into an animated GIF.

# Read the individual maps into a data structure for use with 'magick'
imglayers <- sapply(year.range, function(yr) {
  image_read(paste('../png/landvalues-levels-01-choropleth-', yr, '.png', sep=''))
})

# Generate an animated GIF with the individual maps and write to a file
imganim <- image_animate(image_join(imglayers), fps = 1, dispose = "previous")
image_write(imganim, '../gif/landvalues-levels-01-choropleth-anim.gif')