---
title: "Geostatistics with the R software"
output: rmdformats::readthedown
---

![Fig.: Kriging French Temperatures (as measured by Meteo France August, 31th at 6AM)](https://mribatet.perso.math.cnrs.fr/CentraleNantes/Data/KrigingFrenchTemperatures.png)

## Manipulating spatial objects with R

In this lab we will work with the *sf* library and the *tidyverse*. We thus need to load them (and even install them if not already):

```{r}
library(sf)
library(tidyverse)
```

We are now all set to load a dataset provided with the *sf* library.

```{r}
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"))
plot(nc["SID74"])
```

-   What kind of data it is?

**Write answer here**

Let's have a closer look to the data set.

```{r}
class(nc)
nc
```

From the above output, we can see that it is an object of class *nf* with 100 features (i.e., measurements) and 14 fields (i.e., variables). The geometry of this dataset is a multipolygon and the geophysical dimension are *XY*.

```         
+ Lastly, we does it mean that the Geodetic CRS is set to *NAD27*?

**Write answer here**

+ Run the following code (BIR79 is the number of birth per county for year 1979). Is the plot meaningfull? What was the aim of setting "EPSG:32119"?
```

```{r}
nc2 <- st_transform(nc, "EPSG:32119")
plot(nc2["BIR79"])
```

Polygons are one type of geometry and you can have points, lines, polygons, pixels as well.

If you want to show off a bit you can even get an interactive map by invoking:

```{r}
library(mapview)##install this library first if needed
mapview(nc["AREA"])
```

-   From this interactive plot did you catch why our dataset was named 'nc'?

**Write answer here**

-   Just for fun... Run the following piece of code and guess what we are plotting here.

```{r}
url <- "https://mribatet.perso.math.cnrs.fr/CentraleNantes/Data/track.RData"
if (!file.exists("./track.RData")){
  download.file(url, destfile = "./track.RData")
}
load("./track.RData")
mapview(track)
```

## Empirical variogram

We now change our dataset for $NO_2$ concentration levels in Germany. The dataset is easily retrieved using

```{r}
no2 <- read_csv(system.file("external/no2.csv", package = "gstat"))
```

-   What is the type of the *no2* object? Where are stored the spatial coordinates in this dataset?

**Write answer here**

We thus need to 'convert' it so that one it works easily with the *sp* library:

```{r}
no2.sf <- st_as_sf(no2, crs = "OGC:CRS84",
                   coords = c("station_longitude_deg", "station_latitude_deg")) |>
  st_transform(crs = "EPSG:32632")
```

-   If you tried to plot the no2.sf package, you probably noticed that there were no borders. One could use *mapview* as before to get an interactive map but for static one, and another reason that we will see later, it is good practice to retrieve the boundary of our study region (here Germany). There are plenty of ways to do it, one way it to make use of the *geodata* library. Make sure you understand the following piece of code:

```{r}
library(geodata)
DE <- gadm(country = "DE", level = 0, path = tempdir())
```

-   What is the class of the *DE* object? What is the meaning of *level = 0*?

**Write your answer here.**

-   *SpatVector* is not a supported class of the *sp* library so we need to convert it to the appropriate class. This is done with one simple line of code:

```{r}
DE <- st_as_sf(DE)
```

-   Try to retrieve the boundary of your favorite "région Française"?

```{r}
## Insert code here
```

-   Produce a plot of the study region with weather stations of our $NO_2$ dataset. Make sure that the *crs* are the same!

```{r}
## Insert code here
```

-   Read the documentation of the *variogram* function of the *gstat* library and compute the empirical variogram assuming a (unknown) constant mean. What about is the mean is supposed to change linearly on altitude?

```{r}
## Insert code here
```

-   Play with the *vgm* function of the *gstat* library to see how change the variogram with the family and the parameters. Why one parameter is called *partial sill*?

```{r}
## Insert code here
```

-   Use the *fit.variogram* function (again of the *gstat* library) to fit a parametric variogram to the above empirical variogram from (generalized) least-squares.

```{r}
## Insert code here
```

-   Use the *model* argument to compare the empirical variogram to the fitted one.

```{r}
## Insert code here
```

## Kriging

We are now close to be able to do our first spatial interpolation using kriging. Typically we will perform predictions on a grid within our study region. In our case, a grid over Germany. Fortunately, the *stars* library will help us. Let's see how to do this.

```{r}
library(stars)##install it if needed first
grid <- st_bbox(DE) |> st_as_stars(dx = 0.05) |> st_crop(DE)
```

-   Use the *krige* function of the *gstat* library to perform ordinary kriging prediction on our brand new grid!

```{r}
library(gstat)
## Write code here
```

-   Plot the prediction map as well as the kriging standard error.

```{r}
## Write code here
```

## Conditional simulations

In the lecture we saw that the kriging estimator may provide unrealistically to smooth predictions and that, one way to bypass this hurdle is to get prediction from conditional simulations.

Conditional simulation can be done from the *gstat* library using the *krige* function by passing an extra argument *nsim* (the number of conditional simulation). You might want to set the *nmax* argument to a reasonably small value to lighten the CPU load.

-   Perform a bunch of conditional simulations and plot some of them.

```{r}
## Write code here
```

## Case study: French temperatures

Go to the [ECAD website](https://www.ecad.eu) and download mean temperatures over France. Provide kriging estimator for the annual mean temperature for some years of your choice.
