{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "kernelspec": {
   "name": "ir",
   "display_name": "R",
   "language": "R"
  },
  "language_info": {
   "name": "R"
  },
  "colab": {
   "provenance": [],
   "toc_visible": true
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3 · Digital soil mapping with spectral predictions\n",
    "**Module 4 · Soil spectroscopy for digital soil mapping** · SoilFER Training\n",
    "\n",
    "[Course page](https://training.yigini.net/modules/04-soil-spectroscopy/) · [Manual chapter](https://training.yigini.net/manual/soil-spectroscopy-for-digital-soil-mapping.html) · [Original script](https://github.com/SoilFER/SoilFER-Training-Resources/blob/main/02_scripts/module4/003-Soil%20Spectroscopy%20for%20Digital%20Soil%20Mapping.R)\n",
    "\n",
    "**Before you start**\n",
    "1. Check that the runtime is **R**: *Runtime → Change runtime type → R* (this notebook should open in R automatically).\n",
    "2. Run the **Setup** cell below once per session (≈1–3 minutes). It downloads the training project, installs the R packages and downloads the course rasters and MIR data (≈1.2 GB) and sets the working folder.\n",
    "3. Then run the cells in order with **Shift + Enter**.\n",
    "\n",
    "> Colab resets when you close it or after ~90 minutes without activity. Save results you want to keep with *Files → Download* (left sidebar), or re-run the Setup cell after a reset.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ==== SoilFER · Colab setup (run first, once per session) ==================\n",
    "options(repos = c(CRAN = \"https://cloud.r-project.org\"), timeout = 3600)\n",
    "# Works in Google Colab and in any other Jupyter (JupyterHub, JupyterLab on your computer)\n",
    "on_colab <- nzchar(Sys.getenv(\"COLAB_RELEASE_TAG\")) || dir.exists(\"/content/sample_data\")\n",
    "root <- if (on_colab) \"/content/SoilFER-Training-Resources\" else path.expand(\"~/SoilFER-Training-Resources\")\n",
    "Sys.setenv(SOILFER_ROOT = root)\n",
    "\n",
    "# 1. Training project (scripts, small data, outputs, assignments)\n",
    "if (!dir.exists(root))\n",
    "  system(paste(\"git clone --depth 1 https://github.com/SoilFER/SoilFER-Training-Resources\", root))\n",
    "\n",
    "# 2. R packages (in Colab the R runtime installs ready-made binaries, so this is fast)\n",
    "pkgs <- c(\"aqp\", \"dplyr\", \"matrixStats\", \"prospectr\", \"ranger\", \"readxl\", \"terra\")\n",
    "need <- setdiff(pkgs, rownames(installed.packages()))\n",
    "if (length(need)) install.packages(need)\n",
    "\n",
    "# 3. Course rasters + MIR spectra (Google Drive folder of the SoilFER training, ≈1.2 GB)\n",
    "td <- file.path(root, \"01_data/module1/training_data\")\n",
    "if (!file.exists(file.path(td, \"MIR_KANSAS_data.xlsx\"))) {\n",
    "  system(\"python3 -m pip -q install gdown\")\n",
    "  drv <- file.path(dirname(root), \"soilfer_drive\")\n",
    "  system(paste(\"python3 -m gdown --folder --quiet https://drive.google.com/drive/folders/1K7tq9zX5HsqbqWcNoT27WtfPtehcKBCu -O\", shQuote(drv)))\n",
    "  f <- list.files(drv, recursive = TRUE, full.names = TRUE)\n",
    "  file.copy(f, td, overwrite = FALSE)\n",
    "}\n",
    "\n",
    "setwd(root)\n",
    "cat(\"Ready. Working folder:\", getwd(), \"\\n\")\n",
    "missing <- setdiff(pkgs, rownames(installed.packages()))\n",
    "if (length(missing)) message(\"Not installed: \", paste(missing, collapse = \", \"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fresh runtime? Re-create the objects from the previous session(s)\n",
    "This session uses objects created earlier. If you just opened this notebook (or Colab was reset), run this cell first. It takes a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "for (s in c(\"02_scripts/module4/001-Soil Spectroscopy for Digital Soil Mapping.R\", \"02_scripts/module4/002-Soil Spectroscopy for Digital Soil Mapping.R\")) { message(\"Running \", s); source(file.path(root, s), echo = FALSE) }\n",
    "setwd(root)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CODE FOR SOILFER MODULE 4:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## CODE FOR SOILFER MODULE 4: \n",
    "## SOIL SECTROSCOPY FOR DIGITAL SOIL MAPPING\n",
    "## Alex Wadoux and Leonardo Ramirez-Lopez\n",
    "\n",
    "## Digital Soil Mapping part\n",
    "\n",
    "## Preceeding code: 002-Soil Spectroscopy for Digital Soil Mapping.R\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load covariates and extract"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Load covariates and extract ###################################\n",
    "################################################################################\n",
    "\n",
    "library(terra)\n",
    "\n",
    "covs <- rast(\n",
    "  \"../SoilFER-Training-Resources/01_data/module1/training_data/Environmental_Covariates_250m_KANSAS.tif\"\n",
    ")\n",
    "cov_names <- names(covs)\n",
    "\n",
    "# Reference (lab only)\n",
    "soil_pts_ref <- vect(\n",
    "  soc_ref_0_30_xy,\n",
    "  geom = c(\"Long_Site.x\", \"Lat_Site.x\"),\n",
    "  crs = \"EPSG:4326\"\n",
    ")\n",
    "soil_pts_ref <- project(soil_pts_ref, covs)\n",
    "soil_pts_ref <- cbind(soil_pts_ref, terra::extract(covs, soil_pts_ref)[, -1])\n",
    "\n",
    "# Augmented (lab + DRS)\n",
    "soil_pts_aug <- vect(\n",
    "  soc_aug_0_30_xy,\n",
    "  geom = c(\"Long_Site.x\", \"Lat_Site.x\"),\n",
    "  crs = \"EPSG:4326\"\n",
    ")\n",
    "soil_pts_aug <- project(soil_pts_aug, covs)\n",
    "soil_pts_aug <- cbind(soil_pts_aug, terra::extract(covs, soil_pts_aug)[, -1])\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute weights"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Compute weights ###############################################\n",
    "################################################################################\n",
    "\n",
    "soil_pts_ref$weight <- 1 / soil_pts_ref$SOC_sd_0_30^2\n",
    "soil_pts_aug$weight <- 1 / soil_pts_aug$SOC_sd_0_30^2\n",
    "\n",
    "# Normalise by the maximum weight\n",
    "soil_pts_ref$weight <- soil_pts_ref$weight / max(soil_pts_ref$weight, na.rm = T)\n",
    "soil_pts_aug$weight <- soil_pts_aug$weight / max(soil_pts_aug$weight, na.rm = T)\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baseline model (lab only)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Baseline model (lab only) #####################################\n",
    "################################################################################\n",
    "\n",
    "library(ranger)\n",
    "\n",
    "soil_df_ref <- as.data.frame(soil_pts_ref, geom = \"XY\")[, c(\"x\", \"y\", \"SOC_0_30\", \"SOC_sd_0_30\", \"weight\", cov_names)\n",
    "]\n",
    "soil_df_ref <- soil_df_ref[!is.na(soil_df_ref$SOC_0_30), ]\n",
    "\n",
    "form <- as.formula(paste(\"SOC_0_30 ~\", paste(cov_names, collapse = \" + \")))\n",
    "\n",
    "mod_ref <- ranger(\n",
    "  formula = form, data = soil_df_ref,\n",
    "  case.weights = soil_df_ref$weight,\n",
    "  replace = FALSE, sample.fraction = 0.632,\n",
    "  num.trees = 500, seed = 201909\n",
    ")\n",
    "\n",
    "soc_map_ref <- predict(covs, mod_ref, na.rm = TRUE)\n",
    "terra::writeRaster(soc_map_ref, \"soc_map_ref.tif\", overwrite = TRUE)\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Augmented model (lab + DRS)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Augmented model (lab + DRS) ###################################\n",
    "################################################################################\n",
    "\n",
    "soil_df_aug <- as.data.frame(soil_pts_aug, geom = \"XY\")[, c(\"x\", \"y\", \"SOC_0_30\", \"SOC_sd_0_30\", \"weight\", cov_names)\n",
    "]\n",
    "soil_df_aug <- soil_df_aug[!is.na(soil_df_aug$SOC_0_30), ]\n",
    "\n",
    "mod_aug <- ranger(\n",
    "  formula = form, data = soil_df_aug,\n",
    "  case.weights = soil_df_aug$weight,\n",
    "  replace = FALSE, sample.fraction = 0.632,\n",
    "  num.trees = 500, seed = 201909\n",
    ")\n",
    "\n",
    "soc_map_aug <- predict(covs, mod_aug, na.rm = TRUE)\n",
    "terra::writeRaster(soc_map_aug, \"soc_map_aug.tif\", overwrite = TRUE)\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weighted performance helper functions"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Weighted performance helper functions #########################\n",
    "################################################################################\n",
    "\n",
    "wME <- function(obs, pred, w) {\n",
    "  sum(w * (pred - obs), na.rm = TRUE) / sum(w, na.rm = TRUE)\n",
    "}\n",
    "\n",
    "wRMSE <- function(obs, pred, w) {\n",
    "  sqrt(sum(w * (pred - obs)^2, na.rm = TRUE) / sum(w, na.rm = TRUE))\n",
    "}\n",
    "\n",
    "wNSE <- function(obs, pred, w) {\n",
    "  zbar_w <- sum(w * obs, na.rm = TRUE) / sum(w, na.rm = TRUE)\n",
    "  1 - sum(w * (pred - obs)^2, na.rm = TRUE) /\n",
    "    sum(w * (obs - zbar_w)^2, na.rm = TRUE)\n",
    "}\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10-fold cross-validation of DSM models"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ 10-fold cross-validation of DSM models ########################\n",
    "################################################################################\n",
    "\n",
    "run_cv <- function(data, response, predictors, nfolds = 10, seed = 201909) {\n",
    "  set.seed(seed)\n",
    "  n <- nrow(data)\n",
    "  fold_id <- sample(rep(seq_len(nfolds), length.out = n))\n",
    "  cv_pred <- rep(NA_real_, n)\n",
    "  form <- as.formula(paste(response, \"~\", paste(predictors, collapse = \" + \")))\n",
    "  \n",
    "  for (fold in seq_len(nfolds)) {\n",
    "    train <- data[fold_id != fold, ]\n",
    "    test <- data[fold_id == fold, ]\n",
    "    mod <- ranger(\n",
    "      formula = form, data = train,\n",
    "      case.weights = train$weight,\n",
    "      replace = FALSE, sample.fraction = 0.632,\n",
    "      seed = seed\n",
    "    )\n",
    "    cv_pred[fold_id == fold] <- predict(mod, data = test)$predictions\n",
    "  }\n",
    "  \n",
    "  data.frame(obs = data[[response]], pred = cv_pred, weight = data$weight)\n",
    "}\n",
    "\n",
    "cv_ref <- run_cv(soil_df_ref, \"SOC_0_30\", cov_names)\n",
    "cv_aug <- run_cv(soil_df_aug, \"SOC_0_30\", cov_names)\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare the two DSM models"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Compare the two DSM models ####################################\n",
    "################################################################################\n",
    "\n",
    "compute_metrics <- function(cv) {\n",
    "  data.frame(\n",
    "    wME = round(wME(cv$obs, cv$pred, cv$weight), 3),\n",
    "    wRMSE = round(wRMSE(cv$obs, cv$pred, cv$weight), 3),\n",
    "    wNSE = round(wNSE(cv$obs, cv$pred, cv$weight), 3)\n",
    "  )\n",
    "}\n",
    "\n",
    "metrics <- rbind(\n",
    "  data.frame(Model = \"Reference (lab only)\", compute_metrics(cv_ref)),\n",
    "  data.frame(Model = \"Augmented (lab + DRS)\", compute_metrics(cv_aug))\n",
    ")\n",
    "\n",
    "metrics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore your results\n",
    "This cell finds the maps and layers you created in this session, plots them here, and packs them into **`soilfer_outputs.zip`** (rasters as GeoTIFF, vector layers as GeoJSON).\n",
    "\n",
    "1. Download `soilfer_outputs.zip` from the **Files** panel on the left (⋮ → *Download*).\n",
    "2. Open the **[SoilFER Map viewer](https://training.yigini.net/viewer.html)** and drop the zip on it: change colours and value ranges, switch bands, overlay layers and click to read values. Nothing is uploaded.\n",
    "3. The same files also open in QGIS on your own computer.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ==== Explore your results ====================================================\n",
    "root <- Sys.getenv(\"SOILFER_ROOT\", \"/content/SoilFER-Training-Resources\")\n",
    "library(terra); library(sf)\n",
    "\n",
    "# New or changed files in the training project = what you produced in this session\n",
    "st  <- system(paste(\"git -C\", shQuote(root), \"-c core.quotePath=off ls-files --others --modified --exclude-standard\"), intern = TRUE)\n",
    "new <- unique(file.path(root, st))\n",
    "new <- new[file.exists(new) & grepl(\"\\\\.(tif|tiff|gpkg|shp|geojson)$\", new, ignore.case = TRUE)]\n",
    "rasters <- new[grepl(\"\\\\.tiff?$\", new, ignore.case = TRUE)]\n",
    "vectors <- new[!grepl(\"\\\\.tiff?$\", new, ignore.case = TRUE)]\n",
    "if (!length(new)) stop(\"No new maps found yet: run the notebook cells above first.\")\n",
    "print(data.frame(file = sub(paste0(root, \"/\"), \"\", new), MB = round(file.size(new) / 1e6, 1)))\n",
    "\n",
    "# 1. Quick look (first band of up to 4 rasters)\n",
    "for (f in head(rasters, 4)) plot(rast(f)[[1]], main = basename(f))\n",
    "\n",
    "# 2. Pack everything for the SoilFER Map viewer / QGIS\n",
    "dir.create(pack <- file.path(dirname(root), \"soilfer_outputs\"), showWarnings = FALSE)\n",
    "unlink(list.files(pack, full.names = TRUE))\n",
    "file.copy(rasters, pack, overwrite = TRUE)\n",
    "for (v in vectors) try(st_write(st_transform(st_read(v, quiet = TRUE), 4326),\n",
    "                                file.path(pack, sub(\"\\\\.[^.]+$\", \".geojson\", basename(v))),\n",
    "                                delete_dsn = TRUE, quiet = TRUE), silent = TRUE)\n",
    "zipf <- file.path(dirname(root), \"soilfer_outputs.zip\"); unlink(zipf)\n",
    "zip(zipf, list.files(pack, full.names = TRUE), flags = \"-j9Xq\")\n",
    "cat(sprintf(\"\\nsoilfer_outputs.zip: %d files, %.1f MB\\nDownload it from the Files panel and drop it on https://training.yigini.net/viewer.html\\n\",\n",
    "            length(list.files(pack)), file.size(zipf) / 1e6))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive map\n",
    "Zoom and pan on the first raster (change `rasters[1]` to `rasters[2]`, … to see the others)."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "library(mapview)\n",
    "mapviewOptions(fgb = FALSE, georaster = FALSE)\n",
    "mapview(rast(rasters[1])[[1]], layer.name = basename(rasters[1]), maxpixels = 4e5)\n"
   ]
  }
 ]
}