{
 "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 1 · Read and pre-process spectra\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/001-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(\"dplyr\", \"prospectr\", \"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": [
    "### 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",
    "## Read and pre-process spectra part\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load and plot the DRS data"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Load and plot the DRS data ####################################\n",
    "################################################################################\n",
    "\n",
    "library(readxl)\n",
    "library(dplyr)\n",
    "\n",
    "original_dat <- read_excel(\n",
    "  \"../SoilFER-Training-Resources/01_data/module1/training_data/MIR_KANSAS_data.xlsx\",\n",
    "  .name_repair = \"minimal\",\n",
    "  na = c(\"\", \"NA\", \"N/A\", \"NaN\"),\n",
    "  progress = FALSE\n",
    ")\n",
    "\n",
    "original_dat <- as.data.frame(original_dat)\n",
    "\n",
    "dim(original_dat)\n",
    "sum(duplicated(original_dat$smp_id))\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Separate spectral from non-spectral data"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Separate spectral from non-spectral data ######################\n",
    "################################################################################\n",
    "\n",
    "# Identify spectral columns\n",
    "spc_cols <- grep(\"^X[0-9]{3,}\", colnames(original_dat))\n",
    "\n",
    "# Extract spectral matrix\n",
    "spc_original <- original_dat[, spc_cols]\n",
    "colnames(spc_original) <- gsub(\"^X\", \"\", colnames(spc_original))\n",
    "original_wavs <- as.numeric(colnames(spc_original))\n",
    "\n",
    "range(original_wavs)\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Resample spectra"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Resample spectra ##############################################\n",
    "################################################################################\n",
    "\n",
    "library(prospectr)\n",
    "\n",
    "wavs <- seq(600, 3992, by = 8)\n",
    "\n",
    "spc_resampled <- prospectr::resample(\n",
    "  spc_original,\n",
    "  wav = original_wavs,\n",
    "  new.wav = wavs\n",
    ")\n",
    "\n",
    "dim(spc_resampled)\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Attach spectra and average replicates"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Attach spectra and average replicates #########################\n",
    "################################################################################\n",
    "\n",
    "repeated_dat <- original_dat[, -spc_cols]\n",
    "repeated_dat$spc <- spc_resampled\n",
    "\n",
    "# Average replicate spectra by sample ID\n",
    "spc_avg <- aggregate(\n",
    "  repeated_dat$spc,\n",
    "  by = list(smp_id = repeated_dat$smp_id),\n",
    "  FUN = mean\n",
    ")\n",
    "spc_avg <- spc_avg[order(spc_avg$smp_id), ]\n",
    "rownames(spc_avg) <- spc_avg$smp_id\n",
    "spc_avg <- as.matrix(spc_avg[, -1])\n",
    "\n",
    "dat <- repeated_dat[!duplicated(repeated_dat$smp_id), ]\n",
    "dat <- dat[order(dat$smp_id), ]\n",
    "stopifnot(all(rownames(spc_avg) == dat$smp_id))\n",
    "dat$spc <- spc_avg\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the spectra"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Plot the spectra ##############################################\n",
    "################################################################################\n",
    "\n",
    "matplot(\n",
    "  wavs,\n",
    "  t(dat$spc),\n",
    "  type = \"l\",\n",
    "  lty = 1,\n",
    "  xlim = c(4000, 500),\n",
    "  xlab = expression(Wavenumber~(cm^{-1})),\n",
    "  ylab = \"Absorbance\",\n",
    "  col = rgb(0.8, 0, 0, 0.3)\n",
    ")\n",
    "grid(lty = 1)\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Profile identifier"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Profile identifier ############################################\n",
    "################################################################################\n",
    "\n",
    "coords <- paste(dat$Long_Site.x, dat$Lat_Site.x, sep = \"_\")\n",
    "prof_num <- match(coords, unique(coords))\n",
    "dat$ProfID <- sprintf(\"PROF%04d\", prof_num)\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: one spectrum, raw"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Example: one spectrum, raw ####################################\n",
    "################################################################################\n",
    "\n",
    "a_spectrum <- dat$spc[1, ]\n",
    "plot(\n",
    "  wavs, a_spectrum,\n",
    "  type = \"l\", lty = 1,\n",
    "  xlim = c(4000, 500),\n",
    "  xlab = expression(Wavenumber~(cm^{-1})),\n",
    "  ylab = \"Absorbance\",\n",
    "  col = rgb(0.3, 0, 0.8)\n",
    ")\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First derivative (Savitzky-Golay)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ First derivative (Savitzky-Golay) #############################\n",
    "################################################################################\n",
    "\n",
    "a_spectrum_first_der <- prospectr::savitzkyGolay(\n",
    "  a_spectrum, m = 1, p = 1, w = 25\n",
    ")\n",
    "der_wavs <- as.numeric(names(a_spectrum_first_der))\n",
    "\n",
    "plot(\n",
    "  der_wavs, a_spectrum_first_der,\n",
    "  type = \"l\", lty = 1,\n",
    "  xlim = c(4000, 500),\n",
    "  xlab = expression(Wavenumber~(cm^{-1})),\n",
    "  ylab = \"First derivative absorbance\",\n",
    "  col = rgb(0, 0.3, 0.8)\n",
    ")\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### After SNV correction"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ After SNV correction ##########################################\n",
    "################################################################################\n",
    "\n",
    "a_spectrum_first_der_snv <- prospectr::standardNormalVariate(\n",
    "  matrix(a_spectrum_first_der, nrow = 1)\n",
    ")\n",
    "\n",
    "plot(\n",
    "  der_wavs, a_spectrum_first_der_snv,\n",
    "  type = \"l\", lty = 1,\n",
    "  xlim = c(4000, 500),\n",
    "  xlab = expression(Wavenumber~(cm^{-1})),\n",
    "  ylab = \"SNV of first derivative absorbance\",\n",
    "  col = rgb(0.8, 0.3, 0)\n",
    ")\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Apply preprocessing to all spectra"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Apply preprocessing to all spectra ############################\n",
    "################################################################################\n",
    "\n",
    "dat$spc_processed <- dat$spc |>\n",
    "  prospectr::savitzkyGolay(m = 1, p = 1, w = 25) |>\n",
    "  prospectr::standardNormalVariate()\n",
    "\n",
    "wavs_processed <- as.numeric(colnames(dat$spc_processed))\n",
    "\n",
    "matplot(\n",
    "  wavs_processed, t(dat$spc_processed),\n",
    "  type = \"l\", lty = 1,\n",
    "  xlim = c(4000, 500),\n",
    "  xlab = expression(Wavenumber~(cm^{-1})),\n",
    "  ylab = \"SNV of first derivative absorbance\",\n",
    "  col = rgb(0.3, 0, 0.8, 0.3)\n",
    ")\n",
    "grid(lty = 1)\n"
   ]
  }
 ]
}