{
 "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": [
    "# Helper · Soil property classes (pH and clay × pH)\n",
    "**Module 3 · Digital soil mapping** · SoilFER Training\n",
    "\n",
    "[Course page](https://training.yigini.net/modules/03-digital-soil-mapping/) · [Manual chapter](https://training.yigini.net/manual/digital-soil-mapping.html) · [Original script](https://github.com/SoilFER/SoilFER-Training-Resources/blob/main/02_scripts/module3/Soil_property_classes_generator.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\", \"ggplot2\")\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": [
    "### Subset soil data and create pH and Clay x pH classes"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Subset soil data and create pH and Clay x pH classes\n",
    "# Author: Generated for SoilFER Training\n",
    "# Date: 2025-01-25\n",
    "\n",
    "library(dplyr)\n",
    "library(ggplot2)      # For plotting\n",
    "\n",
    "# Working directory\n",
    "setwd(file.path(Sys.getenv(\"SOILFER_ROOT\"), \"02_scripts/module3\"))  # folder of this script (was rstudioapi::getActiveDocumentContext, RStudio only)\n",
    "\n",
    "# Read the harmonized data\n",
    "data <- read.csv(\"harmonized.csv\")\n",
    "\n",
    "# Subset: keep only 0-30 cm depth and select relevant columns\n",
    "subset_data <- data %>%\n",
    "  filter(top == 0 & bottom == 30) %>%\n",
    "  select(\n",
    "    ProfID,\n",
    "    lon,\n",
    "    lat,\n",
    "    top,\n",
    "    bottom,\n",
    "    Clay = Clay_p.q50,\n",
    "    Silt = Silt_p.q50,\n",
    "    Sand = Sand_p.q50,\n",
    "    SOC = SOC_p.q50,\n",
    "    pH = pH_p.q50\n",
    "  )\n",
    "\n",
    "# --- Create ordinal pH variable ---\n",
    "subset_data$pH_Class <- cut(\n",
    "  subset_data$pH,\n",
    "  breaks = c(4, 5, 6, 7, 8, 9),\n",
    "  labels = c(\"very_acid\", \"acid\", \"neutral\", \"alkaline\", \"very_alkaline\"),\n",
    "  include.lowest = TRUE,\n",
    "  ordered_result = TRUE\n",
    ")\n",
    "\n",
    "# --- Create nominal variable (pH x Clay segmentation) ---\n",
    "# Clay categories: Low (<20%), Medium (20-40%), High (>40%)\n",
    "# pH categories: Acid (<6), Neutral (6-7), Alkaline (>7)\n",
    "\n",
    "subset_data <- subset_data %>%\n",
    "  mutate(\n",
    "    Clay_Category = case_when(\n",
    "      Clay < 20 ~ \"low_clay\",\n",
    "      Clay >= 20 & Clay < 40 ~ \"medium_clay\",\n",
    "      Clay >= 40 ~ \"high_clay\",\n",
    "      TRUE ~ NA_character_\n",
    "    ),\n",
    "    pH_Category = case_when(\n",
    "      pH < 6 ~ \"acid\",\n",
    "      pH >= 6 & pH < 7 ~ \"neutral\",\n",
    "      pH >= 7 ~ \"alkaline\",\n",
    "      TRUE ~ NA_character_\n",
    "    ),\n",
    "    Clay_pH_Class = case_when(\n",
    "      !is.na(Clay_Category) & !is.na(pH_Category) ~ paste(Clay_Category, pH_Category, sep = \"_\"),\n",
    "      TRUE ~ NA_character_\n",
    "    )\n",
    "  ) %>%\n",
    "  select(-Clay_Category, -pH_Category)\n",
    "\n",
    "# Display summary\n",
    "cat(\"=== Subset Data Summary ===\\n\")\n",
    "cat(\"Number of profiles:\", nrow(subset_data), \"\\n\\n\")\n",
    "\n",
    "cat(\"pH class distribution:\\n\")\n",
    "print(table(subset_data$pH_Class, useNA = \"ifany\"))\n",
    "\n",
    "cat(\"\\nClay x pH class distribution:\\n\")\n",
    "print(table(subset_data$Clay_pH_Class, useNA = \"ifany\"))\n",
    "\n",
    "cat(\"\\n=== First 10 rows ===\\n\")\n",
    "print(head(subset_data, 10))\n",
    "\n",
    "# Save to CSV\n",
    "output_file <- \"harmonized_subset_data.csv\"\n",
    "write.csv(subset_data, output_file, row.names = FALSE)\n",
    "\n",
    "cat(\"\\n\\nFile saved as:\", output_file, \"\\n\")\n",
    "\n",
    "# --- Plot pH Class distribution (bar plot) ---\n",
    "ggplot(subset_data %>% filter(!is.na(pH_Class)), aes(x = pH_Class)) +\n",
    "  geom_bar(fill = \"steelblue\") +\n",
    "  labs(\n",
    "    title = \"pH Class Distribution (0-30 cm)\",\n",
    "    x = \"pH Class\",\n",
    "    y = \"Count\"\n",
    "  ) +\n",
    "  theme_minimal()\n",
    "\n",
    "# --- Plot Clay x pH Class (scatter plot with categories) ---\n",
    "ggplot(subset_data %>% filter(!is.na(Clay_pH_Class)), aes(x = Clay, y = pH, color = Clay_pH_Class)) +\n",
    "  geom_point(size = 2, alpha = 0.7) +\n",
    "  geom_vline(xintercept = c(20, 40), linetype = \"dashed\", color = \"gray50\") +\n",
    "  geom_hline(yintercept = c(6, 7), linetype = \"dashed\", color = \"gray50\") +\n",
    "  labs(\n",
    "    title = \"Clay x pH Segmentation (0-30 cm)\",\n",
    "    x = \"Clay (%)\",\n",
    "    y = \"pH\",\n",
    "    color = \"Class\"\n",
    "  ) +\n",
    "  theme_minimal() +\n",
    "  theme(legend.position = \"right\")\n"
   ]
  }
 ]
}