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. You are free to approach this in your own way, but you will need to provide at least the following:

A map of ozone concentrations at the measurement stations

INSTALLING PACKAGES

library(ggplot2)
library(gstat)
library(RColorBrewer)
library(sf)
library(stars)
library(viridis)

DATA PREP

ozone.sf = st_read("../datafiles/mwozone/mwozone/ozone.shp", quiet = TRUE)
st_crs(ozone.sf) = 4326

states <- st_read("../datafiles/ne_50m_admin_1_states_provinces/ne_50m_admin_1_states_provinces/ne_50m_admin_1_states_provinces_shp.shp")
## Reading layer `ne_50m_admin_1_states_provinces_shp' from data source 
##   `N:\Projects\geog6000\datafiles\ne_50m_admin_1_states_provinces\ne_50m_admin_1_states_provinces\ne_50m_admin_1_states_provinces_shp.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 40 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -178.1945 ymin: -43.61946 xmax: 153.617 ymax: 83.11611
## Geodetic CRS:  WGS 84
st_crs(states) = 4326
lakes <- st_read("../datafiles/ne_50m_lakes/ne_50m_lakes/ne_50m_lakes.shp")
## Reading layer `ne_50m_lakes' from data source 
##   `N:\Projects\geog6000\datafiles\ne_50m_lakes\ne_50m_lakes\ne_50m_lakes.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 398 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -165.8985 ymin: -50.62002 xmax: 176.0827 ymax: 81.94033
## Geodetic CRS:  WGS 84
st_crs(lakes) = 4326
places <- st_read("../datafiles/ne_50m_populated_places/ne_50m_populated_places/ne_50m_populated_places.shp")
## Reading layer `ne_50m_populated_places' from data source 
##   `N:\Projects\geog6000\datafiles\ne_50m_populated_places\ne_50m_populated_places\ne_50m_populated_places.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1249 features and 92 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -175.2206 ymin: -90 xmax: 179.2166 ymax: 78.21668
## Geodetic CRS:  WGS 84
st_crs(places) = 4326

Mapping ozone concentrations with ggplot

ggplot() +
  geom_sf(data = lakes, fill = "lightblue") +
  geom_sf(data = states, fill = NA) +
  geom_sf(data = ozone.sf, aes(col = ozone), size = 2) +
  scale_color_viridis_c() +
  geom_label(data = places, aes(LATITUDE, LONGITUDE, label = NAME), size = 2.5) +
  coord_sf(xlim = c(-94, -82), ylim = c(36, 45), expand = FALSE) +
  theme_bw()

ozonecont = ozone.sf$ozone
greaterthan = ozonecont > 100
ozone.sf$ozone100 = greaterthan

ggplot() +
  geom_sf(data = ozone.sf, aes(col = ozone100), size = 2) +
  geom_label(data = places, aes(LATITUDE, LONGITUDE, label = NAME), size = 2.5) +
  coord_sf(xlim = c(-94, -82), ylim = c(36, 45), expand = FALSE) +
  theme_bw()

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

ozone100.var = variogram(ozone100 ~ 1, 
                         data = ozone.sf)
ozone100.vgm = vgm(0.07, "Sph", 120, 0.07)
plot(ozone100.var, ozone100.vgm)

It was hard to find a definite structure in the variogram indicating a geographic trend, so the sharp cutoff and short range is indicative of that.

A variogram model fitted to this (as a figure), and give the model chosen, the sill, range and nugget value (if appropriate)

My sill was .07, model was spherical, range was 120, and nugget was .07

ozone100.vgmfit = fit.variogram(ozone100.var, ozone100.vgm)

plot(ozone100.var, ozone100.vgmfit)

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. See code below for how to do this simply in R

mybbox <- st_bbox(ozone.sf)
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(ozone100 ~ 1,
                    ozone.sf,
                    pred.grid,
                    ozone100.vgmfit,
                    nmax = 40)
## [using ordinary kriging]
plot(ozone100.ik["var1.pred"], reset = FALSE)
plot(st_geometry(lakes), reset = FALSE, add = TRUE, col = "lightblue")
plot(st_geometry(states), reset = FALSE, add = TRUE)
plot(ozone.sf["ozone100"], add = TRUE, pch = 16)

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

Ozone tends to be in high concentration around Lake Michigan, and the dense presence of ozone concentrations above 100ppb heavily influences the probability estimated by the kriging function. That area is heavily industrial. There are areas with moderate probability of > 100ppb ozone despite the absence of any observations in some areas, indicating that the kriging is aware of the lack of data and thus says “it’s kind of a toss-up”