{ "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 2 · Spectroscopy calibration models\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/002-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\")\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\")) { 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", "## Soil spectroscopy models part\n", "\n", "## Preceeding code: 001-Soil Spectroscopy for Digital Soil Mapping.R\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Select samples and split" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Select samples and split ######################################\n", "################################################################################\n", "\n", "n_cal <- 85\n", "\n", "set.seed(201909)\n", "kms <- prospectr::naes(\n", " dat$spc_processed,\n", " k = n_cal,\n", " pc = 10,\n", " iter.max = 1000\n", ")\n", "\n", "# Profile-complete the selection\n", "kms_samples <- which(dat$ProfID %in% unique(dat$ProfID[kms$model]))\n", "\n", "# Split into calibration / prediction subsets\n", "dat_cal <- dat[kms_samples, ]\n", "dat_pred <- dat[-kms_samples, ]\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define the CV folds" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Define the CV folds ###########################################\n", "################################################################################\n", "\n", "set.seed(14092019)\n", "nfolds <- 10\n", "\n", "# Assign folds at the profile level (not at the observation level)\n", "cal_profile_ids <- unique(dat_cal$ProfID)\n", "fold_index <- sample(rep(1:nfolds, length.out = length(cal_profile_ids)))\n", "\n", "# Each observation inherits the fold of its parent profile\n", "dat_cal$fold <- fold_index[match(dat_cal$ProfID, cal_profile_ids)]\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run 10-fold cross-validation" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Run 10-fold cross-validation ##################################\n", "################################################################################\n", "\n", "library(ranger)\n", "library(matrixStats)\n", "\n", "SOC_predRF <- rep(NA, nrow(dat_cal))\n", "SOC_varRF <- rep(NA, nrow(dat_cal))\n", "SOC_sdRF <- rep(NA, nrow(dat_cal))\n", "\n", "for (k in seq_len(nfolds)) {\n", " test_rows <- which(dat_cal$fold == k)\n", " train_rows <- which(dat_cal$fold != k)\n", " \n", " model <- ranger(\n", " x = dat_cal$spc_processed[train_rows, ],\n", " y = dat_cal$SOC[train_rows],\n", " quantreg = TRUE, num.trees = 3500,\n", " sample.fraction = 1, replace = TRUE,\n", " splitrule = \"maxstat\", min.node.size = 10,\n", " seed = 201909\n", " )\n", " \n", " pred_matrix <- predict(\n", " model, data = dat_cal$spc_processed[test_rows, ],\n", " type = \"quantiles\",\n", " what = function(x) sample(x, 100, replace = TRUE)\n", " )$predictions\n", " \n", " SOC_predRF[test_rows] <- rowMeans(pred_matrix)\n", " SOC_varRF[test_rows] <- rowVars(pred_matrix)\n", " SOC_sdRF[test_rows] <- sqrt(SOC_varRF[test_rows])\n", "}\n", "\n", "dat_cal$SOC_predRF <- SOC_predRF\n", "dat_cal$SOC_varRF <- SOC_varRF\n", "dat_cal$SOC_sdRF <- SOC_sdRF\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The validation plot (base R version)" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ The validation plot (base R version) ##########################\n", "################################################################################\n", "\n", "xy_lims <- range(dat_cal$SOC_predRF, dat_cal$SOC)\n", "\n", "plot(\n", " x = dat_cal$SOC_predRF,\n", " y = dat_cal$SOC,\n", " xlab = \"Predicted SOC (%)\",\n", " ylab = \"Measured SOC (%)\",\n", " xlim = xy_lims,\n", " ylim = xy_lims,\n", " col = rgb(0, 0, 0, 0.5),\n", " pch = 16,\n", " cex = 1.5\n", ")\n", "\n", "# Error bars (prediction SD)\n", "arrows(\n", " x0 = dat_cal$SOC_predRF - dat_cal$SOC_sdRF,\n", " y0 = dat_cal$SOC,\n", " x1 = dat_cal$SOC_predRF + dat_cal$SOC_sdRF,\n", " y1 = dat_cal$SOC,\n", " code = 3,\n", " length = 0,\n", " col = rgb(0, 0, 0, 0.5)\n", ")\n", "\n", "abline(0, 1, col = col_ref, lwd = 1.5, lty = 2)\n", "grid(col = rgb(0.8, 0.8, 0.8, 0.6), lty = 1)\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computing the error (RMSE)" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Computing the error (RMSE) ####################################\n", "################################################################################\n", "\n", "drs_rmse <- sqrt(mean((dat_cal$SOC_predRF - dat_cal$SOC)^2, na.rm = TRUE))\n", "cat(\"DRS model RMSE:\", round(drs_rmse, 3), \"% SOC\\n\")\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Final model fitting" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Final model fitting ###########################################\n", "################################################################################\n", "\n", "final_drs_model <- ranger(\n", " x = dat_cal$spc_processed,\n", " y = dat_cal$SOC,\n", " quantreg = TRUE,\n", " num.trees = 3500,\n", " sample.fraction = 1,\n", " replace = TRUE,\n", " splitrule = \"maxstat\",\n", " min.node.size = 10,\n", " seed = 201909\n", ")\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Predict SOC for the remaining 80%" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Predict SOC for the remaining 80% #############################\n", "################################################################################\n", "\n", "drs_soc_preds <- predict(\n", " final_drs_model,\n", " data = dat_pred$spc_processed,\n", " type = \"quantiles\",\n", " what = function(x) sample(x, 100, replace = TRUE)\n", ")$predictions\n", "\n", "dat_pred$SOC_predRF <- rowMeans(drs_soc_preds)\n", "dat_pred$SOC_varRF <- rowVars(drs_soc_preds)\n", "dat_pred$SOC_sdRF <- sqrt(dat_pred$SOC_varRF)\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Assemble the augmented dataset" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Assemble the augmented dataset ################################\n", "################################################################################\n", "\n", "dat_augmented <- dat[, c(\n", " \"smp_id\", \"ProfID\", \"Long_Site.x\", \"Lat_Site.x\",\n", " \"Top_depth_cm.x\", \"Bottom_depth_cm.x\"\n", ")]\n", "\n", "dat_augmented$SOC <- NA\n", "dat_augmented$SOC_sd <- NA\n", "\n", "# Lab SOC for calibration samples\n", "dat_augmented$SOC[kms_samples] <- dat$SOC[kms_samples]\n", "\n", "# DRS-predicted SOC for the rest\n", "dat_augmented$SOC[-kms_samples] <- dat_pred$SOC_predRF\n", "\n", "# DRS prediction SD as uncertainty for predicted samples\n", "dat_augmented$SOC_sd[-kms_samples] <- dat_pred$SOC_sdRF\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Apply lab uncertainty" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Apply lab uncertainty #########################################\n", "################################################################################\n", "\n", "# Lab uncertainty (Stevens et al. 2013)\n", "dat$SOC_sd <- 0.15\n", "\n", "# Same value for calibration samples in the augmented dataset\n", "dat_augmented$SOC_sd[kms_samples] <- 0.15\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Build SoilProfileCollection objects" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Build SoilProfileCollection objects ###########################\n", "################################################################################\n", "\n", "library(aqp)\n", "depths(dat) <- ProfID ~ Top_depth_cm.x + Bottom_depth_cm.x\n", "depths(dat_augmented) <- ProfID ~ Top_depth_cm.x + Bottom_depth_cm.x\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Helper for weighted SD over 0-30 cm" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Helper for weighted SD over 0-30 cm ###########################\n", "################################################################################\n", "\n", "weighted_sd_0_30 <- function(top, bottom, sd, z1 = 0, z2 = 30) {\n", " w <- pmax(0, pmin(bottom, z2) - pmax(top, z1))\n", " keep <- w > 0 & !is.na(sd)\n", " w <- w[keep]\n", " sd <- sd[keep]\n", " if (sum(w) < (z2 - z1)) return(NA)\n", " sqrt(sum((w * sd)^2)) / sum(w)\n", "}\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Harmonise the reference dataset" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Harmonise the reference dataset ###############################\n", "################################################################################\n", "\n", "soc_mean_ref <- profileApply(dat, function(p) {\n", " h <- horizons(p)\n", " w <- pmax(0, pmin(h$Bottom_depth_cm.x, 30) - pmax(h$Top_depth_cm.x, 0))\n", " keep <- w > 0 & !is.na(h$SOC)\n", " if (sum(w[keep]) < 30) return(NA)\n", " sum(w[keep] * h$SOC[keep]) / sum(w[keep])\n", "})\n", "\n", "soc_sd_ref <- horizons(dat) |>\n", " as.data.frame() |>\n", " group_by(ProfID) |>\n", " summarise(\n", " SOC_sd_0_30 = weighted_sd_0_30(Top_depth_cm.x, Bottom_depth_cm.x, SOC_sd),\n", " .groups = \"drop\"\n", " )\n", "\n", "coords_ref <- as.data.frame(dat)[, c(\"ProfID\", \"Lat_Site.x\", \"Long_Site.x\")]\n", "coords_ref <- coords_ref[!duplicated(coords_ref$ProfID), ]\n", "\n", "soc_ref_0_30_xy <- soc_sd_ref |>\n", " mutate(SOC_0_30 = soc_mean_ref) |>\n", " select(ProfID, SOC_0_30, SOC_sd_0_30) |>\n", " merge(coords_ref, by = \"ProfID\", all.x = TRUE)\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Harmonise the augmented dataset" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Harmonise the augmented dataset ###############################\n", "################################################################################\n", "\n", "soc_mean_aug <- profileApply(dat_augmented, function(p) {\n", " h <- horizons(p)\n", " w <- pmax(0, pmin(h$Bottom_depth_cm.x, 30) - pmax(h$Top_depth_cm.x, 0))\n", " keep <- w > 0 & !is.na(h$SOC)\n", " if (sum(w[keep]) < 30) return(NA)\n", " sum(w[keep] * h$SOC[keep]) / sum(w[keep])\n", "})\n", "\n", "soc_sd_aug <- horizons(dat_augmented) |>\n", " as.data.frame() |>\n", " group_by(ProfID) |>\n", " summarise(\n", " SOC_sd_0_30 = weighted_sd_0_30(Top_depth_cm.x, Bottom_depth_cm.x, SOC_sd),\n", " .groups = \"drop\"\n", " )\n", "\n", "coords_aug <- as.data.frame(dat_augmented)[, c(\"ProfID\", \"Lat_Site.x\", \"Long_Site.x\")]\n", "coords_aug <- coords_aug[!duplicated(coords_aug$ProfID), ]\n", "\n", "soc_aug_0_30_xy <- soc_sd_aug |>\n", " mutate(SOC_0_30 = soc_mean_aug) |>\n", " select(ProfID, SOC_0_30, SOC_sd_0_30) |>\n", " merge(coords_aug, by = \"ProfID\", all.x = TRUE)\n", "\n", "\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Caveat on uncertainty propagation" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################ Caveat on uncertainty propagation #############################\n", "################################################################################\n", "\n", "# Lab uncertainty fixed at the profile level (Stevens et al. 2013)\n", "soc_ref_0_30_xy$SOC_sd_0_30 <- 0.15\n", "\n", "soc_aug_0_30_xy$SOC_sd_0_30[\n", " soc_aug_0_30_xy$ProfID %in% unique(dat$ProfID[kms_samples])\n", "] <- 0.15\n" ] } ] }