# 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)
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)
hist(sierraFebpts$ELEVATION) hist(sierraFebpts$LATITUDE)
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)
plotRGB(sierraBase)
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)

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_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) %>%
summarise(
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),
FACILITY_NAME = first(FACILITY_NAME),
COUNTY = first(COUNTY),
LATITUDE = first(LATITUDE),
LONGITUDE = first(LONGITUDE))
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) %>% rowid_to_column("TRI_ID") 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])) + labs(x='',y='') 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")) %>% mutate( 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") 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 : We’ll start with a regional view using a hillshade and terrain shading: tmap_mode("plot") 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) + tm_graticules(lines=F) Then zoom into Upper Sinking Cove to look at some water sample results: library(igisci) 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) library(terra) tmap_mode("plot") 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")) + tm_graticules(lines=F) An interactive map view of Upper Sinking Cove: tmap_mode("view") 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)
tmap_mode("plot")

### References

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.