ch14

14_kriging.Rmd

@@ -2,12 +2,223 @@
 
 **Learning objectives:**
 
-- THESE ARE NICE TO HAVE BUT NOT ABSOLUTELY NECESSARY
+- What is Kriging
+- How it can be used to interpolate spatial data
+- How to perform Kriging in R
 
-## SLIDE 1 {-}
 
-- ADD SLIDES AS SECTIONS (`##`).
-- TRY TO KEEP THEM RELATIVELY SLIDE-LIKE; THESE ARE NOTES, NOT THE BOOK ITSELF.
+
+
+## What is Kriging
+
+
+```{r warning=FALSE, message=FALSE, echo=FALSE}
+library(tidyverse)
+theme_set(theme_minimal())
+```
+
+Kriging is a geostatistical interpolation technique that is used to estimate the value of a variable at an unobserved location based on the values of the variable at nearby locations. Kriging is based on the assumption that the spatial correlation between the values of the variable decreases with distance. Kriging is widely used in to interpolate spatial data when data are sparse or irregularly distributed.
+
+
+## Types of Kriging
+
+For example, Simple Kriging assumes the mean of the random field, μ(s), is known; 
+
+- **Simple Kriging**: Assumes that the mean of the variable is known and constant across the study area.
+
+Formula: $Z(s_0) = \mu + \sum_{i=1}^{n} \lambda_i (Z(s_i) - \mu)$
+
+where $\mu$ is the mean, $\lambda_i$ are the weights, and $Z(s_i)$ are the observed values.
+
+
+Ordinary Kriging assumes a constant unknown mean, μ(s)=μ; 
+
+- **Ordinary Kriging**: Assumes that the mean of the variable is unknown and varies across the study area.
+
+Formula: $Z(s_0) = \sum_{i=1}^{n} \lambda_i Z(s_i)$
+
+where $\lambda_i$ are the weights and $Z(s_i)$ are the observed values.
+
+Universal Kriging can be used for data with an unknown non-stationary mean structure.
+
+- **Universal Kriging**: Assumes that the mean of the variable is unknown and varies across the study area, but can be modeled as a function of covariates.
+
+Formula: $Z(s_0) = \sum_{i=1}^{n} \lambda_i Z(s_i) + \beta X(s_0)$
+
+where $\lambda_i$ are the weights, $Z(s_i)$ are the observed values, $\beta$ is the coefficient for the covariate $X(s_0)$, and $X(s_0)$ is the value of the covariate at the unobserved location.
+
+## Performing Kriging in R
+
+To perform Kriging in R, you can use the `gstat` package, which provides functions for geostatistical analysis. The `gstat` package allows you to create a variogram model, fit the model to the data, and use the model to interpolate values at unobserved locations.
+
+```{r}
+library(gstat)
+```
+
+### Example: Simple Kriging
+
+In this example, we will perform simple Kriging on a dataset of air quality measurements. We will estimate the air quality at unobserved locations based on the measurements at nearby monitoring stations.
+
+Step 1: Load the air quality data and create a variogram model.
+
+```{r warning=FALSE, message=FALSE}
+library(sp)
+# Load the air quality data
+data(meuse)
+coordinates(meuse) <- c("x", "y")
+meuse%>%
+  as_data_frame()%>%
+  select(x,y,zinc)%>%
+  head()
+```
+
+Step 2: Visualize the data.
+```{r}
+meuse%>%
+  as_data_frame()%>%
+  select(x,y,zinc)%>%
+  ggplot(aes(x,y,color=zinc))+
+  geom_point()+
+  scale_color_viridis_c()
+```
+
+What we need is a grid of points where we want to predict the zinc concentration. 
+
+```{r}
+data(meuse.grid)
+meuse.grid %>%
+  as_data_frame()%>%
+  select(x,y)%>%
+  ggplot(aes(x,y))+
+  geom_point(shape=".")
+```
+
+```{r}
+zinc_data <- meuse%>%
+  as_data_frame()%>%
+  select(x,y,zinc)
+
+grid_data <- meuse.grid %>%
+  as_data_frame()%>%
+  select(x,y)
+
+ggplot()+
+  geom_point(data=grid_data,aes(x,y),shape=".")+
+  geom_point(data=zinc_data,aes(x,y,color=zinc))+
+  scale_color_viridis_c()
+```
+
+### Variogram Model
+
+We perform a Variogram analysis to understand the spatial correlation of the zinc concentration in the dataset.
+
+Formula: $V(h) = \frac{1}{2N(h)} \sum_{i=1}^{N(h)} (Z(s_i) - Z(s_i + h))^2$
+
+where $V(h)$ is the semivariance at lag distance $h$, $N(h)$ is the number of pairs of observations at lag distance $h$, $Z(s_i)$ is the value of the variable at location $s_i$, and $Z(s_i + h)$ is the value of the variable at location $s_i$ plus lag distance $h$.
+
+```{r}
+# Create a variogram model
+vc <- variogram(log(zinc) ~ 1, meuse, cloud = TRUE)
+plot(vc)
+v <- variogram(log(zinc) ~ 1, meuse)
+plot(v)
+```
+
+```{r}
+# Fit the variogram model
+fv <- fit.variogram(object = v, 
+                       model = vgm(psill = 0.5,
+                                   model = "Sph", 
+                                   range = 900,
+                                   nugget = 0.1))
+plot(v,fv)
+```
+
+### Perform Simple Kriging
+
+`gstat` function to compute the Kriging predictions:
+
+      `?gstat`
+
+```{r}
+library(sf)
+
+data(meuse)
+data(meuse.grid)
+
+meuse <- st_as_sf(meuse, 
+                  coords = c("x", "y"), 
+                  crs = 28992)
+
+meuse.grid <- st_as_sf(meuse.grid, 
+                       coords = c("x", "y"),
+                       crs = 28992)
+```
+
+```{r}
+v <- variogram(log(zinc) ~ 1, meuse)
+plot(v)
+```
+
+```{r}
+fv <- fit.variogram(object = v,
+                    model = vgm(psill = 0.5, 
+                                model = "Sph",
+                                range = 900, 
+                                nugget = 0.1))
+fv
+plot(v, fv, cex = 1.5)
+```
+
+
+```{r}
+k <- gstat(formula = log(zinc) ~ 1, 
+           data = meuse, 
+           model = fv)
+kpred <- predict(k,meuse.grid)
+```
+
+```{r}
+ggplot() + 
+  geom_sf(data = kpred, 
+          aes(color = var1.pred)) +
+ # geom_sf(data = meuse) +
+  viridis::scale_color_viridis(name = "log(zinc)")
+
+ggplot() + 
+  geom_sf(data = kpred, 
+          aes(color = var1.var)) +
+  geom_sf(data = meuse) +
+  viridis::scale_color_viridis(name = "variance") 
+```
+
+
+```{r}
+# Perform simple Kriging
+kriged <- krige(log(zinc) ~ 1, meuse,
+                kpred,
+                model = fv)
+```
+
+### Plotting the Results
+
+```{r}
+# Plot the results
+plot(kriged)
+```
+
+## Summary
+
+- Kriging is a geostatistical interpolation technique used to estimate the value of a variable at unobserved locations.
+- There are different types of Kriging methods, including simple Kriging, ordinary Kriging, universal Kriging, and co-Kriging.
+- In R, you can perform Kriging using the `gstat` package, which provides functions for geostatistical analysis.
+
+## Additional Resources
+
+- [gstat package documentation](https://cran.r-project.org/web/packages/gstat/gstat.pdf)
+- [Geostatistics with R](https://www.r-bloggers.com/2016/02/geostatistics-with-r-tutorial/)
+- [Introduction to Geostatistics](https://www.youtube.com/watch?v=8v9v7JcOwJc)
+
 
 ## Meeting Videos {-}
 

DESCRIPTION

@@ -19,12 +19,14 @@ Imports:
     nngeo,
     rmarkdown,
     sf,
+    sp,
     SpatialEpi,
     spData,
     spdep,
     terra,
     tibble,
     tidyterra,
-    tidyverse
+    tidyverse,
+    viridis
 Additional_repositories: https://inla.r-inla-download.org/R/stable
 Encoding: UTF-8

map.adj

@@ -0,0 +1,68 @@
+67
+1 3 21 28 67
+2 5 3 4 10 63 65
+3 6 2 10 16 32 33 65
+4 4 2 10 37 63
+5 5 7 11 29 31 56
+6 6 15 36 38 39 46 54
+7 5 5 11 14 17 31
+8 5 41 57 58 59 66
+9 4 39 46 48 51
+10 8 2 3 4 16 37 43 61 65
+11 7 5 7 14 17 32 56 65
+12 5 17 18 24 42 53
+13 5 39 40 45 48 54
+14 7 7 11 17 18 31 44 60
+15 4 6 23 36 46
+16 5 3 10 27 33 61
+17 8 7 11 12 14 18 24 32 33
+18 6 12 14 17 41 53 60
+19 6 40 41 47 49 54 57
+20 4 25 43 61 62
+21 5 1 22 28 50 67
+22 7 21 36 38 49 50 54 67
+23 3 15 46 51
+24 6 12 17 27 33 42 62
+25 2 20 62
+26 4 30 56 63 65
+27 6 16 24 33 42 61 62
+28 6 1 21 29 31 34 50
+29 3 5 28 31
+30 2 26 63
+31 7 5 7 14 28 29 34 44
+32 5 3 11 17 33 65
+33 6 3 16 17 24 27 32
+34 6 28 31 44 49 50 55
+35 5 40 45 58 64 66
+36 5 6 15 22 38 67
+37 3 4 10 43
+38 4 6 22 36 54
+39 6 6 9 13 46 48 54
+40 7 13 19 35 45 54 57 66
+41 9 8 18 19 47 49 53 57 59 60
+42 5 12 24 27 53 62
+43 4 10 20 37 61
+44 5 14 31 34 55 60
+45 6 13 35 40 48 52 64
+46 6 6 9 15 23 39 51
+47 3 19 41 49
+48 4 9 13 39 45
+49 9 19 22 34 41 47 50 54 55 60
+50 6 21 22 28 34 49 55
+51 3 9 23 46
+52 2 45 64
+53 5 12 18 41 42 59
+54 8 6 13 19 22 38 39 40 49
+55 5 34 44 49 50 60
+56 4 5 11 26 65
+57 5 8 19 40 41 66
+58 4 8 35 64 66
+59 3 8 41 53
+60 6 14 18 41 44 49 55
+61 6 10 16 20 27 43 62
+62 6 20 24 25 27 42 61
+63 5 2 4 26 30 65
+64 4 35 45 52 58
+65 8 2 3 10 11 26 32 56 63
+66 5 8 35 40 57 58
+67 4 1 21 22 36