10 Spatial Analysis case studies
10.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)
<- st_as_sf(read_csv(ex("sierra/sierraFeb.csv")),
sierraFebpts 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)
<- sierraFebpts %>%
sierraFebpts2000 filter(ELEVATION >= 2000)
<- get_tiles(sierraFebpts2000)
sierraBase plotRGB(sierraBase)
<- vect(sierraFebpts2000)
pts 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.
10.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)
<- read_csv(ex("TRI/TRI_2017_CA.csv"))
TRI_CA <- TRI_CA %>%
TRI_BySite 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))
<- st_read(ex("BayArea/BayAreaCounties.shp"))
co <- st_read(ex("CA/CAfreeways.shp"))
freeways <- st_make_valid(st_read(ex("BayArea/BayAreaTracts.shp")))
censusBayArea <- read_csv(ex("CA/CA_MdInc.csv")) %>%
incomeByTract ::select(trID, HHinc2016) %>%
dplyrmutate(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"))
<- st_centroid(censusBayArea)
censusCentroids <- censusBayArea
census <- st_bbox(co)
bnd <- TRI_BySite %>%
TRIdata filter(air_releases > 0) %>%
filter((LONGITUDE < bnd[3]) & (LONGITUDE > bnd[1]) & (LATITUDE < bnd[4]) & (LATITUDE > bnd[2]))
<- co %>% dplyr::select(NAME, CNTY_FIPS) %>% st_set_geometry(NULL)
BA_counties <- st_as_sf(TRIdata, coords = c("LONGITUDE", "LATITUDE"), crs=4326) %>%
TRI_sp st_join(censusBayArea) %>%
filter(CNTY_FIPS %in% BA_counties$CNTY_FIPS)
<- TRI_sp %>%
TRI2join 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
<- read_csv(ex("TRI/CDC_health_data_by_Tract_CA_2018_release.csv")) %>%
CDCdata 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)
<- st_as_sf(CDCdata, coords = c('lon', 'lat'), crs=4326)
CDCbay
# D to TRI: create vector of index of nearest TRI facility to each CDC (tract centroid) point
<- st_nearest_feature(CDCbay, TRI_sp)
nearest_TRI_ids # get locations of each NEAREST TRI facility only
<- TRI_sp$geometry[nearest_TRI_ids]
near_TRI_loc
<- 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"))
<- left_join(CDCbay, TRI2join, by = c("nearest_TRI" = "TRI_ID")) %>%
CDCTRIbay 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.)
10.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):
We’ll start with a regional view using a hillshade and terrain shading:
tmap_mode("plot")
<- rast(ex("SinkingCove/DEM_SinkingCoveUTM.tif"))
DEM <- terrain(DEM, v='slope')
slope <- terrain(DEM, v='aspect')
aspect <- shade(slope/180*pi, aspect/180*pi, angle=40, direction=330)
hillsh # Need to crop a bit since grid north != true north
<- st_bbox(DEM)
bbox0 <- bbox0$xmax - bbox0$xmin
xrange <- bbox0$ymax - bbox0$ymin
yrange <- bbox0
bbox1 <- 0.05
crop 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
bbox1[<- bbox1 %>% st_as_sfc() # makes a polygon
bboxPoly #
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)
<- read_excel(ex("SinkingCove/SinkingCoveWaterChem.xlsx")) %>%
wChemData mutate(siteLoc = str_sub(Site,start=1L, end=1L))
<- wChemData %>% filter(siteLoc == "T") %>% mutate(siteType = "trunk")
wChemTrunk <- wChemData %>% filter(siteLoc %in% c("D","S")) %>% mutate(siteType = "dripwater")
wChemDrip <- wChemData %>% filter(siteLoc %in% c("B", "F", "K", "W", "P")) %>%
wChemTrib mutate(siteType = "tributary")
<- bind_rows(wChemTrunk, wChemDrip, wChemTrib)
wChemData <- read_csv(ex("SinkingCove/SinkingCoveSites.csv"))
sites <- wChemData %>%
wChem left_join(sites, by = c("Site" = "site")) %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
library(terra)
tmap_mode("plot")
<- rast(ex("SinkingCove/DEM_SinkingCoveUTM.tif"))
DEM <- terrain(DEM, v='slope')
slope <- terrain(DEM, v='aspect')
aspect <- shade(slope/180*pi, aspect/180*pi, 40, 330)
hillsh <- st_bbox(wChem)
bounds <- bounds$xmax - bounds$xmin
xrange <- bounds$ymax - bounds$ymin
yrange <- as.numeric(bounds$xmin - xrange/10)
xMIN <- as.numeric(bounds$xmax + xrange/10)
xMAX <- as.numeric(bounds$ymin - yrange/10)
yMIN <- as.numeric(bounds$ymax + yrange/10)
yMAX <- st_bbox(c(xmin=xMIN, xmax=xMAX, ymin=yMIN, ymax=yMAX), crs= st_crs(4326))
newbounds 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")
<- st_bbox()
bounds <- filter(wChem, Month == 8)
wChem2map <- min(wChem2map$TH); maxVal <- max(wChem2map$TH)
minVal 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")