R12_Semivariogram.Rmd

---
title: "Semivariogram example - Meuse river sediments"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(sp)
library(gstat)
data(meuse)
```

* This illustration uses the Meuse dataset included in gstat which comprises of heavy metals measured in the top soil in a flood plain along the river Meuse.

* Data set gives locations and topsoil heavy metal concentrations, along with a number of soil and landscape variables at the observation locations, collected in a flood plain of the river Meuse, near the village of Stein (NL). 

* Heavy metal concentrations are from composite samples of an area of approximately 15 m x 15 m. (recorded sample location is a single point, yet each actual sample may cover approx. 15x15 m). 

* Total area covered by sampling fits into a rectangle 2.8 km (x-axis) by 3.9 km (y-axis) approximately.

* Apparently, polluted sediment is carried by the river and mostly deposited close to the river bank.

```{r}
str(meuse)
```

--- 

## Coordinates

Function `coordinates()`: promotes the data.frame meuse into a SpatialPointsDataFrame which knows about its spatial coordinates.

* `x` and `y` coordinates are in meters.

    + `x` a numeric vector; Easting (m) in Rijksdriehoek (RDH) (Netherlands topographical) map coordinates  
    + `y` a numeric vector; Northing (m) in RDH coordinates



```{r}
coordinates(meuse) = ~x+y
class(meuse)
str(meuse)
head(meuse)
```

The function `coordinates()` can retrieve spatial coordinates from a SpatialPointsDataFrame: 

```{r}
coordinates(meuse)[5:15,]
```

Now, we can plot the data using `spplot` and `bubble` as illustrated below (note: the x- and y-axis are the spatial coordinates)


```{r}
spplot(meuse, "zinc",  colorkey = TRUE, main = "zinc concentrations (ppm)")
bubble(meuse, "zinc", col=c("#00ff0088"), main = "zinc concentrations (ppm)")
```

Concentrations seem to follow (inversely) from distances to the river

```{r}
data(meuse.grid) # we read-in a separately provided sample grid (map)
coordinates(meuse.grid) = ~x+y
gridded(meuse.grid) = TRUE
class(meuse.grid)
# plot distances to river
image(meuse.grid["dist"], main="distance to river (red = 0)")
```

---- 

## Semivariogram

* Calculate the sample variogram. This is done with the `variogram()` function.

* Fit a model to the sample variogram using `fit.variogram()` function

The `variogram()` function takes two arguments: the first denotes how one or more variables interact spatially, and the second is an SPDF where those variables reside. 

* In our example, we assume that there is a constant mean (no trend) for the variable `log(zinc)`

```{r}
lzn.vgm <- variogram(log(zinc)~1, meuse) # calculates sample variogram values 
lzn.vgm # 15 groups | # obs per group | repres distances | gamma | ...isotropy assumed
```


For the `fit.variogram()` function, a sample variogram is the first argument. The second is the model, with parameters, to be fit to the sample variogram. For a list of all possible variograms that can be used, call `?vgm`, and to see graphical properties/characteristics of these models, call  `?show.vgms`.

* `vgm` arguments (quick summary)   
    + `psill=NA` - starting value for the sill argument is not provided (generated internally for the estimation algorithm)  
    + `Sph` - model type, e.g. "Exp", "Sph", "Gau"  
    + `range=900` range in meters  
    + `nugget=1` nonzero nugget added (included in variogram calculation) 

```{r}
lzn.fit <- fit.variogram(lzn.vgm, model=vgm(psill= NA, "Sph", range=900, nugget=1)) # fit model
lzn.fit
plot(lzn.vgm, lzn.fit) # plot the fitted variogram and the observed variogram on the same graph
```

----- 

## Fitted semivariograms depend on assumptions used:

**Gaussian** semivariance wrt distance:

```{r}
lzn.gau <- fit.variogram(lzn.vgm, model=vgm(psill= NA, "Gau", range=900, nugget=1)) # fit model
#lzn.gau
plot(lzn.vgm, lzn.gau) 
```

**Exponential** semivariance wrt distance:

```{r}
lzn.exp <- fit.variogram(lzn.vgm, model=vgm(psill= NA, "Exp", range=900, nugget=1)) # fit model
#lzn.exp
plot(lzn.vgm, lzn.exp) 
```


------ 

**Directional (anisotropic)** semivariance wrt distance:

```{r}
vgm3.dir <- variogram(log(zinc)~1, meuse, alpha = c(0, 45, 90, 135))
vgm3.fit.d <- fit.variogram(vgm3.dir, model = vgm(.59, "Sph", 1200, .05, anis = c(45, .4)))
#vgm3.fit.d
plot(vgm3.dir, vgm3.fit.d, as.table = T)
```


------

For spatio-temporal semivariogram, see e.g. [tutorial here](http://r-video-tutorial.blogspot.com/2015/08/spatio-temporal-kriging-in-r.html)