9 Spatial Analysis case studies

9.1 Plotting filtered map data using base plot, terra::plotRGB and SpatVector data

In the Spatial Analysis chapter, we’ve mainly used tmap, but we might also want to see how we’d plot the same thing in the base plot system, which also support’s terra’s SpatVector data. So in this brief example, we’ll create a SpatVector after a data filter, then plot in the base plot system.

library(tidyverse); library(sf); library(igisci)
sierraFebpts <- st_as_sf(read_csv(ex("sierra/sierraFeb.csv")), 
                         coords = c("LONGITUDE", "LATITUDE"), remove=FALSE, crs=4326)

We might also look at how we get information about out data to decide what we might want to plot on the map. Since our goal is to filter for a range of elevations and latitudes, we might look at a couple of histograms, and we’ll use a base R plot’s parameter setting (par) to create a side by side plot (1 row, 2 columns).

par(mfrow=c(1,2))  # lets base plot functions create 2 adjacent plots (1 row, 2 columns)
Histograms of Sierra station elevations and latitudes

Figure 9.1: Histograms of Sierra station elevations and latitudes

par(mfrow=c(1,1)) # need to reset to defaults

Now let’s continue with that information where we’ve gotten an idea about what kinds of elevations occur at the stations, and make a map that shows all stations > 2000 m elevation.

We’ll briefly use terra::plotRGB to add a basemap to a static map. Warning: the get_tiles function goes online to get the basemap data, so if you don’t have a good internet connection or the site goes down, this may fail. We’ll also modify the station names and set the text label position as pos=3, which means centered on top (see ?graphics::text).

library(terra); library(maptiles)
sierraFebpts2000 <- sierraFebpts %>%
  filter(ELEVATION >= 2000)
sierraBase <- get_tiles(sierraFebpts2000)
pts <- vect(sierraFebpts2000)
points(pts, col="red")
text(pts, labels=str_sub(pts$STATION_NAME, 1, str_locate(pts$STATION_NAME, ",")-1), cex = 0.7, pos=3)
Plotting filtered data (ELEVATION >= 2000)

Figure 9.2: Plotting filtered data (ELEVATION >= 2000)

Before we continue, let’s look at some of the parameters just used for sizing and placing text labels, and figure out what they do:

text(pts, labels=str_sub(pts$STATION_NAME, 1, str_locate(pts$STATION_NAME, ",")-1), cex = 0.7, pos=3)

  • What does cex do? (Use ?par or ?text then search within that help)
  • What does pos do? (Use ?text)
  • What did the stringr methods do?

Tweak each of the above to see the effect on the map.

Maybe there’s a way with plotRGB, but it’s limited in not apparently having a way to provide a graticule of latitude and longitude, so its use as a map is limited, and there are generally better options in tmap.

9.2 Distance from feature to closest feature: Air quality / TRI study

In the Spatial Analysis chapter section on Distance and specifically Distance to Closest Feature, we looked at point to closest point distances. Here’s another example of this that attempts to find the output of closest TRI facility to a census tract where we might compare health or socioeconomic variables to that value. The igisci::BayAreaTracts data has census data from 2010 by census tracts in the SF Bay area. In addition to the distance functions, the code below includes a few useful sf functions applied to that data.

library(igisci); library(sf); library(tidyverse)
TRI_CA <- read_csv(ex("TRI/TRI_2017_CA.csv"))
TRI_BySite <- TRI_CA %>%
  mutate(all_air = `5.1_FUGITIVE_AIR` + `5.2_STACK_AIR`) %>%
  mutate(carcin_release = all_air * (CARCINOGEN == "YES")) %>%
  group_by(TRI_FACILITY_ID) %>%
    count = n(),
    air_releases = sum(all_air, na.rm = TRUE),
    fugitive_air = sum(`5.1_FUGITIVE_AIR`, na.rm = TRUE),
    stack_air = sum(`5.2_STACK_AIR`, na.rm = TRUE),
    COUNTY = first(COUNTY),
co <- st_read(ex("BayArea/BayAreaCounties.shp"))        
freeways <- st_read(ex("CA/CAfreeways.shp"))          
censusBayArea <- st_make_valid(st_read(ex("BayArea/BayAreaTracts.shp"))) 
incomeByTract <- read_csv(ex("CA/CA_MdInc.csv")) %>%
  select(trID, HHinc2016) %>%
  mutate(HHinc2016 = as.numeric(str_c(HHinc2016)),
         joinid = str_c("0", trID))
censusBayArea <- censusBayArea %>%
  mutate(whitepct = WHITE / POP2010 * 100) %>%
  mutate(People_of_Color_pct = 100 - whitepct) %>%
  left_join(incomeByTract, by = c("FIPS" = "joinid"))
censusCentroids <- st_centroid(censusBayArea)
census <- censusBayArea
bnd <- st_bbox(co)
TRIdata <- TRI_BySite %>%
  filter(air_releases > 0) %>%
  filter((LONGITUDE < bnd[3]) & (LONGITUDE > bnd[1]) & (LATITUDE < bnd[4]) & (LATITUDE > bnd[2]))
BA_counties <- co %>% select(NAME, CNTY_FIPS) %>% st_set_geometry(NULL)
TRI_sp <- st_as_sf(TRIdata, coords = c("LONGITUDE", "LATITUDE"), crs=4326) %>%
  st_join(censusBayArea) %>%
  filter(CNTY_FIPS %in% BA_counties$CNTY_FIPS)
TRI2join <- TRI_sp %>% 
  st_set_geometry(NULL) %>%
ggplot() +
  geom_sf(data = co, fill = NA) +
  geom_sf(data=TRI_sp) +
  coord_sf(xlim = c(bnd[1], bnd[3]), ylim = c(bnd[2], bnd[4])) +
EPA Toxic Release Inventory (TRI) facilities 2017

Figure 9.3: EPA Toxic Release Inventory (TRI) facilities 2017

Asthma model related to distance to nearest TRI facility.

# CDC data
CDCdata <- read_csv(ex("TRI/CDC_health_data_by_Tract_CA_2018_release.csv")) %>%
    lon = as.numeric(str_sub(Geolocation, 17, 30)),
    lat = as.numeric(str_sub(Geolocation, 2, 14)),
    CDC_CNTY_FIPS = str_sub(TractFIPS, 2, 4)) %>%
  filter(CDC_CNTY_FIPS %in% BA_counties$CNTY_FIPS)
CDCbay <- st_as_sf(CDCdata, coords = c('lon', 'lat'), crs=4326)

# D to TRI:  create vector of index of nearest TRI facility to each CDC (tract centroid) point
nearest_TRI_ids <- st_nearest_feature(CDCbay, TRI_sp)
# get locations of each NEAREST TRI facility only
near_TRI_loc <- TRI_sp$geometry[nearest_TRI_ids]

CDCbay <- CDCbay %>%
  mutate(d2TRI = st_distance(CDCbay, near_TRI_loc, by_element=TRUE),
         d2TRI = units::drop_units(d2TRI),
         nearest_TRI = nearest_TRI_ids) %>%
  left_join(BA_counties, by = c("CDC_CNTY_FIPS" = "CNTY_FIPS"))

CDCTRIbay <- left_join(CDCbay, TRI2join, by = c("nearest_TRI" = "TRI_ID")) %>%  
  filter(d2TRI > 0) # %>%
CDCTRIbay %>%
  ggplot(aes(air_releases, CASTHMA_CrudePrev)) +
  geom_point(aes(color = NAME)) + geom_smooth(method="lm", se=FALSE, color = "black") +
  labs(x = "Air releases at nearest TRI facility", 
       y = "Asthma Prevalence CDC Model, Adults >= 18, 2016")
Asthma prevalence related to nearest TRI total releases

Figure 9.4: Asthma prevalence related to nearest TRI total releases

This is an admittedly crude analysis, lumping all releases together, and so it’s not surprising that there’s no clear relationship over the whole Bay Area. There’s quite a lot you could explore in the TRI data, and you might want to focus on one county, specific toxins, or probably aggregate releases at the total output in the area. And analyzing health data is very challenging given the wide range of factors influencing patterns. You might also consider looking at earlier TRI data (the igisci extdata "TRI" folder has it back to 1987, for California.)

9.3 Exploring terrain and water samples in a karst area

Sinking Cove, Tennessee is a karst valley system carved into the Cumberland Plateau, a nice place to see the use of a hillshade raster created from a digital elevation model using raster or terra functions for slope, aspect, and hillshade. We’ll look at this area and some results from a hydrologic and geochemical study (Jerry D. Davis and Brook 1993):

Cave Cove Cave Entrance

Figure 9.5: Cave Cove Cave Entrance

Sinking Cove Cave

Figure 9.6: Sinking Cove Cave

Sinking Cove Cave spring

Figure 9.7: Sinking Cove Cave spring

We’ll start with a regional view using a hillshade and terrain shading:

DEM <- rast(ex("SinkingCove/DEM_SinkingCoveUTM.tif"))
slope <- terrain(DEM, v='slope')
aspect <- terrain(DEM, v='aspect')
hillsh <- shade(slope/180*pi, aspect/180*pi, angle=40, direction=330)
# Need to crop a bit since grid north != true north 
bbox0 <- st_bbox(DEM)
xrange <- bbox0$xmax - bbox0$xmin
yrange <- bbox0$ymax - bbox0$ymin
bbox1 <- bbox0
crop <- 0.05
bbox1[1] <- bbox0[1] + crop * xrange # xmin
bbox1[3] <- bbox0[3] - crop * xrange # xmax
bbox1[2] <- bbox0[2] + crop * yrange # ymin
bbox1[4] <- bbox0[4] - crop * yrange # ymax
bboxPoly <- bbox1 %>% st_as_sfc() # makes a polygon
tm_shape(hillsh, bbox=bboxPoly) +
  tm_raster(palette="-Greys",legend.show=F,n=20) +
  tm_shape(DEM) +
  tm_raster(palette=terrain.colors(24), alpha=0.5) +
Sinking Cove Area map, using terra rasters in tmap

Figure 9.8: Sinking Cove Area map, using terra rasters in tmap

Then zoom into Upper Sinking Cove to look at some water sample results:

library(sf); library(tidyverse); library(readxl); library(tmap)
wChemData <- read_excel(ex("SinkingCove/SinkingCoveWaterChem.xlsx")) %>%
  mutate(siteLoc = str_sub(Site,start=1L, end=1L))
wChemTrunk <- wChemData %>% filter(siteLoc == "T") %>% mutate(siteType = "trunk")
wChemDrip <- wChemData %>% filter(siteLoc %in% c("D","S")) %>% mutate(siteType = "dripwater")
wChemTrib <- wChemData %>% filter(siteLoc %in% c("B", "F", "K", "W", "P")) %>% 
  mutate(siteType = "tributary")
wChemData <- bind_rows(wChemTrunk, wChemDrip, wChemTrib)
sites <- read_csv(ex("SinkingCove/SinkingCoveSites.csv"))
wChem <- wChemData %>%
  left_join(sites, by = c("Site" = "site")) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
DEM <- rast(ex("SinkingCove/DEM_SinkingCoveUTM.tif"))   
slope <- terrain(DEM, v='slope')
aspect <- terrain(DEM, v='aspect')
hillsh <- shade(slope/180*pi, aspect/180*pi, 40, 330) 
bounds <- st_bbox(wChem)
xrange <- bounds$xmax - bounds$xmin
yrange <- bounds$ymax - bounds$ymin
xMIN <- as.numeric(bounds$xmin - xrange/10)
xMAX <- as.numeric(bounds$xmax + xrange/10)
yMIN <- as.numeric(bounds$ymin - yrange/10)
yMAX <- as.numeric(bounds$ymax + yrange/10)
newbounds <- st_bbox(c(xmin=xMIN, xmax=xMAX, ymin=yMIN, ymax=yMAX), crs= st_crs(4326))
tm_shape(hillsh,bbox=newbounds) +
  tm_raster(palette="-Greys",legend.show=F,n=20) +
  tm_shape(DEM) + tm_raster(palette=terrain.colors(24), alpha=0.5,legend.show=F) +
  tm_shape(wChem) + tm_symbols(size="TH", col="Lithology", scale=2, shape="siteType") +
  tm_layout(legend.position = c("left", "bottom")) +
Upper Sinking Cove water sample results

Figure 9.9: Upper Sinking Cove water sample results

An interactive map view of Upper Sinking Cove:

bounds <- st_bbox()
wChem2map <- filter(wChem, Month == 8)
minVal <- min(wChem2map$TH); maxVal <- max(wChem2map$TH)
tm_basemap(leaflet::providers$Esri.WorldTopoMap) + 
  tm_shape(wChem2map) + tm_symbols(col="siteType", size="TH", scale=2) +
  tm_layout(title=paste("Total Hardness ",as.character(minVal),"-",as.character(maxVal)," mg/L", sep=""))
tm_basemap(leaflet::providers$Esri.WorldTopoMap) + 
  tm_shape(wChem2map) + tm_symbols(col="Lithology", size="TH", scale=2)


Davis, Jerry D, and George A Brook. 1993. “Geomorphology and Hydrology of Upper Sinking Cove, Cumberland Plateau, Tennessee.” Earth Surface Processes and Landforms 18 (4): 339–62.