library(ggplot2)
library(gstat)
library(RColorBrewer)
library(sf)
library(stars)
library(viridis)
Read in data
ozone.sf <- st_read("ozone.shp")
st_crs(ozone.sf) <- 4326

states <- st_read("ne_50m_admin_1_states_provinces_shp.shp")
lakes <- st_read("ne_50m_lakes.shp")
places <- st_read("ne_50m_populated_places.shp")
places <- cbind(places, st_coordinates(st_centroid(places)))
## Warning in st_centroid.sf(places): st_centroid assumes attributes are constant
## over geometries of x
## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for longitude/
## latitude data
Map of ozone concentrations at measurment stations
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(X, Y, label = NAME), size = 2.5) +
  coord_sf(xlim = c(-94, -82), ylim = c(36, 45), expand = FALSE) +
  theme_bw()

Histogram plot
hist(ozone.sf$ozone)

Plotof Ozone greater than 100ppb
#Create new variable
ozone.sf$oz100 <- ozone.sf$ozone > 100
plot(ozone.sf["oz100"], 
     pch = 16, 
       main = "Ozone > 100ppb")

Create sample variogram

For the sample variogram and fitted model for ozone concentrations over 100ppb, a spherical model was used with sill = 0.15, range = 100, and nugget = 0.01. At smaller distances, the semivariogram model rises sharply, then levels off. The range where it flattens out indicates little auto correlation in ozone values.

ppt100.var <- variogram(oz100 ~ 1, ozone.sf)

plot(ppt100.var, plot.numbers =TRUE)

ppt100.vgm <- vgm(0.15, "Sph", 100, 0.01)

plot(ppt100.var, ppt100.vgm)

Fit variogram model
ppt100.vgm2 <- fit.variogram(ppt100.var, ppt100.vgm)

plot(ppt100.var, ppt100.vgm2)

Interpolate using oridnary kriging
#First create a grid
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

ppt100.ik <- krige(oz100 ~ 1, 
                  ozone.sf, 
                  pred.grid, 
                  ppt100.vgm2, 
                  nmax = 40)
## [using ordinary kriging]
# Plot predictions 
plot(ppt100.ik["var1.pred"])

#ppt100.ik$var1.pred[which(ppt100.ik$var1.pred < 0)] <- 0

my.pal <- rev(viridis::magma(10))

plot(ppt100.ik["var1.pred"], 
     col = my.pal, 
     breaks = seq(0, 1, length.out = 11),
     main = "Indicator Kriging Ozone > 100ppbs", reset = FALSE,)
plot(st_geometry(states), reset = FALSE, add = TRUE)
plot(st_geometry(places), size = 2.5, reset = FALSE, add = TRUE)

The Northern region and the small regions in the east and east-central on the map indicate high probabilities of ozone concentrations over 100ppb. These regions fall within metro areas of the Midwest where there’s heavy traffic and likely industrial sites. The southern region of the map indicates having the lowest probability of ozone concentrations of over 100ppbs.