set.seed(436)

Load libraries

library(ggplot2)
library(gstat)
library(RColorBrewer)
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(stars)
## Loading required package: abind
## Registered S3 methods overwritten by 'stars':
##   method             from
##   st_bbox.SpatRaster sf  
##   st_crs.SpatRaster  sf
library(viridis)
## Loading required package: viridisLite
library(tmap)

Read in data

ozone.sf <- st_read("datafiles/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.shp")
## Reading layer `ne_50m_admin_1_states_provinces' from data source 
##   `/Users/jessieeastburn/Documents/Fall 2021/GEOG 6000/datafiles/ne_50m_admin_1_states_provinces/ne_50m_admin_1_states_provinces.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 83 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.shp")
## Reading layer `ne_50m_lakes' from data source 
##   `/Users/jessieeastburn/Documents/Fall 2021/GEOG 6000/datafiles/ne_50m_lakes/ne_50m_lakes.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 275 features and 35 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.shp")
## Reading layer `ne_50m_populated_places' from data source 
##   `/Users/jessieeastburn/Documents/Fall 2021/GEOG 6000/datafiles/ne_50m_populated_places/ne_50m_populated_places.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1249 features and 119 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

Plot data and produce a map of ozone concentrations at the measurement stations

plot(st_geometry(ozone.sf), reset = TRUE)
plot(st_geometry(lakes), reset = TRUE, add = TRUE, col = "lightblue")
plot(st_geometry(states), reset = TRUE, add = TRUE)
plot(ozone.sf["ozone"], add = TRUE, pch = 16)

places <- cbind(places, st_coordinates(st_centroid(places)))

tmap_options(check.and.fix = TRUE)
mybbox <- st_bbox(ozone.sf)
tm_shape(lakes, bbox = mybbox) + 
  tm_fill("lightblue") + 
  tm_shape(states) + 
  tm_borders() +
  tm_shape(ozone.sf) + 
  tm_symbols(col = "ozone", size = 0.5, palette = "viridis") +
  tm_shape(places) + 
  tm_text("NAME", size = 0.75)

Produce 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

ozone.sf$ppb100 <- ozone.sf$ozone > 100

ozone.sf$lppb100 <- log(ozone.sf$ozone + 1e-1)
ozone.sf$logppb100 <- ozone.sf$ozone > 100
plot(ozone.sf["ppb100"], 
     pch = 16, 
     main = "Ozone > 100 ppb")
plot(states,
     add = TRUE, 
     pch = 16, 
     cex = 1.5,
     col = NA,
     fill = NA,
     border = "black")

I observe more ozone values over 100 ppb over Lake Michigan and the Chicago area.

ppb100.var <- variogram(logppb100 ~ 1, ozone.sf)

plot(ppb100.var)

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

modNugget <- 0.07
modRange <- 88
modSill <- 0.070
ppb100.vgm <- vgm(psill = modSill, 
                  model = "Gau", 
                  range = modRange, 
                  nugget = modNugget)
plot(ppb100.var, ppb100.vgm)

ppb100.vgm2 <- fit.variogram(ppb100.var, ppb100.vgm)
plot(ppb100.var, ppb100.vgm2)

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

pred.grid <- st_as_stars(mybbox, 
                         xlim = c(-94, -82), 
                         ylim = c(36, 45), 
                         dx = 0.1, 
                         dy = 0.1,
                         breaks = "quantile")

st_crs(pred.grid) <- 4326
ppb100.ik <- krige(ppb100 ~ 1, 
                  ozone.sf, 
                  pred.grid, 
                  ppb100.vgm2, 
                  nmax = 40)
## [using ordinary kriging]
st_crs (ppb100.ik) <- 4326
plot(ppb100.ik["var1.pred"])

ppb100.ik$var1.pred[which(ppb100.ik$var1.pred < 0)] <- 0
my.pal <- rev(viridis::magma(10))
newppb100.ik <- plot(ppb100.ik["var1.pred"], 
                       col = my.pal, 
                       breaks = seq(0, 1, length.out = 11),
                       main = "Ozone > 100ppb", reset = FALSE,)
plot(st_geometry(states), reset = FALSE, add = TRUE)
plot(st_geometry(places), reset = FALSE, add = TRUE)
plot(st_geometry(lakes), col = "lightblue", reset = FALSE, add = TRUE)

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

-In the interpolated map we can see high concentrations of ozone in the Chicago area and around Lake Michigan. I conducted some research on the high ozone levels around Lake Michigan and it is suspected the source of the ozone is from pollution being blown upwind from the urban city of Chicago.