## ----eval=FALSE---------------------------------------------------------------------------------------------------------
# l.p <- c("terra", "sf", "fields", "ranger", "ggplot2") # optional: "mapview"
# install.packages(l.p)


## -----------------------------------------------------------------------------------------------------------------------
#| label: info_on_dem
library(terra)
library(sf)
library(fields)

# set working directory to folder above "data"
r.dem <- rast("datasets/lareunion/data/terrain/reunion_30m.tif")
r.dem


## -----------------------------------------------------------------------------------------------------------------------
#| label: display_dem
plot(r.dem)


## -----------------------------------------------------------------------------------------------------------------------
#| label: prepare_data
# create a list of file names
l.rast <- paste0("datasets/lareunion/data/terrain/", c("reunion_30m.tif", "slope100m-perc.tif", "tpi_900m.tif"))

# create a raster stack
r.stack <- rast(l.rast)

# create a data frame, each row corresponds to a pixel
d.pixels <- as.data.frame(r.stack, xy = T, na.rm = F)
str(d.pixels)


## -----------------------------------------------------------------------------------------------------------------------
# only use every second pixel, ensures 60 m distance between sampling points
t.sel.x <- d.pixels$x %in% seq(min(d.pixels$x), max(d.pixels$x), by = res(r.stack)[1]*2)
t.sel.y <- d.pixels$y %in% seq(min(d.pixels$y), max(d.pixels$y), by = res(r.stack)[1]*2)

# remove NA, and remove every second pixel
d.pixels <- d.pixels[ complete.cases(d.pixels) & t.sel.x & t.sel.y, ]

str(d.pixels)


## -----------------------------------------------------------------------------------------------------------------------
#| label: clustering

# scale variables and round to 5 digits to make clustering more stable
d.pixels.sc <- round(scale(d.pixels), 5)

# set seed for reproducibility 
set.seed(11)
# k-means clustering, with centers = numbers of samples to plan
m.kmeans <- kmeans(d.pixels.sc, centers = 150, iter.max = 50)

# number of pixels per cluster
m.kmeans$size


## -----------------------------------------------------------------------------------------------------------------------
#| label: extract_clustercenters

# extract real points closest to cluster center
m.dist <- rdist(d.pixels.sc, m.kmeans$centers)
# select row indices 
idx <- apply(m.dist, 2, which.min)

d.sample <- d.pixels[idx, 1:2]

# create sample IDs
d.sample$id <- paste0("train_", 1:nrow(d.sample))

# create a spatial points object and save to Geopackage
v.sample <- st_as_sf(d.sample, coords = c("x", "y"), crs = 32740)
write_sf(v.sample, "datasets/lareunion/results/fsc_design_150.gpkg")

v.sample


## -----------------------------------------------------------------------------------------------------------------------
library(mapview)
mapview(v.sample)


## -----------------------------------------------------------------------------------------------------------------------
# extract cluster for each pixel
d.pixels$clusters <- fitted(m.kmeans, "classes")
# create a raster object 
r.clusters <- rast(d.pixels[, c("x", "y", "clusters")], type="xyz", crs= "epsg:32740")
writeRaster(r.clusters, filename = "datasets/lareunion/results/fsc_design_150_clusters.tif", overwrite = T)
plot(r.clusters)


## -----------------------------------------------------------------------------------------------------------------------
# use elevation as strata, from 3 strata
t.breaks <- quantile(d.pixels$reunion_30m, probs = seq(0,1,1/3))
d.pixels$strata <- cut(d.pixels$reunion_30m, breaks = t.breaks, labels = 1:3)

# check if it worked
table(d.pixels$strata)


## -----------------------------------------------------------------------------------------------------------------------
# save strata map
r.strata <- classify(r.stack[[1]], rcl = t.breaks)
writeRaster(r.strata, filename = "datasets/lareunion/results/strs_design_3_strata.tif", overwrite = T)

plot(r.strata)


## -----------------------------------------------------------------------------------------------------------------------
set.seed(1)

# apply random selection of 10 points for each strata
l.sample.rows <- lapply(c("1", "2", "3"), function(stratum){
  # select relevant row names 
  l.row <- rownames(d.pixels[ d.pixels$strata == stratum, ])
  # randomly select 10 of the row names
  l.rs <- sample(l.row, 10)
  return(l.rs)
})

# select corresponding pixel 
d.sample <- d.pixels[ unlist(l.sample.rows), ]

# check if number per strata are ok
table(d.sample$strata)

# create sample IDs
d.sample$id <- paste0("val_", 1:nrow(d.sample))

v.validation <- st_as_sf(d.sample[, c(1:2,5)], coords = c("x", "y"), crs = 32740)
write_sf(v.validation, "datasets/lareunion/results/strs_design_30.gpkg")


## -----------------------------------------------------------------------------------------------------------------------
mapview(v.validation)


## -----------------------------------------------------------------------------------------------------------------------
r.soc <- rast("datasets/lareunion/data/soil/soilgrids_SOCstock_0-30_mean30mSmooth.tif")
plot(r.soc)


## -----------------------------------------------------------------------------------------------------------------------
library(ranger)
library(ggplot2)

l.rast <- list.files("datasets/lareunion/data/terrain/", pattern = ".tif$", full.names = T)
r.stack <- rast(l.rast)

v.samples <- st_read("datasets/lareunion/results/fsc_design_150.gpkg")

# read response value 
# (in a real project, we would read the values from the lab measured on soil samples)
d.SOCstock <- extract(r.soc, v.samples)

str(d.SOCstock)


## -----------------------------------------------------------------------------------------------------------------------
# read predictor values at sampling locations 
d.predictors <- extract(r.stack, v.samples)

str(d.predictors)


## -----------------------------------------------------------------------------------------------------------------------
## fit RF model
m.rf <- ranger(y = d.SOCstock$`soilgrids_SOCstock_0-30_mean30mSmooth`,
               x = d.predictors[,-c(1)], # do not use ID
               importance ="permutation")
m.rf


## ----fig.width=6, fig.height=8------------------------------------------------------------------------------------------
t.importance <- importance(m.rf)

par(oma = c(2,10,2,2))
barplot(t.importance[order(t.importance)], horiz= T, las = 1, cex.names = 0.6)


## -----------------------------------------------------------------------------------------------------------------------
## Compute validation stats 
v.validation <- st_read("datasets/lareunion/results/strs_design_30.gpkg")

d.predictors.val <- extract(r.stack, v.validation)
d.SOCstock.val <- extract(r.soc, v.validation)

d.SOCstock.val$predictions <- predict(m.rf, d.predictors.val)$prediction
summary(d.SOCstock.val$predictions)


## -----------------------------------------------------------------------------------------------------------------------
# compute marginal bias (positive values: overestimation, negative values: underestimation)
error <- d.SOCstock.val$`soilgrids_SOCstock_0-30_mean30mSmooth` - d.SOCstock.val$predictions
mean(-error)

# RMSE
sqrt(mean(error^2))

# R2, computed as skill score / model efficiency coefficient 
obs <- d.SOCstock.val$`soilgrids_SOCstock_0-30_mean30mSmooth`
1 - ( sum( error^2 ) / sum( (obs - mean(obs))^2 ) )


## ----fig.width=4, fig.height=4------------------------------------------------------------------------------------------
library(ggplot2)
ggplot(d.SOCstock.val, aes(predictions, `soilgrids_SOCstock_0-30_mean30mSmooth`)) +
  geom_point(size = 2) + ylim(52,108) + xlim(52,108) +
  xlab(expression("predicted SOC stocks [t" ~ha^-1~"]")) +
  ylab(expression("observed SOC stocks [t" ~ha^-1~"]")) +
  stat_smooth(se = FALSE, span = 1.5) + 
  geom_abline (slope=1, linetype = "dashed", color="black") +
  theme_minimal()


## -----------------------------------------------------------------------------------------------------------------------
# create a data frame, each row corresponds to a pixel
d.pixels <- as.data.frame(r.stack, xy = T, na.rm = F)

# remove NA
d.pixels <- d.pixels[ complete.cases(d.pixels), ]

# check if we have all predictors, list missing columns 
names(importance(m.rf))[ !names(importance(m.rf)) %in% names(d.pixels) ]

# predict using ranger object 
d.pixels$SOCstock_predictions <- predict(m.rf, d.pixels)$predictions

# create a raster object and save 
r.pred <- rast(d.pixels[, c("x", "y", "SOCstock_predictions")], type="xyz", crs= "epsg:32740")
writeRaster(r.pred, filename = "datasets/lareunion/results/predicted_SOC_stocks.tif", overwrite = T)


## -----------------------------------------------------------------------------------------------------------------------
plot(r.pred)