1 Abstract

Most data has a spatial dimension to it. Knowing ‘where’ the data is coming from is often as crucial as knowing the ‘what’, ‘when’ and ‘who’ dimensions of a given dataset. Therefore it should be no surprise that R has a rich suite of packages for constructing maps and analyzing spatial data. R’s capability has grown so much over the years that it’s functionality rivals many dedicated geographic information systems (GIS). During this Meetup the basics for managing and mapping spatial data with be introduced, using the following packages: sf, ggplot2, tmap, mapview and leaflet.

2 Example datasets

How you plot your data in R depends on what format it’s in. R has several different formats for managing spatial data (e.g. sp vs sf), and different formats are only compatible with certain plotting systems. Thankfully converting between different spatial formats is not difficult.

2.1 maps R package

library(maps)

map(database = "state")

map(database = "state", region = "indiana")

2.2 County Boundaries

library(USAboundaries)
library(sf)

cnty <- us_counties() 
cnty <- subset(cnty, !state_name %in% c("Alaska", "Hawaii"))

vars <- c("statefp", "countyfp")
plot(cnty[vars])

2.3 USDA-NASS Corn Yield Data

https://github.com/potterzot/rnassqs

devtools::install_github('potterzot/rnassqs')

library(rnassqs)

source("C:/Users/steph/Nextcloud/code/api_keys.R")

data(state)
st <- state.abb

corn_us <- lapply(st, function(x) {
  cat("getting", x, as.character(Sys.time()), "\n")
  tryCatch({
    corn = nassqs_yield(
      list("commodity_desc"="CORN",
           "agg_level_desc"="COUNTY",
           "state_alpha"=x
           ),
      key = nass_key
      )},
    error = function(err) {
      print(paste("Error occured:  ",err))
      return(NULL)
      }
    )
  })
corn_us <- do.call("rbind", corn_us)

save(corn_us, file = "C:/workspace2/corn_us.RData")
write.csv(corn_us, file = "nass_corn_us.csv", row.names = FALSE)
load(file = "C:/Users/Stephen.Roecker/Nextcloud/data/corn_us.RData")
corn_yield <- subset(corn_us, short_desc == "CORN, GRAIN - YIELD, MEASURED IN BU / ACRE")
corn_yield <- within(corn_yield, {
  Value      = as.numeric(Value)
  year       = as.numeric(year)
  state_name = NULL
  state      = state_alpha 
  })


cnty_corn <- merge(cnty, corn_yield, 
                   by.x = c("state_abbr", "countyfp"), 
                   by.y = c("state_alpha", "county_code"), 
                   all.x = TRUE
                   )

corn_states <- c("IL", "IA", "IN", "MI", "MN", "MO", "NE", "OH", "SD", "ND", "WI")

library(dplyr)
library(ggplot2)

group_by(corn_yield, state_alpha, year) %>%
  summarize(
    yield_low    = min(Value, na.rm = TRUE),
    yield_median = median(Value, na.rm = TRUE),
    yield_max    = max(Value, na.rm = TRUE)
    ) %>%
  filter(state_alpha %in% corn_states) %>%
  mutate(state = state_alpha,
         source = "NASS"
         ) %>%
  ggplot() +
  geom_line(aes(x = year, y = yield_median, col = source)) +
  geom_ribbon(aes(x = year, ymin = yield_low, ymax = yield_max), alpha = 0.25) +
  facet_wrap(~ state) +
  ylab("yield per county (bu/acre)") +
  ggtitle("USDA-NASS Corn Yields")

  # geom_point(data = yld_sum[yld_sum$state %in% corn_states, ], aes(x = 2018, y = yield_med, col = "NASIS"), size = 1) +
  # geom_ribbon(data = yld_sum2[yld_sum2$state %in% corn_states, ], aes(x = year, ymin = yield_low2, ymax = yield_max2, col = "NASIS"), alpha = 0.25) +
  # geom_pointrange(data = yld_sum[yld_sum$state %in% corn_states, ], aes(x = 2018, y = yield_med, ymin = yield_low2, ymax = yield_max2, col = source))

2.4 Indiana General Soil Map

IN <- sf::read_sf(dsn = "D:/geodata/soils/soils_GSMCLIP_mbr_2599033_03/wss_gsmsoil_IN_[2006-07-06]/spatial/gsmsoilmu_a_in.shp", layer = "gsmsoilmu_a_in")

# simplify polygons
IN <- rmapshaper::ms_simplify(IN)

2.5 Kellogg Soil Survey Laboratory Data

library(aqp)
library(soilDB)

# download lab locations for the Miami soil series
miami <- fetchKSSL("Miami")
miami <- site(miami)
miami <- subset(miami, complete.cases(x, y))
miami <- within(miami, {
  lon = x
  lat = y
  })
head(miami)
##   pedon_key pedlabsampnum      pedon_id pscs_top pscs_bottom
## 2      7153       79P0060    78IN177007       15          65
## 3      7154       79P0061    78IN177008       30          80
## 4      8836       81P0725 S1981IN011002       20          70
## 5      8840       81P0729 S1981IN011006       20          70
## 6      8846       81P0735 S1981IN011012       25          75
## 7      8850       81P0739 S1981IN011016       15          65
##   pedon_completeness_index   state mlra geomposhill geomposmntn
## 2                        1 Indiana 111A                      NA
## 3                        1 Indiana 111A                      NA
## 4                        7 Indiana 111A                      NA
## 5                        7 Indiana 111A                      NA
## 6                        6 Indiana 111A                      NA
## 7                        7 Indiana 111A                      NA
##   geompostrce geomposflats hillslopeprof geomslopeseg pmgroupname
## 2          NA           NA                                       
## 3          NA           NA                                       
## 4          NA           NA                                       
## 5          NA           NA                                       
## 6          NA           NA                                       
## 7          NA           NA                                       
##     drainagecl taxonname taxorder taxsuborder taxgrtgroup
## 2                  Miami Alfisols      Udalfs  Hapludalfs
## 3                  Miami Alfisols      Udalfs  Hapludalfs
## 4 Well drained     Miami Alfisols      Udalfs  Hapludalfs
## 5 Well drained     Miami Alfisols      Udalfs  Hapludalfs
## 6 Well drained     Miami Alfisols      Udalfs  Hapludalfs
## 7 Well drained     Miami Alfisols      Udalfs  Hapludalfs
##             taxsubgrp taxpartsize taxpartsizemod taxceactcl taxreaction
## 2    Typic Hapludalfs  fine-loamy                                      
## 3    Typic Hapludalfs  fine-loamy                              not used
## 4 Oxyaquic Hapludalfs  fine-loamy                    active            
## 5 Oxyaquic Hapludalfs  fine-loamy                    active            
## 6 Oxyaquic Hapludalfs  fine-loamy                    active            
## 7 Oxyaquic Hapludalfs  fine-loamy                    active            
##   taxtempcl taxmoistscl taxtempregime soiltaxedition osdtypelocflag
## 2     mesic                     mesic                             0
## 3                               mesic                             0
## 4     mesic    Oxyaquic         mesic  tenth edition              0
## 5     mesic    Oxyaquic         mesic  tenth edition              0
## 6     mesic    Oxyaquic         mesic  tenth edition              0
## 7     mesic    Oxyaquic         mesic  tenth edition              0
##           x        y      lat       lon
## 2 -85.02861 39.78639 39.78639 -85.02861
## 3 -85.01667 39.77750 39.77750 -85.01667
## 4 -86.27920 39.98251 39.98251 -86.27920
## 5 -86.27920 39.98251 39.98251 -86.27920
## 6 -86.27920 39.98251 39.98251 -86.27920
## 7 -86.27920 39.98251 39.98251 -86.27920

3 Construct and convert spatial objects

3.1 point objects

# construct

# sp object

library(sp)

miami_sp <- SpatialPointsDataFrame(
  data   = miami,
  coords = cbind(miami$lon, miami$lat),
  proj4string = CRS("+init=epsg:4326")
  )


#  sf object

library(sf)

miami_sf <- st_as_sf(
  miami,
  coords = c("lon", "lat"),
  crs = 4326
  )


# data structures
str(miami_sp, 2)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  56 obs. of  34 variables:
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:56, 1:2] -85 -85 -86.3 -86.3 -86.3 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ bbox       : num [1:2, 1:2] -89.3 39 -83.9 42.4
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
str(miami_sf)
## Classes 'sf' and 'data.frame':   56 obs. of  33 variables:
##  $ pedon_key               : chr  "7153" "7154" "8836" "8840" ...
##  $ pedlabsampnum           : chr  "79P0060" "79P0061" "81P0725" "81P0729" ...
##  $ pedon_id                : chr  "78IN177007" "78IN177008" "S1981IN011002" "S1981IN011006" ...
##  $ pscs_top                : int  15 30 20 20 25 15 NA NA NA NA ...
##  $ pscs_bottom             : int  65 80 70 70 75 65 NA NA NA NA ...
##  $ pedon_completeness_index: int  1 1 7 7 6 7 6 6 6 6 ...
##  $ state                   : chr  "Indiana" "Indiana" "Indiana" "Indiana" ...
##  $ mlra                    : chr  "111A" "111A" "111A" "111A" ...
##  $ geomposhill             : chr  "" "" "" "" ...
##  $ geomposmntn             : logi  NA NA NA NA NA NA ...
##  $ geompostrce             : logi  NA NA NA NA NA NA ...
##  $ geomposflats            : logi  NA NA NA NA NA NA ...
##  $ hillslopeprof           : chr  "" "" "" "" ...
##  $ geomslopeseg            : chr  "" "" "" "" ...
##  $ pmgroupname             : chr  "" "" "" "" ...
##  $ drainagecl              : chr  "" "" "Well drained" "Well drained" ...
##  $ taxonname               : chr  "Miami" "Miami" "Miami" "Miami" ...
##  $ taxorder                : chr  "Alfisols" "Alfisols" "Alfisols" "Alfisols" ...
##  $ taxsuborder             : chr  "Udalfs" "Udalfs" "Udalfs" "Udalfs" ...
##  $ taxgrtgroup             : chr  "Hapludalfs" "Hapludalfs" "Hapludalfs" "Hapludalfs" ...
##  $ taxsubgrp               : chr  "Typic Hapludalfs" "Typic Hapludalfs" "Oxyaquic Hapludalfs" "Oxyaquic Hapludalfs" ...
##  $ taxpartsize             : chr  "fine-loamy" "fine-loamy" "fine-loamy" "fine-loamy" ...
##  $ taxpartsizemod          : chr  "" "" "" "" ...
##  $ taxceactcl              : chr  "" "" "active" "active" ...
##  $ taxreaction             : chr  "" "not used" "" "" ...
##  $ taxtempcl               : chr  "mesic" "" "mesic" "mesic" ...
##  $ taxmoistscl             : chr  "" "" "Oxyaquic" "Oxyaquic" ...
##  $ taxtempregime           : chr  "mesic" "mesic" "mesic" "mesic" ...
##  $ soiltaxedition          : chr  "" "" "tenth edition" "tenth edition" ...
##  $ osdtypelocflag          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ x                       : num  -85 -85 -86.3 -86.3 -86.3 ...
##  $ y                       : num  39.8 39.8 40 40 40 ...
##  $ geometry                :sfc_POINT of length 56; first list element:  'XY' num  -85 39.8
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  "pedon_key" "pedlabsampnum" "pedon_id" "pscs_top" ...
# convert

miami_sf <- st_as_sf(miami_sp)
miami_sp <- as(miami_sf, "Spatial")

3.2 Convert map objects

library(maptools)

st <- map("state", fill = TRUE, plot = FALSE)

# convert to sp object
st_sp <- map2SpatialPolygons(st, IDs = st$names)
proj4string(st_sp) <- CRS("+init=epsg:4326")
st_sp$state <- st$names


# loads data with ggplot2 package
st <- ggplot2::map_data("state")

3.3 Convert with broom package

library(broom)

st_tidy <- tidy(st_sp, region = "state") 
IN_tidy <- tidy(as(IN, "Spatial"), region = "MUSYM")

4 Example maps with several R packages

4.1 ggplot2

ggplot2 plots data frames, but can also use sf objects (which are a special case of data frames). It can plot rasters, but only if they are converted to a data frame.

library(ggplot2)


# Lines
ggplot() +
  geom_point(data = miami, aes(x = lon, y = lat)) +
  geom_path(data = st, aes(x = long, y = lat, group = group)) +
  xlim(range(miami$lon)) +
  ylim(range(miami$lat)) +
  ggtitle("Location of Miami Lab Pedons")

# Polygons
ggplot() +
  geom_polygon(data = st_tidy, aes(x = long, y = lat, group = group, fill = id)) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  # remove legend
  guides(fill = FALSE)

# sf objects

# Polygons
ggplot() +
  geom_sf(data = cnty, aes(fill = statefp, lty = NA)) +
  geom_sf(data = miami_sf) +
  coord_sf(crs = "+init=epsg:5070") +
  guides(fill = FALSE)

# Facets

test <- subset(cnty_corn, year %in% 2012:2017)

ggplot() +
  geom_sf(data = test, aes(fill = Value, lty = NA)) +
  scale_fill_viridis_c(na.value = "transparent") +
  facet_wrap(~ year) +
  geom_path(data = st, aes(x = long, y = lat, group = group)) +
  ggtitle(corn_yield$short_desc[1])

4.2 ggmap

ggmap expands ggplot2 to download and plot base maps.

library(ggmap)

# build bound box and get base map via ggmap
bb <- sf::st_bbox(IN)
bb <- make_bbox(lon = bb[c(3, 1)], lat = bb[c(2, 4)])
gmap <- get_map(bb, maptype = "terrain", source = "osm")


# Lines
ggmap(gmap) +
  geom_path(data = IN_tidy, aes(x = long, y = lat, group = group))

  # geom_sf(data = IN, fill = NA, inherit.aes = FALSE) +
  # guides(fill = FALSE)

# geom_sf() doesn't work with ggmp, their is a systematic shift https://github.com/r-spatial/sf/issues/336

4.3 tmap

The “t” in tmap stands for thematic, but tmap can also plot rasters natively. tmap’s syntax is very similar to ggplot2, but with a few twists.

library(tmap)

tm_shape(IN) + tm_polygons("MUSYM", border.col = NULL) +
  tm_shape(cnty) + tm_borders() +
  tm_shape(miami_sf) + tm_dots() +
  tm_legend(legend.outside = TRUE)

# interactive web map 
tmap_mode("view")

tm_basemap("OpenStreetMap") +
  tm_shape(IN) + tm_borders()

4.4 mapview

library(mapview)

cols <- RColorBrewer::brewer.pal(50, "Paired")

test <- mapview(IN, zcol = "MUSYM", lwd = 0, col.regions = cols) +
  mapview(cnty, type = "l")

test
# export to html
mapshot(test, url = "C:/workspace2/test.html", selfcontained = FALSE)

4.5 leaflet

library(leaflet)

test <- leaflet() %>%
  addProviderTiles("Esri.WorldImagery", group = "Imagery") %>%
  addPolygons(data = IN, fill = FALSE, color = "black", weight = 2)

test
# export to html
htmlwidgets::saveWidget(test, file = "C:/workspace2/test.html", selfcontained = FALSE)

5 Additional Reading

Healy, K., 2018. Data Visualization: a practical introduction. Princeton University Press. http://socviz.co/

Gimond, M., 2019. Intro to GIS and Spatial Analysis. https://mgimond.github.io/Spatial/

Hijmans, R.J., 2019. Spatial Data Science with R. https://rspatial.org/

Lovelace, R., J. Nowosad, and J. Muenchow, 2019. Geocomputation with R. CRC Press. https://bookdown.org/robinlovelace/geocompr/

Pebesma, E., and R. Bivand, 2019. Spatial Data Science. https://keen-swartz-3146c4.netlify.com/