{ "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": [ "# Session 3 · Soil data preparation with the KSSL dataset (part 2)\n", "**Module 1 · Introduction to R, spatial data and soil data preparation** · SoilFER Training\n", "\n", "[Course page](https://training.yigini.net/modules/01-r-soil-data/) · [Manual chapter](https://training.yigini.net/manual/introduction-to-soil-data-preparation-spatial-data.html) · [Original script](https://github.com/SoilFER/SoilFER-Training-Resources/blob/main/02_scripts/module1/Session3_Soil_Data_Preparation_Part2.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\", \"knitr\", \"readxl\", \"sf\", \"terra\", \"tidyverse\", \"writexl\")\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/module1/Session2_Soil_Data_Preparation_Part1.R\")) { message(\"Running \", s); source(file.path(root, s), echo = FALSE) }\n", "setwd(root)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### SoilFER Online Training Programme — Module 1" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "###############################################################################\n", "# SoilFER Online Training Programme — Module 1\n", "# SESSION 3: Data Preparation with the KSSL Dataset — Part 2 (1.5 hours)\n", "# Sections: Lab data, Duplicates, Harmonization, Spectroscopy prep, Export\n", "###############################################################################\n", "#\n", "# LEARNING OBJECTIVES\n", "# -------------------\n", "# By the end of this session, participants will be able to:\n", "# 1. Attach laboratory (wet chemistry) data to the cleaned site structure\n", "# 2. Validate analytical values against feasible thresholds\n", "# 3. Report and document out-of-bounds issues\n", "# 4. Validate texture data (particle-size fractions)\n", "# 5. Apply targeted corrections or NA replacement for erroneous values\n", "# 6. Detect and resolve duplicated horizons and competing depth sequences\n", "# 7. Harmonize data to standard depth intervals (0–30, 30–60 cm)\n", "# 8. Export cleaned and standardized datasets for DSM and spectroscopy\n", "#\n", "# PREREQUISITE\n", "# ------------\n", "# This session continues directly from Session 2.\n", "#\n", "# The `raw_data` and `site` objects produced in Session 2 must be available\n", "# in the R environment. If they are not available, re-run Session 2 before\n", "# executing this script.\n", "#\n", "# IMPORTANT:\n", "# Session 3 does NOT reconstruct ProfID, coordinates, or horizon depths from\n", "# raw_data. The cleaned `site` object is the structural basis of the workflow.\n", "# Laboratory values are attached to it using the persistent `rowID`.\n", "#\n", "# TIMING GUIDE (approximate)\n", "# ---------------------------\n", "# 0:00 – 0:25 Lab data extraction, joining to site data, threshold validation\n", "# 0:25 – 0:45 Texture validation; targeted correction and NA replacement\n", "# 0:45 – 1:10 Duplicate detection and resolution (average horizons,\n", "# chain_horizons function, surface coverage check)\n", "# 1:10 – 1:30 Depth standardization with aqp::slab(); wide-format output;\n", "# DSM subset; spectroscopy merge; export\n", "###############################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### PART 1 — PREPARING LAB DATA" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# =============================================================================\n", "# PART 1 — PREPARING LAB DATA\n", "# =============================================================================\n", "# NOTE: `raw_data` and `site` from Session 2 must already be in the environment.\n", "#\n", "# NEW WORKFLOW:\n", "# raw_data -> site -> site_lab\n", "#\n", "# The `site` object already contains the cleaned coordinates, sampling dates,\n", "# profile IDs, horizon IDs, and validated depth structure. Therefore these\n", "# steps are not repeated here.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### library(readxl) # Read Excel files" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# =============================================================================\n", "\n", "library(readxl) # Read Excel files\n", "library(tidyverse) # Data manipulation and visualization\n", "library(writexl) # Write Excel files\n", "\n", "# Define the folder to store the results of the exercise\n", "output_dir <- \"03_outputs/module1/\"\n", "\n", "# Define the relative path to the folder with the MIR data\n", "training_dir <- \"01_data/module1/training_data/\"\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 1.1 Extract and Standardize Laboratory Columns\n", "# -----------------------------------------------------------------------------\n", "# NOTE: All analytical parameters must be numeric.\n", "#\n", "# APPROACH:\n", "# 1. Extract analytical data from the original raw_data object\n", "# 2. Preserve rowID as the link to the original observations\n", "# 3. Rename analytical columns to consistent names\n", "# 4. Convert analytical parameters to numeric\n", "# 5. Attach laboratory measurements to the already-cleaned `site` object\n", "#\n", "# This prevents the profile/site cleaning performed in Session 2 from being\n", "# duplicated.\n", "# -----------------------------------------------------------------------------\n", "\n", "# Extract laboratory measurements from the original data.\n", "# rowID links each laboratory record back to the cleaned site record.\n", "lab_data <- raw_data %>%\n", " select(\n", " rowID,\n", " `Estimated Organic Carbon`, `Carbon, Total`, # Soil Organic Carbon and Total Carbon (%)\n", " `Bulk Density, <2mm Fraction, 1/3 Bar`, `Bulk Density, <2mm Fraction, Ovendry`, # Bulk density at 1/3 bar and oven dry (g/cm³)\n", " `Sand, Total`, `Silt, Total`, `Clay`, # Texture (%)\n", " `pH, 1:1 Soil-Water Suspension`, # pH H2O\n", " `CEC, NH4OAc, pH 7.0, 2M KCl displacement`, # CEC in cmol(+)/kg\n", " `Nitrogen, Total`, # Total nitrogen (%)\n", " `Phosphorus, Mehlich3 Extractable`, `Phosphorus, Olsen Extractable`, # Available P (mg/kg)\n", " `Potassium, NH4OAc Extractable, 2M KCl displacement`, # Extractable K (cmol(+)/kg)\n", " `Calcium, NH4OAc Extractable, 2M KCl displacement` # Extractable Ca (cmol(+)/kg)\n", " )\n", "\n", "# Rename laboratory columns to standard, consistent names\n", "names(lab_data) <- c(\n", " \"rowID\",\n", " \"SOC\", # Soil Organic Carbon (%)\n", " \"Carbon_Total\", # Total carbon (%)\n", " \"Bulk.Density_1_3.BAR\", # BD at 1/3 bar (g/cm³)\n", " \"Bulk.Density_ovendry\", # BD oven dry (g/cm³)\n", " \"Sand\", # Sand content (%)\n", " \"Silt\", # Silt content (%)\n", " \"Clay\", # Clay content (%)\n", " \"pH\", # Soil pH (H2O)\n", " \"CEC\", # Cation exchange capacity (cmol(+)/kg)\n", " \"Nitrogen_Total\", # Total nitrogen (%)\n", " \"Phosphorus_Mehlich3\", # Available P (mg/kg)\n", " \"Phosphorus_Olsen\", # Available P (mg/kg)\n", " \"Potassium\", # Exchangeable K (cmol(+)/kg)\n", " \"Calcium\" # Extractable Ca (cmol(+)/kg)\n", ")\n", "\n", "# Ensure numeric type for all analytical parameters\n", "lab_data <- lab_data %>%\n", " mutate(across(-rowID, as.numeric))\n", "\n", "# Attach laboratory measurements to the cleaned site records.\n", "# Only rows retained in `site` are carried forward.\n", "site_lab <- site %>%\n", " left_join(lab_data, by = \"rowID\")\n", "\n", "# Explore the combined site + laboratory dataset\n", "site_lab\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### PART 2 — LABORATORY DATA VALIDATION" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# =============================================================================\n", "# PART 2 — LABORATORY DATA VALIDATION\n", "# WHY: Soil properties have known valid ranges. Values outside these ranges\n", "# may indicate measurement errors, unit mistakes, or data-entry errors.\n", "# All analytical parameters must be numeric.\n", "#\n", "# APPROACH:\n", "# 1. Load thresholds for analytical soil properties\n", "# 2. Find values outside these thresholds\n", "# 3. Generate a detailed report of issues\n", "# 4. Apply corrections where possible\n", "#\n", "# SOURCES for valid ranges:\n", "# NOTE: The analytical thresholds used in this tutorial are based on global\n", "# soil datasets and literature and use the same measurement units as\n", "# the KSSL dataset. Adjust them for your region, soil types, methods,\n", "# and measurement units.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Check 1: Load Property Thresholds and Identify Out-of-Bounds Values" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# =============================================================================\n", "\n", "# -----------------------------------------------------------------------------\n", "# 2.1 Check 1: Load Property Thresholds and Identify Out-of-Bounds Values\n", "# -----------------------------------------------------------------------------\n", "# Analytical thresholds are stored in a CSV file for transparency.\n", "\n", "property_thresholds <- read_csv(\"01_data/module1/kssl/property_thresholds.csv\")\n", "property_thresholds\n", "\n", "# Identify out-of-bounds values; create a list to store issues\n", "out_of_bounds_issues <- list()\n", "\n", "for (i in seq_len(nrow(property_thresholds))) {\n", " prop <- property_thresholds$property[i]\n", " prop_desc <- property_thresholds$description[i]\n", " min_val <- property_thresholds$min_valid[i]\n", " max_val <- property_thresholds$max_valid[i]\n", " \n", " # Check property exists in the dataset\n", " if (prop %in% names(site_lab)) {\n", " x <- site_lab[[prop]]\n", " \n", " # Detect out-of-bounds: non-missing values outside [min_val, max_val]\n", " idx <- which(!is.na(x) & (x < min_val | x > max_val))\n", " \n", " if (length(idx) > 0) {\n", " out_of_bounds_issues[[prop]] <- tibble(\n", " rowID = site_lab$rowID[idx],\n", " property = prop,\n", " description = prop_desc,\n", " value = x[idx],\n", " min_valid = min_val,\n", " max_valid = max_val,\n", " issue = ifelse(\n", " x[idx] < min_val,\n", " paste0(\"Below minimum: \", round(x[idx], 2), \" < \", min_val),\n", " paste0(\"Above maximum: \", round(x[idx], 2), \" > \", max_val)\n", " )\n", " )\n", " }\n", " }\n", "}\n", "\n", "# We can easily detect potential issues with site_lab data\n", "out_of_bounds_issues\n", "\n", "# Remove temporary objects created by the loop, if present\n", "rm(i, max_val, min_val, prop, prop_desc, x)\n", "if (exists(\"idx\")) rm(idx)\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 2.2 Reporting Out-of-Bounds Values and Creating an Audit Trail\n", "# -----------------------------------------------------------------------------\n", "# Export a QC report for review and documentation before correcting data.\n", "\n", "if (length(out_of_bounds_issues) > 0) {\n", " all_issues <- bind_rows(out_of_bounds_issues)\n", " cat(\"\\n Out-of-bounds properties found\\n\")\n", " \n", " # Summary by property\n", " issue_summary <- all_issues %>%\n", " group_by(property, description) %>%\n", " summarise(\n", " count = n(),\n", " min_value_found = min(value, na.rm = TRUE),\n", " max_value_found = max(value, na.rm = TRUE),\n", " min_valid = first(min_valid),\n", " max_valid = first(max_valid),\n", " .groups = \"drop\"\n", " ) %>%\n", " arrange(desc(count))\n", " \n", " cat(\"Issues by property:\\n\")\n", " print(issue_summary)\n", " \n", " # Rows with multiple issues\n", " rows_with_multiple_issues <- all_issues %>%\n", " group_by(rowID) %>%\n", " summarise(\n", " n_issues = n(),\n", " properties = paste(property, collapse = \", \"),\n", " .groups = \"drop\"\n", " ) %>%\n", " filter(n_issues > 1) %>%\n", " arrange(desc(n_issues))\n", " \n", " if (nrow(rows_with_multiple_issues) > 0) {\n", " cat(\"\\n Records with MULTIPLE property issues:\\n\")\n", " print(head(rows_with_multiple_issues, 10))\n", " cat(\"\\nThese records likely have data entry errors and should be reviewed.\\n\")\n", " }\n", " \n", " # Export QC report\n", " write_xlsx(\n", " list(\n", " Summary = issue_summary,\n", " Issues_by_record = rows_with_multiple_issues,\n", " All_issues = all_issues\n", " ),\n", " paste0(output_dir, \"soil_property_validation_report.xlsx\")\n", " )\n", " \n", " cat(\"\\n Detailed report saved to: 03_outputs/module1/soil_property_validation_report.xlsx\\n\")\n", " \n", " rm(all_issues, issue_summary, rows_with_multiple_issues)\n", " \n", "} else {\n", " cat(\"\\n All soil properties within valid ranges!\\n\")\n", "}\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 2.3 Check 2: Texture Validation\n", "# -----------------------------------------------------------------------------\n", "# Particle-size fractions (Clay + Silt + Sand) should sum to approximately 100%.\n", "# Values failing this check are flagged for review — NOT automatically removed.\n", "\n", "texture_problems <- site_lab %>%\n", " mutate(\n", " texture_sum = Clay + Silt + Sand,\n", " texture_valid = abs(texture_sum - 100) < 2\n", " ) %>%\n", " filter(!texture_valid)\n", "\n", "# View texture issues\n", "if (nrow(texture_problems) > 0) {\n", " cat(\" Found\", nrow(texture_problems),\n", " \"records with invalid texture sums\\n\\n\")\n", " print(\n", " texture_problems %>%\n", " select(rowID, ProfID, Clay, Silt, Sand, texture_sum)\n", " )\n", " # Flag for review (do not automatically remove)\n", "} else {\n", " cat(\" Texture problems not found\")\n", "}\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 2.4 Check 3: Correction of Out-of-Bounds Laboratory Values\n", "# -----------------------------------------------------------------------------\n", "# Two options:\n", "# Option 1 (preferred): Targeted correction when the error mechanism is known.\n", "# Use this option if the true value can be recovered from the source or if\n", "# the problem is clearly attributable to an identifiable error.\n", "#\n", "# In this database, SOC is negative in some rows while\n", "# Phosphorus_Mehlich3 contains an extreme value in one row.\n", "# If a Phosphorus_Mehlich3 value is known to use the wrong units\n", "# (e.g. ppb instead of mg/kg), it can be corrected.\n", "#\n", "# Option 2: Replace suspect values with NA when the true value cannot be\n", "# reliably reconstructed.\n", "#\n", "# WARNING:\n", "# - Other datasets may contain errors in different properties.\n", "# - Always inspect potential issues before applying corrections.\n", "# - Only apply automatic corrections when the cause is known and justified.\n", "# -----------------------------------------------------------------------------\n", "\n", "# --- Option 1: Targeted corrections (when error mechanism is known) ---\n", "\n", "# Inspect the nature of each issue before correcting\n", "for (property in names(out_of_bounds_issues)) {\n", " cat(\n", " \"Total errors in\", property, \":\",\n", " n_distinct(out_of_bounds_issues[[property]]$rowID), \"\\n\"\n", " )\n", " print(summary(data.frame(out_of_bounds_issues[property])[4]))\n", "}\n", "\n", "# Correction: Phosphorus Mehlich 3 > 2000 mg/kg\n", "# Example of a likely 1000x unit error (ppb instead of ppm/mg kg-1).\n", "idx <- !is.na(site_lab$Phosphorus_Mehlich3) &\n", " site_lab$Phosphorus_Mehlich3 > 2000\n", "\n", "n_idx <- sum(idx)\n", "\n", "if (n_idx > 0) {\n", " site_lab$Phosphorus_Mehlich3[idx] <-\n", " site_lab$Phosphorus_Mehlich3[idx] / 1000\n", "}\n", "\n", "rm(idx, n_idx)\n", "\n", "\n", "# --- Option 2: Replace out-of-bounds values with NA ---\n", "# Use when the true value cannot be reliably reconstructed.\n", "#\n", "# NOTE:\n", "# The issue list was generated before the targeted correction above.\n", "# Consequently, values corrected under Option 1 should not subsequently be\n", "# replaced with NA. Here the correction is re-checked against the thresholds\n", "# before replacement.\n", "\n", "for (property in names(out_of_bounds_issues)) {\n", " \n", " min_valid <- property_thresholds %>%\n", " filter(.data$property == property) %>%\n", " pull(min_valid)\n", " \n", " max_valid <- property_thresholds %>%\n", " filter(.data$property == property) %>%\n", " pull(max_valid)\n", " \n", " # Replace only values that remain outside the valid range\n", " site_lab <- site_lab %>%\n", " mutate(\n", " \"{property}\" := if_else(\n", " !is.na(.data[[property]]) &\n", " (.data[[property]] < min_valid | .data[[property]] > max_valid),\n", " NA_real_,\n", " .data[[property]]\n", " )\n", " )\n", "}\n", "\n", "if (exists(\"min_valid\")) rm(min_valid)\n", "if (exists(\"max_valid\")) rm(max_valid)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### PART 3 — RESOLVING DUPLICATED DATA IN SOIL PROFILES" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# =============================================================================\n", "# PART 3 — RESOLVING DUPLICATED DATA IN SOIL PROFILES\n", "# =============================================================================\n", "# Repeated or apparently duplicated profiles can arise when:\n", "# - The same location is sampled on different dates\n", "# - Multiple laboratory analyses exist for the same horizon\n", "# - Multiple depth sequences have been associated with one initial profile\n", "# - Identifiers are reused across merged surveys\n", "#\n", "# Because ProfID already includes sampling date, true temporal observations are\n", "# separated before this stage.\n", "#\n", "# Resolution order:\n", "# 1. Average duplicated horizons with identical depth intervals\n", "# 2. Resolve competing depth sequences\n", "# 3. Remove newly separated sequences that do not start at the surface\n", "#\n", "# IMPORTANT:\n", "# Basic coordinate and depth validation was already performed on `site` in\n", "# Session 2 and inherited by `site_lab`. Those checks are NOT repeated here.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Check 1: Detect Potential Horizon Duplicates Within Profiles" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# =============================================================================\n", "\n", "# -----------------------------------------------------------------------------\n", "# 3.1 Check 1: Detect Potential Horizon Duplicates Within Profiles\n", "# -----------------------------------------------------------------------------\n", "\n", "profile_analysis <- site_lab %>%\n", " group_by(ProfID) %>%\n", " summarise(\n", " n_horizons = n(),\n", " n_unique_tops = n_distinct(top),\n", " n_unique_bottoms = n_distinct(bottom),\n", " max_depth = max(bottom, na.rm = TRUE),\n", " .groups = \"drop\"\n", " ) %>%\n", " mutate(\n", " # If all horizons have unique top/bottom values,\n", " # depths are consistent (no repeated intervals)\n", " consistent = (\n", " n_unique_tops == n_horizons &\n", " n_unique_bottoms == n_horizons\n", " ),\n", " likely_duplicates = !consistent\n", " )\n", "\n", "# Find profiles with likely duplicates\n", "duplicates <- profile_analysis %>%\n", " filter(likely_duplicates)\n", "\n", "if (nrow(duplicates) > 0) {\n", " cat(\n", " \" Found\", nrow(duplicates),\n", " \"profiles with likely duplicate measurement sequences\\n\\n\"\n", " )\n", " print(duplicates)\n", "}\n", "\n", "# Select all profiles presenting repeated horizons\n", "duplicates <- site_lab %>%\n", " filter(ProfID %in% duplicates$ProfID)\n", "\n", "# Explore duplicates\n", "duplicates\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 3.2 Resolution Step 1: Average Duplicated Horizons (Same Depth Intervals)\n", "# -----------------------------------------------------------------------------\n", "# When multiple records share identical ProfID + top + bottom, they represent\n", "# repeated measurements of the same layer.\n", "#\n", "# Average numeric analytical properties and retain the first occurrence of\n", "# identifiers and site metadata.\n", "# -----------------------------------------------------------------------------\n", "\n", "site_lab <- site_lab %>%\n", " group_by(ProfID, top, bottom) %>%\n", " summarise(\n", " # Keep identifiers and site metadata as the first value in each group\n", " across(c(rowID, HorID, date, lon, lat), ~ first(.x)),\n", " \n", " # Compute mean for all remaining numeric analytical columns (NA-safe)\n", " across(\n", " where(is.numeric) &\n", " !any_of(c(\"rowID\", \"HorID\", \"lon\", \"lat\", \"top\", \"bottom\")),\n", " ~ if (all(is.na(.x))) NA_real_ else mean(.x, na.rm = TRUE)\n", " ),\n", " .groups = \"drop\"\n", " ) %>%\n", " select(names(site_lab)) # Restore original column order\n", "\n", "site_lab\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 3.3 Resolution Step 2: Resolve ProfID for Multiple Depth Sequences\n", "# -----------------------------------------------------------------------------\n", "# This step separates records that share an initial ProfID but form different,\n", "# non-continuous depth sequences.\n", "#\n", "# The chain_horizons() function identifies consecutive depth chains and assigns\n", "# a numeric suffix to each sequence.\n", "#\n", "# NOTE:\n", "# Missing/negative depths, zero thickness, invalid depth logic, and the initial\n", "# surface-horizon check were already applied to `site` in Session 2.\n", "# -----------------------------------------------------------------------------\n", "\n", "# Create a function to identify sequences of horizons for each profile\n", "chain_horizons <- function(top, bottom) {\n", " n <- length(top)\n", " remaining <- seq_len(n)\n", " chain_id <- integer(n)\n", " cid <- 1\n", " \n", " while (length(remaining) > 0) {\n", " # Start a new chain at the smallest top depth\n", " cur <- remaining[which.min(top[remaining])]\n", " \n", " repeat {\n", " chain_id[cur] <- cid\n", " remaining <- setdiff(remaining, cur)\n", " \n", " # Find the next horizon starting where the current horizon ends\n", " nxt <- remaining[top[remaining] == bottom[cur]]\n", " \n", " if (length(nxt) == 0) break\n", " cur <- nxt[1]\n", " }\n", " \n", " cid <- cid + 1\n", " }\n", " \n", " chain_id\n", "}\n", "\n", "site_lab <- site_lab %>%\n", " group_by(lon, lat, ProfID) %>%\n", " mutate(chain = chain_horizons(top, bottom)) %>% # Detect sequences\n", " arrange(chain, top, .by_group = TRUE) %>% # Sort within each chain\n", " mutate(\n", " ProfID = paste0(ProfID, \"_\", chain) # Add numeric suffix\n", " ) %>%\n", " ungroup()\n", "\n", "if (max(site_lab$chain, na.rm = TRUE) > 1) {\n", " corrected_profiles <- unique(site_lab$ProfID[site_lab$chain >= 2])\n", " \n", " cat(\n", " \"→ Corrected depth continuity in\",\n", " length(corrected_profiles), \"profiles\\n\"\n", " )\n", " \n", " cat(\n", " \" Corrected Profiles:\",\n", " paste(sub(\"_2$\", \"\", corrected_profiles), collapse = \", \"),\n", " \"\\n\"\n", " )\n", "} else {\n", " cat(\"→ No depth continuity corrections were needed\\n\")\n", "}\n", "\n", "# Delete the chain column\n", "site_lab <- site_lab %>%\n", " select(-chain)\n", "\n", "# Delete temporary objects\n", "rm(chain_horizons)\n", "if (exists(\"corrected_profiles\")) rm(corrected_profiles)\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 3.4 Resolution Step 3: Remove Profiles Not Starting at the Surface\n", "# -----------------------------------------------------------------------------\n", "# After competing depth sequences are separated, a new sequence may start below\n", "# 0 cm even though the original ProfID contained a surface horizon.\n", "#\n", "# Therefore the surface check is repeated ONLY for the newly separated profile\n", "# sequences. Keep profiles whose first horizon starts at 0 cm.\n", "# -----------------------------------------------------------------------------\n", "\n", "site_lab <- site_lab %>%\n", " group_by(ProfID) %>%\n", " filter(min(top, na.rm = TRUE) == 0) %>%\n", " arrange(ProfID, top, bottom, HorID) %>%\n", " ungroup()\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 3.5 Export Cleaned Horizon-Level Dataset\n", "# -----------------------------------------------------------------------------\n", "\n", "# Save to CSV\n", "output <- paste0(output_dir, \"KSSL_cleaned.csv\")\n", "write.csv(site_lab, output, row.names = FALSE)\n", "\n", "# Save to Excel\n", "output <- paste0(output_dir, \"KSSL_cleaned.xlsx\")\n", "write_xlsx(site_lab, output)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### PART 4 — STANDARDIZING DATA FOR DSM" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# =============================================================================\n", "# PART 4 — STANDARDIZING DATA FOR DSM\n", "# DSM requires data from one profile per location.\n", "#\n", "# PURPOSE:\n", "# Convert variable-depth horizon data to fixed standard depths\n", "# (0–30 cm and 30–60 cm) for Digital Soil Mapping applications.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.1 Select One Profile per Location (here, Most Complete)" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# =============================================================================\n", "\n", "# -----------------------------------------------------------------------------\n", "# 4.1 Select One Profile per Location (here, Most Complete)\n", "# -----------------------------------------------------------------------------\n", "# When multiple valid profiles exist at the same coordinates, select one using\n", "# an appropriate criterion:\n", "# - Most complete: most horizons (most depth detail) <- used here\n", "# - Best coverage: deepest profile\n", "# - Best quality: fewest missing values\n", "# - Monitoring: profile from period of interest\n", "# -----------------------------------------------------------------------------\n", "\n", "# Keep the most complete profile at each location\n", "horizons <- site_lab %>%\n", " group_by(lon, lat, ProfID) %>%\n", " summarise(\n", " n_hz = n_distinct(paste(top, bottom)),\n", " .groups = \"drop\"\n", " ) %>%\n", " group_by(lon, lat) %>%\n", " dplyr::slice_max(n_hz, n = 1, with_ties = FALSE) %>%\n", " select(lon, lat, ProfID) %>%\n", " inner_join(site_lab, by = c(\"lon\", \"lat\", \"ProfID\")) %>%\n", " ungroup()\n", "\n", "# Explore horizons\n", "horizons\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 4.2 Depth Standardization with aqp::slab()\n", "# -----------------------------------------------------------------------------\n", "# Conceptual explanation: What is slab() doing?\n", "#\n", "# PROBLEM: Profiles have different horizon depths\n", "# Profile 1: 0-10 cm (SOC=3.0%)\n", "# 10-25 cm (SOC=2.5%)\n", "# 25-50 cm (SOC=2.0%)\n", "#\n", "# Profile 2: 0-15 cm (SOC=2.8%)\n", "# 15-40 cm (SOC=2.2%)\n", "# 40-100 cm (SOC=1.5%)\n", "#\n", "# GOAL: Obtain values at standard depths (0-30 cm, 30-60 cm)\n", "#\n", "# DSM requires analytical data at the same depth intervals for every profile.\n", "# The slab() function summarizes horizon data over specified depth intervals.\n", "# -----------------------------------------------------------------------------\n", "\n", "library(aqp)\n", "\n", "# Define standard depth intervals\n", "standard_depths <- c(0, 30, 60) # 0-30, 30-60 cm\n", "\n", "# Select properties to standardize\n", "properties_to_standardize <- c(\n", " \"SOC\",\n", " \"Carbon_Total\",\n", " \"Bulk.Density_1_3.BAR\",\n", " \"Bulk.Density_ovendry\",\n", " \"Sand\",\n", " \"Silt\",\n", " \"Clay\",\n", " \"pH\",\n", " \"CEC\",\n", " \"Nitrogen_Total\",\n", " \"Phosphorus_Mehlich3\",\n", " \"Phosphorus_Olsen\",\n", " \"Potassium\",\n", " \"Calcium\"\n", ")\n", "\n", "# Create SoilProfileCollection object\n", "# aqp needs profile IDs and horizon depth structure\n", "depths(horizons) <- ProfID ~ top + bottom\n", "\n", "# Add spatial information\n", "initSpatial(horizons, crs = \"EPSG:4326\") <- ~ lon + lat\n", "\n", "\n", "# Visual check of the first profiles ===========================================\n", "\n", "# Empty horizons\n", "plotSPC(horizons[1:5])\n", "\n", "# Colour horizons by selected soil properties\n", "plotSPC(horizons[1:5], color = \"SOC\")\n", "plotSPC(horizons[1:5], color = \"pH\")\n", "plotSPC(horizons[1:5], color = \"Clay\")\n", "\n", "# Plot first 10 profiles over the interval 0-30 cm\n", "clods <- profileApply(horizons[1:10], glom, z1 = 0, z2 = 30)\n", "clods <- combine(clods)\n", "\n", "plotSPC(clods, name = \"rowID\", color = \"SOC\")\n", "rect(\n", " xleft = 0.1,\n", " ybottom = 30,\n", " xright = length(horizons[1:10]) + 0.5,\n", " ytop = 0,\n", " border = \"red\",\n", " lty = \"dashed\"\n", ")\n", "\n", "# Density plots of soil parameters\n", "plot(density(horizons$SOC, na.rm = TRUE), main = \"Density plot of SOC\")\n", "plot(density(horizons$pH, na.rm = TRUE), main = \"Density plot of pH\")\n", "\n", "\n", "# Standardize Properties to Fixed Depths ======================================\n", "#\n", "# slab() returns one row per profile, variable, and target depth interval.\n", "# Quantile columns summarize the values contributing to each slab:\n", "# p.q50 = median\n", "# p.q5 / p.q95 = 5th / 95th percentiles\n", "# p.q25 / p.q75 = interquartile range\n", "\n", "# Build the standardization formula\n", "fml <- as.formula(\n", " paste(\"ProfID ~\", paste(properties_to_standardize, collapse = \" + \"))\n", ")\n", "\n", "# View formula\n", "fml\n", "\n", "# Apply slab() to standard depths\n", "KSSL_standardized <- slab(\n", " horizons,\n", " fml,\n", " slab.structure = standard_depths,\n", " na.rm = TRUE\n", ")\n", "\n", "# The output is in long format\n", "KSSL_standardized\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 4.3 Add Percentile-Range Column\n", "# -----------------------------------------------------------------------------\n", "# slab() outputs p.q5, p.q50, and p.q95 among other quantiles.\n", "# Here the 5th–95th percentile range is stored as a compact text field.\n", "# -----------------------------------------------------------------------------\n", "\n", "KSSL_standardized <- KSSL_standardized %>%\n", " mutate(\n", " CI = paste0(\n", " round(p.q5, 3),\n", " \"-\",\n", " round(p.q95, 3)\n", " )\n", " )\n", "\n", "KSSL_standardized\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 4.4 Reshape from Long to Wide Format\n", "# -----------------------------------------------------------------------------\n", "# DSM models typically require one row per profile and depth interval.\n", "# Transform variable names into columns and retain the median plus the\n", "# 5th–95th percentile range.\n", "# -----------------------------------------------------------------------------\n", "\n", "KSSL_standardized <- KSSL_standardized %>%\n", " pivot_wider(\n", " id_cols = c(ProfID, top, bottom),\n", " names_from = variable,\n", " values_from = c(p.q50, CI),\n", " names_glue = \"{variable}_{.value}\"\n", " )\n", "\n", "KSSL_standardized\n", "\n", "# Add geographic coordinates back\n", "KSSL_standardized <- KSSL_standardized %>%\n", " left_join(\n", " site_lab %>%\n", " distinct(ProfID, .keep_all = TRUE) %>%\n", " select(ProfID, lon, lat),\n", " by = \"ProfID\"\n", " ) %>%\n", " relocate(lon, lat, .after = ProfID)\n", "\n", "KSSL_standardized\n", "\n", "# Remove the chain suffix added during depth-sequence resolution.\n", "# At this stage, one selected profile is retained per location.\n", "KSSL_standardized$ProfID <- sub(\"_[0-9]+$\", \"\", KSSL_standardized$ProfID)\n", "\n", "# Result: one row per profile-depth interval with standardized soil properties\n", "head(KSSL_standardized)\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 4.5 Export Standardized Dataset\n", "# -----------------------------------------------------------------------------\n", "\n", "# Save to CSV\n", "output <- paste0(output_dir, \"KSSL_standardized.csv\")\n", "write.csv(KSSL_standardized, output, row.names = FALSE)\n", "\n", "# Save to Excel\n", "output <- paste0(output_dir, \"KSSL_standardized.xlsx\")\n", "write_xlsx(KSSL_standardized, output)\n", "\n", "\n", "# -----------------------------------------------------------------------------\n", "# 4.6 Create DSM Subset (0–30 cm, Key Properties)\n", "# -----------------------------------------------------------------------------\n", "# For Digital Soil Mapping: topsoil (0–30 cm), median estimates (p.q50),\n", "# five properties: Clay, Silt, Sand, SOC, and pH.\n", "\n", "subset_data <- KSSL_standardized %>%\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", "subset_data\n", "\n", "# Save to CSV\n", "output_csv <- paste0(output_dir, \"KSSL_DSM_0-30.csv\")\n", "write.csv(subset_data, output_csv, row.names = FALSE)\n", "cat(\" Saved to:\", output_csv, \"\\n\")\n", "\n", "# Save to Excel\n", "output_xlsx <- paste0(output_dir, \"KSSL_DSM_0-30.xlsx\")\n", "write_xlsx(subset_data, output_xlsx)\n", "cat(\" Saved to:\", output_xlsx, \"\\n\")\n", "\n", "cat(\" Subset data ready for Digital Soil Mapping\\n\")\n", "cat(\" Output file: KSSL_DSM_0-30\\n\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### PART 5 — PREPARING DATA FOR SPECTROSCOPY ANALYSES" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# =============================================================================\n", "# PART 5 — PREPARING DATA FOR SPECTROSCOPY ANALYSES\n", "# =============================================================================\n", "# PURPOSE:\n", "# Create a clean horizon-level dataset with consistent profile structure,\n", "# corrected analytical parameters, and related spectral information for\n", "# estimation of soil properties by spectroscopy.\n", "#\n", "# WHY THIS MATTERS:\n", "# Clean, depth-consistent site and laboratory data improve the reliability of\n", "# relationships developed with spectral measurements.\n", "#\n", "# Merge the cleaned horizon dataset (`site_lab`) with the MIR spectral data.\n", "# Join key: HorID (site_lab) = smp_id (spec).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Read and subset spectral data" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# =============================================================================\n", "\n", "# Read and subset spectral data\n", "spectral_data <- read_excel(\n", " paste0(training_dir, \"/MIR_KANSAS_data.xlsx\"),\n", " sheet = 1\n", ")\n", "\n", "spec <- spectral_data[, -c(1, 3:22)]\n", "\n", "# Merge cleaned site + laboratory data with spectral data by sample identifier\n", "site_lab_spec <- left_join(\n", " site_lab,\n", " spec,\n", " by = c(\"HorID\" = \"smp_id\")\n", ")\n", "\n", "# Save to CSV\n", "output <- paste0(output_dir, \"KSSL_spectral_cleaned.csv\")\n", "write.csv(site_lab_spec, output, row.names = FALSE)\n", "\n", "# Save to Excel\n", "output <- paste0(output_dir, \"KSSL_spectral_cleaned.xlsx\")\n", "write_xlsx(site_lab_spec, output)\n", "\n", "# Remove spectral data object\n", "rm(spec)\n", "\n", "\n", "###############################################################################\n", "# END OF SESSION 3\n", "#\n", "# Summary of exported files:\n", "# KSSL_cleaned — Validated horizon-level dataset\n", "# KSSL_spectral_cleaned — Horizon-level dataset + MIR spectra\n", "# KSSL_standardized — Depth-harmonized (0–30 cm; 30–60 cm)\n", "# KSSL_DSM_0-30 — DSM-ready topsoil dataset (0–30 cm)\n", "# soil_property_validation_report — QC report of out-of-range values\n", "#\n", "# Next session: Session 4 — Spatial Analysis + Covariate Preparation\n", "###############################################################################\n" ] } ] }