02_semi-spatial-gam-climate.Rmd

---
editor_options: 
  chunk_output_type: console
---

# Building and validating semi-spatial GAMs for climate

The overall aim here is to model historical climate variables --- mean temperature, standard deviation in temperature, and total precipitation --- as functions of geographical characteristics of the locations with which the data are associated.

We fit three candidate models of increasing complexity.
Then, the models are used to predict historical climate over the entire study area; we refer to these as semi-spatial climate models, as they include additional predictors over the coordinates of the locations alone.

The climate data extend into the modern period, allowing predictions from the candidate models to be compared against data from remote sensing, in order to select the best model formulation for each climate variable in each season. 

## Load libraries

```{r}
# Load libraries
library(tidyverse)
library(sf)
library(readxl)
library(mgcv) # for GAMs
library(terra) # for spatial rasters
```

## Prepare linked climate-spatial data and spatial data for modelling

We load the historical climate data with linked geographical covariates (hereon, climate data), and the spatial data raster stack.

We prepare the covariates for modelling: distance from the coast is converted from metres to kilometres, and a small amount of random noise is added to all values. This is necessary to allow for fitting as otherwise all predictor values are identical.

```{r}
# read data
data <- read_csv("results/climate/data_for_gam.csv")
# `stack` is the raster stack of geographical predictors, and not the climate
# model predictions.
stack <- terra::rast("results/climate/raster_gam_stack.tif")
```

```{r}
# prepare data for gam fitting
data <- mutate(
  data,

  # divide distance by 1000 for km
  coast = coast / 1000,

  # add some error to all physical variables
  # hopefully prevents model issues
  coast = coast + rnorm(length(coast), 0.1, sd = 0.01),
  elev = elev + rnorm(length(coast), 1, sd = 0.01),
  lat = lat + rnorm(length(lat), 0, sd = 0.01)
)

# set distance to coast raster in km
values(stack[["coast"]]) <- values(stack[["coast"]]) / 1000
```

We pivot the data to long format, and then nest the data by season and decade, gving us smaller datatsets for which models will be fit separately.

While there is no specific reason to believe that the relationship between climate and geographical predictors has changed over the past 150 years, splitting the data by decade may help adjust for long-term climate system effects within certain time-periods that could be drowned out by the overall effect of geographical predictors.

We subset the data for the decades between 1970 and 2010 for which reliable comparator spatial climate datasets are available from remote sensing.

```{r}
# Pivot to long format, remove unnecessary variables,
data <- pivot_longer(
  data,
  cols = c("ppt", "t_mean", "t_sd"),
  names_to = "climvar"
)
data <- dplyr::select(data, !month)

# Select period 1970 -- 2010 where we also have satellite data
data <- filter(data, between(year, 1970, 2010))
```

```{r}
data <- nest(
  data,
  data = !c(climvar, season, year_bin)
)
```

## Fit candidate spatial-climate models

We prepare three candidate spatial-climate models for each climate variable:

- A model where the variable is only a smooth function of elevation;

- A model where the variable is a smooth function of elevation, and a linear function of distance from the coast and latitude; and

- A model where the variable is a smooth function of elevation, and of the interaction between latitude and distance from the coast.

We create combinations of each climate variable dataset (further split by decade and season) and each model formulation, and fit the variable to candidate predictors using generalised additive models (GAMs) from [the _mgcv_ package](https://cran.r-project.org/package=mgcv).

```{r}
# elevation with 3 knots
form_elev_model <- value ~ s(elev, k = 3)

# elevation and coast
form_elev_coast <- value ~ s(elev, k = 3) + coast + lat

# elevation and coast with lat
form_elev_colat <- value ~ s(elev, k = 3) + s(coast, lat, k = 5)
```

```{r}
# make combinations
data <- crossing(
  data,
  forms = c(form_elev_model, form_elev_coast, form_elev_colat)
)
```

```{r}
# fit models
data <- mutate(
  data,
  mod = map2(data, forms, function(df, form) {
    gam(
      formula = form,
      data = df
    )
  })
)
```

```{r}
# Save the fitted model data as an Rdata object
save(data, file = "results/climate/model_gam_climate.Rds")
```

## Model predictions for modern climate

We use the fitted models to predict climatic variables over the remainder of the study area.

```{r}
# ensure that NAs are removed from the stack object before running the model
# projections
terra::global(stack, fun = "isNA")

# this tells us that two raster layers within the stack object have NA values

stack$coast[is.na(stack$coast)] <- 0
stack$elev[stack$elev < 0] <- NA
stack$elev[is.na(stack$elev)] <- 0
stack$elev <- as.numeric(stack$elev)

# save model predictions
model_pred <- map(
  data$mod, function(g) {
    predict(stack, g, type = "response") # order matters here
  }
)

# `Reduce()` with function `c()` combines model predictions into a
# single raster stack, which is then saved
model_pred <- Reduce(model_pred, f = c)
terra::writeRaster(model_pred, "results/climate/data_gam_pred.tif")
```

```{r}
# Read in the saved model predictions
# this allows us to start at this stage without needing to run
# previous steps
model_pred <- terra::rast("results/climate/data_gam_pred.tif")
```

## Average predictions over 1970 to 2010

We combine the model predictions with the appropriate identifier variables (decade, season, climate variable, and model formulation) in the nested climate data dataframe, and get the average predicted values for each variable in each season.

```{r}
# assign as list column
data <- mutate(
  data,
  pred = as.list(model_pred)
)

# Summarise rasters using a reduce over the list
# basically gets the average over 1970 -- 2010
data_pred <- ungroup(data) %>%
  group_by(season, climvar, forms) %>%
  summarise(
    mean_pred = list(Reduce(pred, f = mean))
  )

# make raster stack
gam_pred_avg <- data_pred$mean_pred
gam_pred_avg <- Reduce(f = c, gam_pred_avg)

# save averaged predictions
terra::writeRaster(
  gam_pred_avg,
  filename = "results/climate/gam_pred_valid_avg.tif"
)
```

We save the dataframe of season, climate variable, and model formula for future reference.

```{r}
# load saved pred
gam_pred_avg <- terra::rast("results/climate/gam_pred_valid_avg.tif")

# write data for linking
data_pred |>
  # ungroup the data
  ungroup() |>
  # get unique season, variable, and model formula
  distinct(season, climvar, forms) |>
  # then get the model formula as a character
  mutate(
    forms = map_chr(forms, function(l) {
      str_flatten(as.character(l)[c(2, 1, 3)])
    })
  ) |>
  # save to file
  write_csv(
    file = "results/climate/gam_model_formulas.csv"
  )
```

## Get BIOCLIM data

In this section, we acquire remotely sensed cliamte variables for the study area, which will serve as a comparator for the values generated from fitted model predictions, and will help to pick the most appropriate model formulation to reconstruct historical climate for which remotely sensed data are not available.

First, we obtained bio-climatic data from CHELSA. We do not upload the rasters to GitHub as they are very large in size and these can be downloaded from this link: https://envicloud.wsl.ch/#/?prefix=chelsa%2Fchelsa_V1

## Temperature data

We process monthly minimum and maximum CHELSA temperature data to get the mean monthly temperature in the dry and wet season as previously defined.

```{r}
# Define the patterns for data to search local data files
patterns <- c("tmin", "tmax")

# List the filepaths matching the pattern in local data (relative to the proj.)
tkAvg <- map(patterns, function(pattern) {
  # list the paths
  files <- list.files(
    path = "data/climate/chelsa",
    full.names = TRUE,
    recursive = TRUE,
    pattern = pattern
  )
})

# Read found rasters and crop them to the extent of other rasters representing
# the study area; this is the object `stack`
tkAvg <- map(tkAvg, function(paths) {
  # going over the file paths, read them in as rasters, convert CRS and crop
  tempData <- map(paths, function(path) {
    a <- terra::rast(path)
    a <- terra::crop(a, terra::ext(stack))
    a
  })

  # The data are in a native format that needs to be converted before
  # processing.
  # Convert each to Kelvin, first dividing by 10 to get celsius
  tempData <- map(tempData, function(tmpRaster) {
    tmpRaster <- (tmpRaster / 10) + 273.15
    tmpRaster
  })
})

# Iteratively get the mean temperature for each month, and assign
# names to the resulting raster stacks
names(tkAvg) <- patterns

# go over the tmin and tmax and get the average monthly temp
tkAvg <- map2(tkAvg[["tmin"]], tkAvg[["tmax"]], function(tmin, tmax) {
  # return the mean of the corresponding tmin and tmax
  # still in Kelvin
  terra::mean(c(tmin, tmax))
})

# Error if the there are fewer than 12 output layers
assertthat::assert_that(
  length(tkAvg) == 12,
  msg = "temp raster list has fewer than 12 months"
)

# Assign names to identify months
names(tkAvg) <- sprintf("month_%i", seq(12))

# Separate data for the rainy and dry seasons to get the
# rainy and dry season mean monthyl temperature
temp_rainy <- Reduce(tkAvg[seq(6, 11)], f = `c`) |>
  terra::mean()
temp_dry <- Reduce(tkAvg[c(12, seq(5))], f = `c`) |>
  terra::mean()

# Convert values back to celsius
chelsa_t_mean <- c(temp_rainy, temp_dry) - 273.15
names(chelsa_t_mean) <- c("chelsa_temp_rainy_6_11", "chelsa_temp_dry_12_5")

# Save stack
terra::writeRaster(
  chelsa_t_mean,
  filename = "results/climate/chelsa_temp_stack.tif",
  overwrite = TRUE
)
```

## Precipitation data

We process total monthly precipitation data from CHELSA to get the mean monthly temperature in the dry and wet seasons in millimetres.

```{r}
# List precipitation rasters stored locally (relative to this project)
ppt <- list.files(
  path = "data/climate/chelsa",
  full.names = TRUE,
  recursive = TRUE,
  pattern = "prec"
)

# Read each as rasters and crop by extent of study area
ppt <- map(ppt, function(path) {
  a <- terra::rast(path)
  terra::crop(a, terra::ext(stack))
})

# Separate rainy and dry season and get the mean total monthly precipitation
ppt_rainy <- Reduce(ppt[seq(6, 11)], f = `c`) |>
  terra::mean()
ppt_dry <- Reduce(ppt[c(12, seq(5))], f = `c`) |>
  terra::mean()

# Make and save stack
chelsa_ppt_sum <- c(ppt_rainy, ppt_dry)
names(chelsa_ppt_sum) <- c("chelsa_ppt_rainy_6_11", "chelsa_ppt_dry_12_5")

terra::writeRaster(
  chelsa_ppt_sum,
  filename = "results/climate/chelsa_ppt_stack.tif",
  overwrite = TRUE
)
```

## Compare spatial-climate model predictions against BIOCLIM data

## Sample locations

First we draw 10,000 randomly selected coordinates from the extent of the study area, and use these to compare the model predictions against CHELSA-derived data.
We will further stratify these into 10 arbitrary groups at a later stage.

```{r}
# get coordinates from terra
coords <- terra::xyFromCell(
  stack,
  cell = seq(length(values(stack[[1]])))
)

# Sample 1000 locations from the raster data coordinates
# Use `withr::with_local_seed()` to preserve which coordinates are selected
# in each analysis run
coords <- withr::with_seed(
  1,
  coords[sample(1e4, replace = FALSE), ]
) |>
  as_tibble()

# Extract geographical predictors and CHELSA data at locations
# NOTE: recall that `stack` is the raster stack of geographical predictors
sample_locations <-
  mutate(
    coords,
    terra::extract(stack, coords)
  ) |>
  mutate(
    terra::extract(chelsa_ppt_sum, coords)
  ) |>
  mutate(
    terra::extract(chelsa_t_mean, coords)
  )
```

```{r}
# save data
write_csv(
  sample_locations,
  file = "results/climate/data_sample_coords_gam_validation.csv"
)
```

## Link BIOCLIM data and model predictions

Next we link the remote-sensing-derived BIOCLIM data at the sampled locations with the predicted values from our spatial-climate GAMs.

```{r}
# Read in the saved data; this allows continuing from this point without
# re-running the previous code
sample_locations <- read_csv(
  "results/climate/data_sample_coords_gam_validation.csv"
)
```

```{r}
# Pivot the data to long format, retaining coordinates, ID, and geographical
# predictors as identifying variables
sample_locations <-
  sample_locations |>
  pivot_longer(
    cols = !c("x", "y", "ID", "coast", "elev", "lat")
  )

# Assign the season and identify the variable from the variable name
sample_locations <- mutate(
  sample_locations,
  season = str_extract(
    name,
    pattern = "dry|rainy"
  ),
  climvar = str_extract(
    name,
    pattern = "temp|ppt"
  ),
  climvar = if_else(
    climvar == "temp", "t_mean", "ppt"
  )
)
```

We link the spatial-climate models' predictions. Recall that there are three candidate models for each variable, in each season, and the aim is to select a model formulation for the remaining historical data based on minimising deviation from the BIOCLIM data.

```{r}
# Filter the climate-model predictions to remove the SD of temperature
data_pred <- filter(
  data_pred,
  climvar != "t_sd"
)

# Nest model predictions by season and variable
sample_locations <- nest(
  sample_locations,
  chelsa = !c("season", "climvar")
)

# Link predictions for sampled locations with CHELSA data rasters and
# geographical predictors
data_pred <- left_join(
  data_pred,
  sample_locations
)
```

Since the model prediction rasters are stored as-is, we need to extract the values at the coordinates from those rasters.

```{r}
# For each climate variable, extract the predicted values at the sampled
# locations
data_pred <- mutate(
  data_pred,
  chelsa = map2(chelsa, mean_pred, function(ch, pr) {
    ch |>
      rename(
        bioclim_val = "value"
      ) |>
      mutate(
        pred_val = terra::extract(pr, ch[, c("x", "y")])
        # the naming doesn't quite work
      )
  })
)

# Make a copy for further operations
data_gam_validate <- data_pred |>
  select(season, climvar, forms, chelsa)

# Unnest data
data_gam_validate <- unnest(
  data_gam_validate,
  cols = chelsa
)

# Convert model formula from formula to characters
data_gam_validate <- mutate(
  data_gam_validate,
  forms = as.character(forms)
)

# Extract the predicted values
# TODO: check why this is being done again
data_gam_validate <- data_gam_validate |>
  mutate(
    pred_val = pred_val$lyr1
  )

# Drop some variables and save data
write_csv(
  data_gam_validate,
  file = "results/climate/data_gam_validate_compare.csv"
)
```

## Error between model prediction and CHELSA data

We use mean absolute error (MAE) as a measure of how much the model predictions deviate from the CHELSA data. In order to calculate the mean, we group the data into ten arbitrary groups of 1,000 samples each.

We calculate MAE as the mean of the absolute differences between the CHELSA data and the model predicted data, and save the data to a local file.
MAE is calculated for each climate variable, model formulation, and season separately (in addition to the sample group, which is arbitrary).

```{r}
# Group samples into 10 chunks of 1000 coordinates each
data_gam_validate <- group_by(
  data_gam_validate,
  season, climvar, forms
) |>
  mutate(
    group = rep(seq(10), each = 1e3L)
  )

# Grouping by season, climate variable, and model formula, calculate MAE
data_mae <-
  group_by(
    data_gam_validate,
    season, climvar, forms, group
  ) |>
  summarise(
    mae = mean(
      abs(
        bioclim_val - pred_val
      ),
      na.rm = TRUE
    )
  )

# Save data
write_csv(
  data_mae,
  file = "results/climate/data_gam_comparison_mae.csv"
)
```

## GAM predictions at survey sites

In this section, we get the spatial-climate model predictions at the modern survey site locations (called 'survey sites').
This is a repitition of the steps above for the much smaller subset of survery sites.
Recall that for each survey site, there will be three candidate model predictions for each climate variable in each season.

```{r}
# Read survey sites from a local file
survey_sites <- read.csv("data/list-of-resurvey-locations.csv")
names(survey_sites)[4] <- "x"
names(survey_sites)[5] <- "y"

# Sample BIOCLIM rasters at survey sites
survey_sites <- mutate(
  survey_sites,
  terra::extract(
    c(chelsa_t_mean, chelsa_ppt_sum),
    survey_sites[, c("x", "y")]
  )
)

# Remove survery site id, keeping only the 'modern_site_code'
survey_sites <- select(
  survey_sites,
  -c(site_name, historical_site_code, modern_landCover_type, ID)
)

# Pivot the data to long format and identify the season and climate
# variable from the variable name (this is autogenerated by `terra::extract()`)
survey_sites <- pivot_longer(
  survey_sites,
  cols = !c("modern_site_code", "x", "y")
) |>
  mutate(
    season = str_extract(
      name,
      pattern = "dry|rainy"
    ),
    climvar = str_extract(
      name,
      pattern = "temp|ppt"
    ),
    climvar = if_else(
      climvar == "temp", "t_mean", "ppt"
    )
  )

# Nest data, and remove the linked rasters
survey_sites <- nest(
  survey_sites,
  chelsa = !c("season", "climvar")
)
data_pred <- select(data_pred, -chelsa) # remove rasters as data extracted

# Link survey data coordinates with model prediction data
survey_sites <- left_join(
  survey_sites,
  data_pred
)

# For each survey site location, extract model-predicted values from
# the model-prediction rasters
survey_sites <- mutate(
  survey_sites,
  chelsa = map2(chelsa, mean_pred, function(ch, pr) {
    ch |>
      rename(
        bioclim_val = "value"
      ) |>
      mutate(
        pred_val = terra::extract(pr, ch[, c("x", "y")])$lyr1
        # the naming doesn't quite work
      )
  })
)

# Convert the model forumlae in the data to characters and unnest the
# CHELSA data
survey_sites <- survey_sites |>
  select(season, climvar, chelsa, forms) |>
  mutate(
    forms = as.character(forms)
  ) |>
  unnest(
    cols = chelsa
  )

# Save the CHELSA and predicted measures at survey sites
write_csv(
  survey_sites,
  file = "results/climate/data_gam_compare_survey_sites.csv"
)
```