Chloropleth map example

Preparation

Load packages and set options

Setup the default knitr options for the R code sections below.

Load the packages used in this example

library(rgdal)    # for reading the shape file
library(maptools) # dependency 
library(viridis)  # for a nice palette
library(foreign)  # for reading the .dbf file
library(ggplot2)  # for plotting
library(scales)   # for the custom gradient scale
library(magrittr) # for piping data
library(dplyr)    # for filtering data
library(ggrepel)  # for non-overlapping labels
library(rgeos)    # for coordinate conversion
library(geosphere)# for the centroid function

Hawaii ahupuaa

Read the data

Data was downloaded from http://planning.hawaii.gov/gis/download-gis-data/ section 016 - CULTURAL AND DEMOGRAPHIC and the .zip archives extracted in the subdirectories ahupuaa and zipcodes

shape <- readOGR("ahupuaa/ahupuaa.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "ahupuaa/ahupuaa.shp", layer: "ahupuaa"
## with 725 features
## It has 18 fields

shape <- spTransform(shape, CRS("+proj=longlat"))
shape1 <- fortify(shape)
dbf <- read.dbf('ahupuaa/ahupuaa.dbf')

Plot the acreage of each ahupua

With the viridis scale

ggplot() + 
  geom_map(data = dbf, 
           aes(map_id = OBJECTID, 
               fill = as.numeric(gisacres)), 
           map = shape1) + 
  expand_limits(x = shape1$long, y = shape1$lat) + 
  guides(fill='none') + 
  scale_fill_viridis() + 
  theme(axis.title = element_blank(), 
        axis.text = element_blank(), 
        axis.ticks = element_blank())

With a custom gradient scale

ggplot() + 
  geom_map(data = dbf, 
           aes(map_id = OBJECTID, 
               fill = as.numeric(gisacres)), 
           map = shape1) + 
  expand_limits(x = shape1$long, y = shape1$lat) + 
  guides(fill='none') + 
  scale_fill_gradient2(
    limits = c(0, max(as.numeric(dbf$gisacres))), 
    low = muted("blue"), 
    midpoint = median(as.numeric(dbf$gisacres)), 
    high = muted("red")) + 
  theme(axis.title = element_blank(), 
        axis.text = element_blank(), 
        axis.ticks = element_blank())

Select just the Oʻahu ahupuaa

shape2 <- fortify(shape, region='mokupuni') %>%
  filter(id == 'Oʻahu')
dbf2 <- dbf %>% 
  filter( mokupuni == 'Oʻahu')
oahu_limits <- shape2 %>% 
  summarise(min_lat=min(lat), max_lat=max(lat),
            min_long=min(long), max_long=max(long))
oahu_ahupuaa <- fortify(shape, region='ahupuaa') %>% 
  filter(between(lat,  oahu_limits$min_lat,  oahu_limits$max_lat) & 
         between(long, oahu_limits$min_long, oahu_limits$max_long)) %>%
  group_by(id) %>%
  do(clong = centroid(cbind(.$long, .$lat))[1],
      clat = centroid(cbind(.$long, .$lat))[2]) %>%
  mutate(clong = unlist(clong), clat = unlist(clat))
ggplot(data = dbf2) + 
  geom_map(aes(map_id = OBJECTID, 
               fill = as.numeric(gisacres)), 
           map = shape1) + 
  expand_limits(x = shape2$long, y = shape2$lat) + 
  guides(fill='none') + 
  geom_text_repel(data = oahu_ahupuaa, 
                  aes(x = clong, y = clat, 
                      label = id, size = 0.05, 
                      vjust='middle', hjust='center')) +
  guides(size = 'none') +
  scale_fill_gradient2(
    limits = c(0, max(as.numeric(dbf2$gisacres))), 
    low = muted("red"), 
    midpoint = median(as.numeric(dbf2$gisacres)), 
    high = muted("blue")) + 
  theme(axis.title = element_blank(), 
        axis.text = element_blank(), 
        axis.ticks = element_blank())

Hawaii zipcodes

Read the data

shape <- readOGR("zipcodes/zcta15.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "zipcodes/zcta15.shp", layer: "zcta15"
## with 94 features
## It has 23 fields

shape <- spTransform(shape, CRS("+proj=longlat"))
shape1 <- fortify(shape, region = 'zcta5ce10')
dbf <- read.dbf('zipcodes/zcta15.dbf')

Plot the population of each zipcode

ggplot() + 
  geom_map(data = dbf, 
           aes(map_id = zcta5ce10, 
               fill = as.numeric(Pop15)), 
           map = shape1) + 
  expand_limits(x = shape1$long, y = shape1$lat) + 
  scale_fill_gradient(
    name = 'Population',
    limits = c(0, max(as.numeric(dbf$Pop15))), 
    low = 'white',
    high = muted('red')) + 
  theme(axis.title = element_blank(), 
        axis.text = element_blank(), 
        axis.ticks = element_blank())

Plot the population for Oʻahu

shape2 <- fortify(shape, region = 'zcta5ce10') %>%
  filter(between(lat,  oahu_limits$min_lat,  oahu_limits$max_lat) & 
         between(long, oahu_limits$min_long, oahu_limits$max_long))
oahu_zipcodes <- shape2 %>% 
  group_by(id) %>%
  do(clong = centroid(cbind(.$long, .$lat))[1],
     clat  = centroid(cbind(.$long, .$lat))[2]) %>%
  mutate(clong = unlist(clong), clat=unlist(clat))
dbf2 <- dbf %>% 
  right_join(oahu_zipcodes, by = c('zcta5ce10' = 'id'))
ggplot() + 
  geom_map(data = dbf2, 
           aes(map_id = zcta5ce10, 
               fill = as.numeric(Pop15)),
           color = 'black',
           map = shape2) + 
  expand_limits(x = shape2$long, y = shape2$lat) + 
  scale_fill_gradient(
    name = 'Population',
    limits = c(0, max(as.numeric(dbf2$Pop15))), 
    low = 'white',
    high = muted('red')) + 
  geom_text_repel(data = oahu_zipcodes, aes(x = clong, y = clat, label = id, size = 0.05)) + 
  guides(size = 'none') +
  theme(axis.title = element_blank(), 
        axis.text = element_blank(), 
        axis.ticks = element_blank())