setwd("~/Documents/geog6000/lab13")

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.6.2
library(gstat)
## Warning: package 'gstat' was built under R version 3.6.2
library(RColorBrewer)
library(sf)
## Warning: package 'sf' was built under R version 3.6.2
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(stars)
## Loading required package: abind
library(viridis)
## Warning: package 'viridis' was built under R version 3.6.2
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 3.6.2

Read in data

ozone.sf <- st_read("../datafiles/mwozone/ozone.shp")
## Reading layer `ozone' from data source 
##   `/Users/serenamadsen/Documents/geog6000/datafiles/mwozone/ozone.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 147 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -93.572 ymin: 36.791 xmax: -82.96 ymax: 44.453
## CRS:           NA
st_crs(ozone.sf) <- 4326
## Read data
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/serenamadsen/Documents/geog6000/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
## CRS:           4326
lakes <- st_read("../datafiles/ne_50m_lakes/ne_50m_lakes.shp")
## Reading layer `ne_50m_lakes' from data source 
##   `/Users/serenamadsen/Documents/geog6000/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
## CRS:           4326
places <- st_read("../datafiles/ne_50m_populated_places/ne_50m_populated_places.shp")
## Reading layer `ne_50m_populated_places' from data source 
##   `/Users/serenamadsen/Documents/geog6000/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
## CRS:           4326

A map of ozone concentrations at the measurement stations

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

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

plot(ozone.sf["ppb100"], 
     pch = 16, 
       main = "Ozone > 100ppb")

The yellow points indicate high concentration of ozone (above 100 ppb). Many of these points are concentrated in the top center (My best guess is SW Wisconcin and NE Illinios).

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(ppb100 ~ 1, ozone.sf)

plot(ozone100.var)

modNugget <- 0.065
modRange <- 90
modSill <- 0.085

ppb100.vgm <- vgm(psill = modSill, 
                  model = "Cir", 
                  range = modRange, 
                  nugget = modNugget)
plot(ozone100.var, ppb100.vgm, main = "Ozone > 100ppb Variogram")

Fitted variogram:

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

plot(ozone100.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

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
ppb100.ok <- krige(ppb100 ~ 1, 
                  ozone.sf, 
                  pred.grid,
                  ppb100.vgm2, 
                  nmax = 40)
## [using ordinary kriging]
plot(ppb100.ok["var1.pred"])

ppb100.ok$var1.pred[which(ppb100.ok$var1.pred < 0)] <- 0

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

colorppb100.ok <- plot(ppb100.ok["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? The spatial pattern we see here is higher ozone in and north of Chicago. I don’t know much about pollutants in that area, so I did some research to make an educated guess as to why ozone pollution is so high. Some main contributing factors are onshore winds that blow pollution into the area, vehicle emissions, industry, and coal-fired power plants. (https://www.cleanwateractioncouncil.org/issues/air-quality/ozone-pollution/)

The lake polygons cover up important information on this map. With a bit more knowledge, I would make the lakes somewhat transparent to avoid this, but I felt that the lakes were crucial to understanding the reasoning behind such high ozone levels.