---
title: "Demonstrate Anisotropy"
author: "D G Rossiter"
date: "April 17, 2024"
output: html_document
---

This shows how anisotropy is detected, modelled, and used for interpolation in Ordinary Kriging, using the Meuse dataset as an example. It uses the older `sp` spatial representation; this could be updated to use `sf`.

Load the dataset, log-transform the Zn concentration, display it as a post-plot:

```{r}
library(gstat)
library(sp)
data(meuse)
coordinates(meuse) <- c("x", "y")
meuse$logZn <- log(meuse$zinc)
spplot(meuse, zcol="logZn")
```


The spatial dependence looks longer-range along the NNE-SSW  axis (parallel to most of the river course) and short-range at right angles to this (moving away from the river).

Confirm with a *variogram map*. This shows directional empircal variograms from the centre; notice the axes are marked `dx` and `dy`.

```{r}
v.map <- variogram(logZn ~ 1, meuse, cutoff=1300, width=90, map=TRUE)
plot(v.map)
```

The major axis appears to be at about 30 deg azimuth from N.
 
Compute and display an anisotropic variogram. The new arguments to `variogram` specify the directions `alpha` and angle tolerance `tol.hor` (angle on each side of the major axis). We choose for a fairly narrow angle to get a good idea of the anisotropy, however, this limits the number of point-pairs.


```{r}
(v.a <- variogram(logZn ~ 1, meuse, cutoff=1300, width=90, alpha=c(30,120), tol.hor=25))
```


Note how the `dir.hor` field is now filled with either 30 or 120 (the two azimuths that were requested), and the average separations in each bin (field `dist`) are different for the two azimuths, because of the different point-pair sets in each bin.

Plot the variogram:

```{r}
print(plot(v.a, plot.numbers=T))
```

Note that there is now one variogram for each direction. The variogram on the long axis (1) has more point-pairs per bin; (2) reaches a definite sill around 1200 m range. The variogram on the short axis is more difficult: (1) it has fewer point pairs, so harder to model; (2) it increases without a definite sill, but note the very small number of pairs after about 600 m.

The simplest kind of anisotropy is *geometric*: same model, same sill, but different ranges. Another kind is *zonal*, which is quite difficult to model. Here we have a small observation set, and the highest kriging weights are at short range anyway, so despite the unbounded variogram at azimuth 120 we choose to fit a geometric anisotropic model.

Estimate a geometric anisotropy model to the anisotropic empirical variogram. The new information is the anisotropy parameter `anis`; this has an azimuth (degrees from N) and an *anisotropy ratio*, that is, the ratio of ranges of the short to the long axis. For example, if the range of the long axis is 850 m, and the ratio is 0.5, the range of the short axis is 425 m.  It must be estimated by eye. Here I estimate the long range as 1000 m and the short as 500 m, so the anisotropy ratio is 0.5.


```{r}
vm.a <- vgm(psill=0.5,model="Sph",range=1000,nugget=0.005, anis=c(30, 0.5))
print(plot(v.a, pl=T, model=vm.a))
```

Adjust the fit with `fit.variogram`. Note that the anisotropy ratio is _not_ adjusted by `fit.variogram`, although the range itself is.

```{r}

(vmf.a <- fit.variogram(v.a, vm.a))
print(plot(v.a, pl=T, model=vmf.a))
```

The fit reduced the range a bit.

Predict over the grid with the anisotropic model:

```{r}
data(meuse.grid)
coordinates(meuse.grid) <- ~ x + y
k40.a <- krige(logZn ~ 1, locations=meuse, newdata=meuse.grid, model=vmf.a)
```

Display the predictions:

```{r}
print(spplot(k40.a, "var1.pred", asp=1, col.regions=bpy.colors(64),
             main="anisotropic OK prediction, log-ppm Zn"))
```

The effect of the anisotropic model is obvious; similar predictions are "stretched" along the major axis.

Display the prediction variances:

```{r}
print(spplot(k40.a, "var1.var",
             col.regions=cm.colors(64),
             asp=1,
             main="anisotropic OK prediction variance, log-ppm Zn^2"))
```

For comparison, compute the omnidirectional variogram, fit a model to it, and use this to predict over the grid:

```{r}
v <- variogram(logZn ~ 1, meuse, cutoff=1300, width=90)
vm <- vgm(psill=0.12,model="Sph",range=850,nugget=0.01)
(vmf <- fit.variogram(v, vm))
print(plot(v, pl=T, model=vmf))
k40 <- krige(logZn ~ 1, locations=meuse, newdata=meuse.grid, model=vmf)
print(spplot(k40, "var1.pred", asp=1, col.regions=bpy.colors(64),
             main="OK prediction, log-ppm Zn"))
print(spplot(k40, "var1.var",
             col.regions=cm.colors(64),
             asp=1,
             main="OK prediction variance, log-ppm Zn^2"))
```


Compute the difference between the anisotropic and isotropic predictions and their variances; add new fields to the kriging grid, and display these differences as maps.

```{r}
k40$pred.diff <- k40.a$var1.pred - k40$var1.pred
summary(k40$pred.diff)
print(spplot(k40, "pred.diff", asp=1, col.regions=topo.colors(64),
             main="difference aniso- vs. iso-tropic prediction, log-ppm Zn"))

k40$var.diff <- k40.a$var1.var - k40$var1.var
summary(k40$var.diff)
print(spplot(k40, "var.diff", asp=1, col.regions=terrain.colors(64),
             main="difference aniso- vs. iso-tropic prediction variances, log-ppm Zn^2"))
```

The effect of taking anisotropy of spatial dependence into account is clear.

This example is somewhat contrived, because the anisotropy is here an artefact of the main course of the river. It does not apply to the southern section of the area. A  better approach is to use distance from river as a covariable in KED.

But anyway we have demonstrated how to detect, specify, and predict using anisotropic OK.