## Exercise  1

* use the dataset found in `system.file("shape/nc.shp", package="sf")`
* Read this dataset with `sf::st_read` (or `sf::read_sf`)

```{r}
library(sf)
nc = read_sf(system.file("shape/nc.shp", package="sf"))
```

* plot the dataset, then plot the second attribute
```{r}
plot(nc)
plot(nc[,2])
# or, e.g. 
# library(tidyverse)
# nc %>% select(2) %>% plot
```

* try to manipulate the color scheme and color breaks

```{r}
plot(nc[,2], pal = viridis::viridis) # note that this is a function, not a set of colors
plot(nc[,2], nbreaks = 5, pal = viridis::viridis)  # note that we don't get exactly 5, because of "pretty" style
```

* plot the first 20 features of the dataset

```{r}
plot(nc[1:20,2])
```

* create a plot with axes

```{r}
plot(nc[1:20,2], axes = TRUE)
```

* what is the coordinate reference system of the dataset?

```{r}
st_crs(nc)
```

* For the `nc` dataset, compute the area of the polygons with `st_area`
* compare the areas with the `AREA` field in the dataset; is there a discrepancy?
```{r}
a = st_area(nc)
plot(a, nc$AREA)
```

* select the polygon that intersects with `POINT (-81.49826 36.4314)`
```{r}
pt = st_as_sfc("POINT (-81.49826 36.4314)", crs = st_crs(nc))
nc[pt,]
```

## Exercises (2)

* Read the file `system.file("tif/L7_ETMs.tif", package = "stars")` with `stars::read_stars`
* print the dataset, and describe what it contains
* plot the dataset; how are the color breaks chosen?
```{r}
library(stars)
r = read_stars(system.file("tif/L7_ETMs.tif", package = "stars"))
r
plot(r)
```

The color breaks are chosen by the default method for `breaks`, which is "quantile"; this
means that (regularly chosen) quantiles of the pixel values are chosen; this means that the
map shows equal amounts of each color.

* find the maximum pixel value for each of the bands (`st_apply`)

```{r}
st_apply(r, 3, max)
```
* find and plot the maximum reflectance over all bands for all pixels
```{r}
plot(st_apply(r, 1:2, max))
```
* find and plot the mean reflectance over all bands for all pixels
```{r}
plot(st_apply(r, 1:2, mean))
```
* plot the subset of these data with the first 10 rows, the second 10 colums, and the last 3 bands
```{r}
r[,1:10, 11:20, 4:6]
```
* use `st_crop` to crop the dataset to the area 300 m around `POINT (293749.5 9115745)`
```{r}
pt = st_as_sfc("POINT (293749.5 9115745)", crs = st_crs(r))
cr = st_crop(r, st_buffer(pt, units::set_units(300, m)))
plot(cr)
```
* is this cropped dataset still a (rectangular) raster?
	
yes: the values outside the circle are simply set to `NA`, and are not plotted.

## Exercise 3

```{r}
data(DE_RB_2005, package = "gstat")
library(spacetime)
library(stars)
s = st_as_stars(as(DE_RB_2005, "STFDF"))
# select a single day, omit NA values:
s1 = na.omit(st_as_sf(s[,,100]))
names(s1)[1] = "PM10"
library(gstat)
variogram(PM10~1, s1)
plot(variogram(PM10~1, s1))
```

Building on the above dataset, create out plot where:

* setting the `cutoff` to a much larger, and a much smaller value

```{r}
plot(variogram(PM10~1, s1, cutoff = 1e5))
```

* setting the `width` argument to a larger and smaller value

```{r}
plot(variogram(PM10~1, s1, width =  10000, cutoff = 5e5))
plot(variogram(PM10~1, s1, width = 100000, cutoff = 5e5))
```

* the number of point pairs involved is also shown in the plot

```{r}
plot(variogram(PM10~1, s1, width = 100000, cutoff = 5e5, plot.numbers = TRUE))
```

* setting `cloud=TRUE` - what do the points in this plot represent? (hint: plot the number of point pairs)

```{r}
plot(variogram(PM10~1, s1, cutoff = 5e5, cloud = TRUE)) # width ignored
```

* Save a variogram you like to an object, and fit a model to it using `fit.variogram`; use `vgm()` to list the possible models and `show.vgms()` to get an idea how they look.

```{r}
v = variogram(PM10~1, s1)
init = vgm(20, "Sph", 300000, 5)
f = fit.variogram(v, init)
plot(v, f)
```

* If you get errors on fit, don't worry: 
    * plot the variogram with fitted model (like: `plot(v, fit)`, and examine how it looks
	* try to choose better initial values for nugget (offset), sill (final level), range (distance where correlation no longer increases)


## Exercises 4

Compare the ranges of the kriging and idw map, and the range of the PM10 data. Explain differences and correspondences.

```{r}
# load German boundaries
data(air, package = "spacetime")
de <- st_transform(st_as_sf(DE_NUTS1), 32632)
bb = st_bbox(de)
dx = seq(bb[1], bb[3], 10000)
dy = seq(bb[4], bb[2], -10000) # decreases!
st_as_stars(matrix(0., length(dx), length(dy))) %>%
  st_set_dimensions(1, dx) %>%
  st_set_dimensions(2, dy) %>%
  st_set_crs(32632) -> grd

i = idw(PM10~1, s1, grd)
v = variogram(PM10~1, s1)
f = fit.variogram(v, vgm(20, "Lin", 0, 4))
kr = krige(PM10~1, s1, grd, f)
range(i[["var1.pred"]])
range(kr[["var1.pred"]])
range(s1$PM10)
```
Kriging smoothes much more, mostly because the variogram model used (here) has a considerable relative nugget (nugget component as fraction of the total variance)

Examine the kriging variance error map.

```{r}
plot(kr[2])
```

Compute a standard error map by taking the square root of the kriging variance

```{r}
kr$stderr = sqrt(kr$var1.var)
```

Plot this map with the data points overlayed: can you explain the pattern in the errors?


```{r}
plot(kr["stderr"], reset = FALSE)
plot(st_geometry(de), add = TRUE, col = NA, border = 'red')
plot(s1, add = TRUE, col = 'green', pch = 3)
```

Repeat the exercise with the same linear model but with a zero nugget effect, and compare the ranges again.

```{r}
# brute-force: setting the nugget to zero:
f$psill[1] = 0
# somewhat "nicer": fitting the same model, but without a nugget:
f = fit.variogram(v, vgm(20, "Lin", 0))
kr0 = krige(PM10~1, s1, grd, f)
range(kr0[[1]])
```

Repeat the exercise with the a spherical model with a short range and zero nugget, and compare the ranges again.

```{r}
# use a short range spherical model, range 100 km, zero nugget:
f = vgm(20, "Sph", 100000)
kr.sph = krige(PM10~1, s1, grd, f)
range(kr.sph[[1]])
```

Repeat the exercise with the a pure nugget effect model (e.g., vgm(20, "Nug", 0)), and compare the ranges again.

```{r}
# fit a pure nugget model:
(f = fit.variogram(v, vgm(1, "Nug", 0)))
kr.nug = krige(PM10~1, s1, grd, f)
range(kr.nug[["var1.pred"]])
mean(s1$PM10)
range(kr.nug[["var1.var"]])
var(s1$PM10)
mean(variogram(PM10~1, s1, cloud=TRUE)$gamma)
mean(variogram(PM10~1, s1, cloud=TRUE, cutoff = 1e6)$gamma)
```

Compute (approximate) 95% prediction intervals by adding +/- 2 standard errors to the kriging


```{r}
kr$stderr = sqrt(kr$var1.var)
kr$upper = kr$var1.pred + 2 * kr$stderr
kr$lower = kr$var1.pred - 2 * kr$stderr
plot(merge(kr[c("lower", "upper")]))
```

Create a map with locations where the threshold of 25 has been exceeded, or might have been exceeded (is inside the prediction interval)


```{r}
kr$classify = ifelse(kr$lower > 25, 1, ifelse(kr$upper < 25, 3, 2))
plot(kr["classify"]) # legend 1: above 25, 2: ?, 3: below 25
```

## Exercises 5

* Carry out a block kriging, setting `block=c(10000,10000)` to estimate 10 km x 10 km block mean values
* Compare the kriged values to those obtained by point kriging, compare the ranges
* Compare the standard errors between point and block kriging

```{r}
f = fit.variogram(v, vgm(20, "Lin", 0, 4))
kr = krige(PM10~1, s1, grd, f) # point kriging
bl = krige(PM10~1, s1, grd, f, block = c(10000, 10000)) # block kriging
kr$point = kr$var1.pred
kr$block = bl$var1.pred
plot(merge(kr[c("point", "block")]))
kr$point.se = sqrt(kr$var1.var)
kr$block.se = sqrt(bl$var1.var)
plot(merge(kr[c("point.se", "block.se")]))
range(kr$point)
range(kr$block)
range(kr$point.se)
range(kr$block.se)
```



