Lab 13

ozone <- st_read("./mwozone/ozone.shp", quiet = TRUE)
st_crs(ozone) <- 4326

states <- st_read("./ne_50m_admin_1_states_provinces/ne_50m_admin_1_states_provinces_shp.shp", quiet = TRUE)
lakes <- st_read("./ne_50m_lakes/ne_50m_lakes.shp", quiet = TRUE)
places <- st_read("./ne_50m_populated_places/ne_50m_populated_places.shp", quiet = TRUE)

Instructions: The compressed file mwozone.zip contains measurements of ground ozone concentration in the midwestern US for June 20, 1987 in the shapefile ozone.shp. The goal of this exercise is use this data to produce a map showing the probability that a certain threshold (100 ppb) was passed on that date. You can do this using either indicator kriging or indicator simulation, using the examples provided in the lab. - A map of ozone concentrations at the measurement stations

head(ozone)

## Simple feature collection with 6 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -91.404 ymin: 39.933 xmax: -87.546 ymax: 41.978
## Geodetic CRS:  WGS 84
##       lon    lat    ozone ozone100 coords_x1 coords_x2               geometry
## 1 -91.404 39.933  75.0000        F   -91.404    39.933 POINT (-91.404 39.933)
## 2 -88.230 40.124  84.2500        F   -88.230    40.124  POINT (-88.23 40.124)
## 3 -87.546 41.757  90.8750        F   -87.546    41.757 POINT (-87.546 41.757)
## 4 -87.671 41.978 127.4286        T   -87.671    41.978 POINT (-87.671 41.978)
## 5 -87.568 41.708 104.5000        T   -87.568    41.708 POINT (-87.568 41.708)
## 6 -87.726 41.739  86.2500        F   -87.726    41.739 POINT (-87.726 41.739)

ozone$threshold <- ozone$ozone > 100

par(mfrow = c(1, 2))
plot(st_geometry(ozone), reset = FALSE, main = "Ozone Concentrations (ppb)")
plot(st_geometry(lakes), reset = FALSE, add = TRUE)
plot(st_geometry(states), reset = FALSE, add = TRUE)
plot(ozone["ozone"], border = "black", add = TRUE, pch = 16, cex = 1.5)

plot(st_geometry(ozone), reset = FALSE, main = "Ozone Concentrations\n Greater Than 100 ppb")
plot(st_geometry(lakes), reset = FALSE, add = TRUE)
plot(st_geometry(states), reset = FALSE, add = TRUE)
plot(ozone["threshold"], add = TRUE, pch = 16, cex = 1.5)

• A plot of the sample variogram describing the spatial structure of the condition you want to interpolate (ozone > 100 ppb), and a brief description of the spatial structure you observe
• A variogram model fitted to this (as a figure), and give the model chosen, the sill, range and nugget value (if appropriate)

ozone100.var <- variogram(threshold ~ 1, ozone)

modNugget <- 0.095
modRange <- 125
modSill <- 0.07

ozone100.vgm <- vgm(modNugget, "Sph", modRange, modSill)
plot(ozone100.var, ozone100.vgm, main = "Variogram Model")

• A map of predicted probabilities of passing 100 ppb ozone. In order to do this, you will need to create a longitude/latitude grid for interpolation.

mybbox <- st_bbox(ozone)
pred.grid <- st_as_stars(mybbox, 
                         xlim = c(-94, -82), 
                         ylim = c(36, 45), 
                         dx = 0.1, 
                         dy = 0.1)
st_crs(pred.grid) <- 4326

ozone100.ik <- krige(threshold ~ 1,
                     ozone,
                     pred.grid,
                     ozone100.vgm,
                     nmax = 50)

## [using ordinary kriging]

plot(ozone100.ik["var1.pred"])

ozone100.ik$var1.pred[which(ozone100.ik$var1.pred < 0)] <- 0
my.pal <- viridis::viridis(10)

par(oma = c(0,0,2,0))
plot(ozone100.ik["var1.pred"],
     col = my.pal,
     breaks = seq(0, 1, length.out = 11),
     main = NA)
title("Probability of Ozone Levels Exceeding \n 100 ppb in the U.S. Midwest", outer = TRUE)

• A brief description of the spatial pattern shown on the interpolated map — what do you think the higher probabilities relate to?

Answer: Likelihood of exceeding the ozone concentration threshold of 100 ppb is higher near Columbus, OH; Cincinnati, OH; Peoria, IL; Grand Rapids, MI; and along the Gold Coast Illinois and Wisconsin. Ozone levels are likely related to human activity, as these are all urban areas with high concentrations of people (in 1990, Peoria was the second-most populous city in Illinois.)