{
 "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 4 · Spatial analysis and covariates\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/Session4_Spatial_Analysis_Covariates.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\", \"sf\", \"terra\")\n",
    "need <- setdiff(pkgs, rownames(installed.packages()))\n",
    "if (length(need)) install.packages(need)\n",
    "# packages that are only on GitHub\n",
    "if (!requireNamespace(\"remotes\", quietly = TRUE)) install.packages(\"remotes\")\n",
    "for (r in c(\"ncss-tech/jNSMR\"))\n",
    "  if (!requireNamespace(basename(r), quietly = TRUE)) remotes::install_github(r, upgrade = \"never\")\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": [
    "### SoilFER Online Training Programme — Module 1"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "###############################################################################\n",
    "# SoilFER Online Training Programme — Module 1\n",
    "# SESSION 4: Spatial Analysis + Preparation of Covariates for Sampling Design\n",
    "#            and Digital Soil Mapping (1.5 hours)\n",
    "# Sections: Spatial Analysis, Vector data with sf,Raster data with terra,\n",
    "#           Newhall NSM\n",
    "###############################################################################\n",
    "#\n",
    "# LEARNING OBJECTIVES\n",
    "# -------------------\n",
    "# By the end of this session, participants will be able to:\n",
    "#   1. Understand the difference between vector and raster spatial data\n",
    "#   2. Install and load {sf} and {terra} packages\n",
    "#   3. Create spatial point objects from CSV data\n",
    "#   4. Import, inspect, and export shapefiles and spatial objects\n",
    "#   5. Set and reproject coordinate reference systems (CRS)\n",
    "#   6. Perform geometry operations: buffer, clip, dissolve, spatial join\n",
    "#   7. Visualize spatial data with base R and {ggplot2}\n",
    "#   8. Import, inspect, and visualize raster data\n",
    "#   9. Perform raster operations: mosaic, crop, mask, resample, algebra\n",
    "#  10. Extract raster covariate values at soil sample locations\n",
    "#  11. Run the Newhall Simulation Model (NSM) to derive soil-climate covariates\n",
    "#\n",
    "# PREREQUISITE\n",
    "# ------------\n",
    "# Session 3 outputs must be available. Specifically:\n",
    "#   - KSSL_DSM_0-30.csv (or .xlsx) in 03_outputs/module1/\n",
    "#   - Raster files in 01_data/module2/rasters/ (downloaded from Google Drive)\n",
    "#   - Shapefile Tiger_2020_Counties.shp in 01_data/module1/shapes\n",
    "#\n",
    "# TIMING GUIDE (approximate)\n",
    "# ---------------------------\n",
    "#  0:00 – 0:30  Vector data: sf object creation, CRS, geometry operations,\n",
    "#               buffering, clipping, dissolving, spatial join, visualization\n",
    "#  0:30 – 1:10  Raster data: read/inspect/plot rasters, mosaic, crop, mask,\n",
    "#               resample, combine stacks, raster algebra, extract at points\n",
    "#  1:10 – 1:30  Newhall Simulation Model (NSM): climate averaging,\n",
    "#               input stack preparation, running NSM, exporting results\n",
    "###############################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PART 1 — SPATIAL PACKAGES: INSTALLATION AND LOADING"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# =============================================================================\n",
    "# PART 1 — SPATIAL PACKAGES: INSTALLATION AND LOADING\n",
    "# =============================================================================\n",
    "\n",
    "  install.packages(\"sf\")\n",
    "  install.packages(\"terra\")\n",
    "\n",
    "  library(sf)\n",
    "  library(terra)\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PART 2 — WORKING WITH VECTOR DATA USING {sf}"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# =============================================================================\n",
    "# PART 2 — WORKING WITH VECTOR DATA USING {sf}\n",
    "# =============================================================================\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 2.1  Creating Spatial Vector Objects from CSV Files\n",
    "# -----------------------------------------------------------------------------\n",
    "# Convert the DSM-ready table (with lon/lat) to an sf POINT object.\n",
    "\n",
    "# Import the previously standardized dataset for 0-30 cm depth \n",
    "kdata <- read.csv (paste0(output_dir,\"KSSL_DSM_0-30.csv\"))\n",
    "# Convert to sf: create POINT geometry from lon/lat, set CRS = WGS84 (EPSG:4326)\n",
    "soil_sf <- st_as_sf(kdata, coords = c(\"lon\", \"lat\"), crs = 4326)\n",
    "# Print object geometry\n",
    "soil_sf$geometry\n",
    "# Print CRS details\n",
    "st_crs(soil_sf)\n",
    "\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 2.2  Exporting sf Objects\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Export as shapefile\n",
    "st_write(soil_sf, paste0(output_dir,\"soil_profiles.shp\"), delete_layer = TRUE) # overwrites shp\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 2.3  Importing Shapefiles\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Read the previously created shapefile containing soil profile points\n",
    "soil_sf <- st_read(paste0(output_dir,\"soil_profiles.shp\"))\n",
    "# Inspect the first rows\n",
    "head(soil_sf)\n",
    "# Quick visualization\n",
    "plot(soil_sf[\"pH\"], pch = 16, cex = 0.6, key.pos = 1)  # key.pos controls legend position\n",
    "\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 2.4  Setting and Reprojecting CRS for Vector Data\n",
    "# -----------------------------------------------------------------------------\n",
    "# Always confirm the CRS of spatial data before operations.\n",
    "# st_crs()      → inspect CRS\n",
    "# st_set_crs()  → assign a missing CRS (no coordinate transformation)\n",
    "# st_transform() → reproject to a different CRS (changes coordinate values)\n",
    "\n",
    "sf_3857 <- st_transform(soil_sf, 3857)\n",
    "\n",
    "# Preview of the projected object\n",
    "sf_3857\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 2.5  Inspecting and Manipulating {sf} Objects\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Check structure\n",
    "str(soil_sf)\n",
    "# Summary of attributes + geometry\n",
    "summary(soil_sf)\n",
    "# Extract the geometry column\n",
    "st_geometry(soil_sf)\n",
    "# Check the geometry type. Use head to avoid long printing\n",
    "head(st_geometry_type(soil_sf))\n",
    "# Check the CRS:\n",
    "st_crs(soil_sf)\n",
    "# Check the spatial extent of the object:\n",
    "st_bbox(soil_sf)\n",
    "\n",
    "# Subset sf rows and columns using the pipping (%>%) operator with {dplyr}\n",
    "alkaline <- soil_sf %>%\n",
    "  filter(pH > 7.4) %>%\n",
    "  select(ProfID, pH) %>%\n",
    "  mutate(alkaline = TRUE)\n",
    "# Preview alkaline data\n",
    "alkaline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PART 3 — GEOMETRY OPERATIONS"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# =============================================================================\n",
    "# PART 3 — GEOMETRY OPERATIONS\n",
    "# =============================================================================\n",
    "# NOTE: Distance and area-based operations require a Projected CRS (metres).\n",
    "#       Always ensure both layers share the same CRS before operations.\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 3.1  Buffering Plots (Influence Area)\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Original KSSL data is in WGS84 EPSG:4326  geographic coordinates(lon/lat)\n",
    "# Transform the subset of alkaline soils to NAD83 / UTM 14N (EPSG: 26914)\n",
    "alkaline_utm <- st_transform(alkaline, crs = 26914)  # UTM 14N, distance in metres\n",
    "alkaline_utm\n",
    "\n",
    "# Create a 1 km buffer around each plot\n",
    "plots_buffer_1k <- st_buffer(alkaline_utm, dist = 1000)\n",
    "\n",
    "# Plot buffers (first 3 plots) and then overlay the same points\n",
    "plot(st_geometry(plots_buffer_1k[1:3, ]), \n",
    "     col = rgb(0, 0, 1, 0.3), border = \"blue\",\n",
    "     main = \"Plot of first 3 alkaline soils with 1 km buffers\")\n",
    "# Overlay the original sample locations\n",
    "plot(st_geometry(alkaline_utm[1:3, ]), \n",
    "     add = TRUE, pch = 16, col = \"red\")\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 3.2  Clip Soil Data to Defined Boundaries (Example 1: Riley County)\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Read administrative boundaries (example file)\n",
    "admin <- st_read(\"01_data/module1/shapes/Tiger_2020_Counties.shp\")\n",
    "admin\n",
    "\n",
    "# Ensure both layers share the same CRS\n",
    "admin <- st_transform(admin, st_crs(soil_sf))\n",
    "admin\n",
    "\n",
    "# Example: select one unit (replace with your column/value)\n",
    "study_area <- admin[admin$NAME == \"Riley\", ]\n",
    "# Keep only points inside the study area\n",
    "soil_clip <- st_intersection(soil_sf, study_area)\n",
    "# Plot result\n",
    "plot(st_geometry(study_area), col = NA, border = \"black\", main=\"Sampled soils in Riley County\", cex.main = 1 )\n",
    "plot(st_geometry(soil_clip), add = TRUE, pch = 16, col = \"red\", cex = 0.6)\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 3.3  Dissolve Polygons to Create a Single Study Region Boundary (Example 2)\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Dissolve all counties into one Kansas boundary\n",
    "admin_dissolved <- st_union(admin)\n",
    "# Plot comparison: before vs after dissolve\n",
    "# Adjust graphics settings to 1 row x 2 columns layout with tight margins\n",
    "op <- par(\n",
    "  mfrow = c(1, 2),\n",
    "  mar = c(0.5, 0.5, 0.1, 0.5),  # very small top margin\n",
    "  xaxs = \"i\", yaxs = \"i\"\n",
    ")\n",
    "plot(st_geometry(admin[1]),\n",
    "     col = NA, border = \"black\",\n",
    "     axes = FALSE, asp = 1)\n",
    "title(\"Before dissolve\", line = -0.6, cex.main = 0.95)\n",
    "plot(st_geometry(admin_dissolved),\n",
    "     col = NA, border = \"black\",\n",
    "     axes = FALSE, asp = 1)\n",
    "title(\"After dissolve\", line = -0.6, cex.main = 0.95)\n",
    "# restore graphics settings\n",
    "par(op)\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 3.4  Spatial Relationships and Joins\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Ensure both layers share the same CRS\n",
    "alkaline <- st_transform(alkaline, st_crs(admin))\n",
    "# Join admin attributes to points (points inherit polygon attributes)\n",
    "alkaline_with_admin <- st_join(alkaline, admin, join = st_within)\n",
    "# Preview result\n",
    "head(alkaline_with_admin)\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 3.5  Visualizing {sf} Objects with {ggplot2}\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "library(ggplot2)\n",
    "ggplot(data = soil_sf) +\n",
    "  geom_sf(aes(color = pH), size = 2) +\n",
    "  scale_color_viridis_c() +\n",
    "  theme_minimal() +\n",
    "  labs(\n",
    "    title = \"Soil sampling sites – pH\",\n",
    "    color = \"pH\"\n",
    "  )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PART 4 — WORKING WITH RASTER DATA USING {terra}"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# =============================================================================\n",
    "# PART 4 — WORKING WITH RASTER DATA USING {terra}\n",
    "# =============================================================================\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 4.1  Importing and Exporting Raster Data\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "library(terra)\n",
    "# Read a single-layer raster (e.g., cropland mask or DEM)\n",
    "crops <- rast(paste0(rasters_dir,\"Cropland_Mask_KANSAS.tif\"))\n",
    "# Read a multi-layer raster (e.g., environmental covariates stack)\n",
    "covs <- rast(paste0(rasters_dir,\"Environmental_Covariates_250m_KANSAS.tif\"))\n",
    "# Export a raster as GeoTIFF\n",
    "writeRaster(crops, \"03_outputs/module1/Cropland_Mask_KANSAS.tif\", overwrite = TRUE) # overwrites tif\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 4.2  Inspecting and Exploring SpatRaster Objects\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Basic metadata (prints summary to console)\n",
    "covs\n",
    "# Number of layers\n",
    "nlyr(covs)\n",
    "# Spatial properties\n",
    "ext(covs)\n",
    "res(covs)\n",
    "crs(covs)\n",
    "# Layer names (multi-layer rasters)\n",
    "names(covs)\n",
    "# Value ranges for the first 4 layers\n",
    "global(covs[[1:4]], fun = \"range\", na.rm = TRUE)\n",
    "# Histograms of the first 4 layers\n",
    "# Adjust graphics settings to 2 row x 2 columns layout with tight margins\n",
    "op <- par(\n",
    "  mfrow = c(2, 2),\n",
    "  mar  = c(3.5, 3.5, 2.2, 1.2),  # more space around each panel\n",
    "  oma  = c(1, 1, 1, 1),          # extra space around the whole figure\n",
    "  mgp  = c(2.2, 0.7, 0)          # axis title/labels spacing\n",
    ")\n",
    "\n",
    "for (i in 1:4) {\n",
    "  hist(covs[[i]],\n",
    "       main = paste(\"Histogram of\", names(covs)[i]),\n",
    "       xlab = names(covs)[i])\n",
    "}\n",
    "# restore previous graphics settings\n",
    "par(op)\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 4.3  Raster Visualization\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Plot single-layer crops raster mask\n",
    "plot(crops, main=\"Cropland Mask\")\n",
    "\n",
    "# Plot a multi-layer raster (e.g., environmental covariates stack)\n",
    "plot(covs)\n",
    "\n",
    "# Plot Elevation and Temperature from the covs stack\n",
    "# Adjust graphics settings to 1 row x 2 columns layout with tight margins\n",
    "op <- par(\n",
    "  mfrow = c(1, 2),\n",
    "  mar = c(0.5, 0.5, 0.5, 0.5),\n",
    "  xaxs = \"i\", yaxs = \"i\"\n",
    ")\n",
    "plot(covs[[\"dtm_elevation_250m\"]], main = \"Elevation (DEM)\", cex.main = .8)\n",
    "plot(covs[[\"bio1\"]], main = \"Annual Mean Temperature\", cex.main = .8)\n",
    "par(op) # restore previous graphics settings\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 4.4  Setting and Reprojecting CRS for Raster Data\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Inspect CRS\n",
    "crs(covs)\n",
    "# Reproject raster (example: to EPSG:3857)\n",
    "covs_3857 <- project(covs, \"EPSG:3857\")\n",
    "crs(covs_3857)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PART 5 — SPATIAL OPERATIONS ON RASTERS"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# =============================================================================\n",
    "# PART 5 — SPATIAL OPERATIONS ON RASTERS\n",
    "# =============================================================================\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 5.1  Merge Rasters (Mosaic)\n",
    "# -----------------------------------------------------------------------------\n",
    "# Join adjacent raster tiles into a single continuous layer.\n",
    "\n",
    "# Read input rasters (adjacent tiles)\n",
    "r1 <- rast(paste0(rasters_dir,\"Cropland_Mask_KANSAS-0000000000-0000000000.tif\"))\n",
    "r2 <- rast(paste0(rasters_dir,\"Cropland_Mask_KANSAS-0000000000-0000023296.tif\"))\n",
    "# Mosaic into a single raster\n",
    "m <- mosaic(r1, r2, fun = \"max\")\n",
    "# Cast to integer to optimize file size (not valid for non-integer continuous data)\n",
    "m <- as.int(m)\n",
    "# Write output with compression to reduce file size\n",
    "writeRaster(m, \"03_outputs/module1/Cropland_Mask_KANSAS.tif\", overwrite = TRUE,\n",
    "  wopt = list(datatype = \"INT1U\",  # Byte (0–255); NA retained as NoData\n",
    "    gdal = c(\"COMPRESS=DEFLATE\", \"PREDICTOR=2\", \"ZLEVEL=9\", \"TILED=YES\")\n",
    "  )\n",
    ")\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 5.2  Crop and Mask to Study Area\n",
    "# -----------------------------------------------------------------------------\n",
    "# crop() reduces raster to bounding box of a polygon.\n",
    "# mask() sets values outside the polygon to NA.\n",
    "# Use both together: crop() first, then mask().\n",
    "\n",
    "# Crop and mask covariates to the extent of Riley\n",
    "# Align CRS of both datasets\n",
    "admin <- st_transform(admin, crs(covs))\n",
    "# Crop covariates to the extent of the county\n",
    "covs_crop <- crop(covs, admin[admin$NAME==\"Riley\",])\n",
    "# Mask covariates\n",
    "covs_mask <- mask(covs_crop, admin[admin$NAME==\"Riley\",])\n",
    "\n",
    "# Plot covariate (bio1 = temperature) in Riley County\n",
    "# Adjust graphics settings to 1 row x 2 columns layout with tight margins\n",
    "op <- par(\n",
    "  mfrow = c(1, 2),\n",
    "  mar = c(0.5, 0.5, 0.5, 0.5),\n",
    "  xaxs = \"i\", yaxs = \"i\"\n",
    ")\n",
    "plot(covs_crop[[\"bio1\"]], main= \"Crop\", cex.main = 1)\n",
    "plot(covs_mask[[\"bio1\"]], main= \"Mask\", cex.main = 1)\n",
    "# restore previous graphics settings\n",
    "par(op) \n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 5.3  Resampling and Alignment\n",
    "# -----------------------------------------------------------------------------\n",
    "# resample() transfers values from one SpatRaster to the grid of another.\n",
    "# Use method = \"near\" for categorical rasters; \"bilinear\" for continuous.\n",
    "\n",
    "# Align raster of crops to the grid of covariates\n",
    "crops_aligned <- resample(crops, covs, method = \"near\")\n",
    "# Spatial properties of aligned crops and covariates are equal\n",
    "ext(covs)\n",
    "ext(crops_aligned)\n",
    "res(covs)\n",
    "res(crops_aligned)\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 5.4  Combine and Remove Rasters in a Raster Stack\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Set a clear layer name for the aligned cropland mask\n",
    "names(crops_aligned) <- \"crops\"\n",
    "# Add the cropland layer to the covariate stack\n",
    "covs <- c(covs,crops_aligned)\n",
    "# Check layer names\n",
    "names(covs)\n",
    "\n",
    "# Remove \"crops\" from the stack\n",
    "covs <- covs[[names(covs) != \"crops\"]]\n",
    "# Check layer names\n",
    "names(covs)\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 5.5  Raster Algebra\n",
    "# -----------------------------------------------------------------------------\n",
    "# Apply cell-by-cell mathematical or logical operations on rasters.\n",
    "# Rasters must share the same CRS, origin, resolution, and extent.\n",
    "\n",
    "# Example: classify a continuous covariate into two classes\n",
    "high_temperature <- covs[[1]] > 12\n",
    "# Adjust graphics settings to 1 row x 2 columns layout with tight margins\n",
    "op <- par(\n",
    "  mfrow = c(1, 2),\n",
    "  mar = c(0.5, 0.5, 0.5, 0.5),\n",
    "  xaxs = \"i\", yaxs = \"i\"\n",
    ")\n",
    "plot(covs[[\"bio1\"]], main= \"Annual Mean Temperature\", cex.main = .8)\n",
    "plot(high_temperature, main= \"High Temperature (>12ºC)\", cex.main = .8)\n",
    "par(op) # restore previous graphics settings\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 5.6  Merging GEE-Exported Tiles (Practical Example)\n",
    "# -----------------------------------------------------------------------------\n",
    "# Google Earth Engine exports large rasters as tiles.\n",
    "# Merge Crop mask rasters (two adjacent tiles)\n",
    "# Read input rasters\n",
    "r1 <- rast(paste0(rasters_dir,\"Cropland_Mask_KANSAS-0000000000-0000000000.tif\"))\n",
    "r2 <- rast(paste0(rasters_dir,\"Cropland_Mask_KANSAS-0000000000-0000023296.tif\"))\n",
    "# Mosaic into a single raster\n",
    "m <- mosaic(r1, r2, fun = \"max\")\n",
    "# Cast to integer to optimize file size\n",
    "# This ensures the output is written as integer values (1 and NA),\n",
    "m <- as.int(m)\n",
    "# Write output with compression to reduce file size\n",
    "writeRaster(m, \"03_outputs/module1/Cropland_Mask_KANSAS.tif\", overwrite = TRUE,\n",
    "  wopt = list(datatype = \"INT1U\",  # Byte (0–255); NA retained as NoData\n",
    "    gdal = c(\"COMPRESS=DEFLATE\", \"PREDICTOR=2\", \"ZLEVEL=9\", \"TILED=YES\")\n",
    "  )\n",
    ")\n",
    "\n",
    "# Merge ESA rasters (two adjacent tiles)\n",
    "# Read input rasters\n",
    "r1 <- rast(paste0(rasters_dir,\"ESA_WorldCover_2021_KANSAS-0000000000-0000000000.tif\"))\n",
    "r2 <- rast(paste0(rasters_dir,\"ESA_WorldCover_2021_KANSAS-0000000000-0000065536.tif\"))\n",
    "# Mosaic into a single raster\n",
    "m <- mosaic(r1, r2, fun = \"max\")\n",
    "# Cast to integer to optimize file size\n",
    "# This ensures the output is written as integer values\n",
    "m <- as.int(m)\n",
    "# Write output with compression to reduce file size\n",
    "writeRaster(m, \"03_outputs/module1/ESA_WorldCover_2021_KANSAS.tif\", overwrite = TRUE,\n",
    "  wopt = list(datatype = \"INT1U\",  # Byte (0–255); NA retained as NoData\n",
    "    gdal = c(\"COMPRESS=DEFLATE\", \"PREDICTOR=2\", \"ZLEVEL=9\", \"TILED=YES\")\n",
    "  )\n",
    ")\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 5.7  Combining Environmental (Raster) and Soil (Vector) Information\n",
    "# -----------------------------------------------------------------------------\n",
    "# Extract raster values at sample locations to link soil observations\n",
    "# (points) with environmental covariates (rasters) for DSM calibration.\n",
    "\n",
    "# Extract covariate values at point locations (soils_sf)\n",
    "cov_values <- extract(covs, soil_sf)\n",
    "cov_values\n",
    "\n",
    "# cov_values includes an ID column in the first column\n",
    "soil_sf <- cbind(soil_sf, cov_values[ , -1, drop = FALSE])\n",
    "soil_sf\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PART 6 — NEWHALL SIMULATION MODEL (NSM) FOR SOIL-CLIMATE COVARIATES"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# =============================================================================\n",
    "# PART 6 — NEWHALL SIMULATION MODEL (NSM) FOR SOIL-CLIMATE COVARIATES\n",
    "# =============================================================================\n",
    "# The NSM estimates soil moisture and temperature regimes from climate and\n",
    "# site characteristics. Outputs are useful as environmental covariates in\n",
    "# DSM and sampling design workflows.\n",
    "#\n",
    "# Required inputs: monthly precipitation (pJan…pDec), monthly temperature\n",
    "#   (tJan…tDec), elevation (DEM), available water capacity (AWC), lon/lat.\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 6.1  Environment Setup\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Packages\n",
    "library(terra)\n",
    "library(sf)\n",
    "library(dplyr)\n",
    "library(jNSMR)\n",
    "\n",
    "# Target CRS for the workflow (use a projected CRS for analysis)\n",
    "epsg <- \"EPSG:3857\"   # Web Mercator (meters). Replace if needed.\n",
    "agg.factor <- 1       # Optional aggregation to speed up Newhall calculations (1 = no aggregation)\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 6.2  Prepare Monthly Average Rasters of Climatic Data\n",
    "# -----------------------------------------------------------------------------\n",
    "# TerraClimate files contain large time series (1981–2020). Monthly\n",
    "# temperature and precipitation averages must be calculated as independent\n",
    "# raster files. Temperature is stored scaled (×10) and must be adjusted.\n",
    "\n",
    "# Monthly temperature time series (multi-year, monthly layers)\n",
    "tmp <- rast(paste0(rasters_dir,\"TerraClimate_AvgTemp_1981_2023_KANSAS.tif\"))\n",
    "\n",
    "# Number of years in the stack (assuming complete years)\n",
    "nyears <- nlyr(tmp) / 12\n",
    "\n",
    "# Monthly means across all years\n",
    "temp_list <- vector(\"list\", 12)\n",
    "tmp_ref <- tmp[[1]]  # reference layer for initializing\n",
    "for(i in 1:12){\n",
    "  var_sum <- tmp_ref * 0\n",
    "  k <- i\n",
    "  for(j in 1:nyears){\n",
    "    var_sum <- var_sum + tmp[[k]]\n",
    "    k <- k + 12\n",
    "  }\n",
    "  temp_list[[i]] <- var_sum / nyears\n",
    "}\n",
    "# Convert to stack\n",
    "Temp_Stack <- rast(temp_list) * 0.1  # convert scaled values to °C (×0.1)\n",
    "# Write to disk\n",
    "writeRaster(Temp_Stack, \"03_outputs/module1/Temp_Stack_monthlyMean.tif\", overwrite = TRUE)\n",
    "\n",
    "\n",
    "# Monthly precipitation time series (multi-year, monthly layers)\n",
    "prec_mm <- rast(paste0(rasters_dir,\"TerraClimate_Precip_1981_2023_KANSAS.tif\"))\n",
    "# Number of years in the stack (assuming complete years\n",
    "nyears <- nlyr(prec_mm) / 12\n",
    "# Monthly means across all years\n",
    "prec_list <- vector(\"list\", 12)\n",
    "pre_ref <- prec_mm[[1]] # reference layer for initializing\n",
    "for(i in 1:12){\n",
    "  var_sum <- pre_ref * 0\n",
    "  k <- i\n",
    "  for(j in 1:nyears){\n",
    "    var_sum <- var_sum + prec_mm[[k]]\n",
    "    k <- k + 12\n",
    "  }\n",
    "  prec_list[[i]] <- var_sum / nyears\n",
    "}\n",
    "# Convert to stack\n",
    "Prec_Stack <- rast(prec_list)\n",
    "# Write to disk\n",
    "writeRaster(Prec_Stack, \"03_outputs/module1/Prec_Stack_monthlyMean.tif\", overwrite = TRUE)\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 6.3  Prepare Monthly Average Stack of Climatic Data\n",
    "# -----------------------------------------------------------------------------\n",
    "# NSM expects monthly layers with specific, predefined names.\n",
    "\n",
    "# Load outputs (or reuse objects from above)\n",
    "Prec_Stack <- rast(\"03_outputs/module1/Prec_Stack_monthlyMean.tif\")\n",
    "Temp_Stack <- rast(\"03_outputs/module1/Temp_Stack_monthlyMean.tif\")\n",
    "\n",
    "names(Prec_Stack) <- c(\"pJan\",\"pFeb\",\"pMar\",\"pApr\",\"pMay\",\"pJun\",\"pJul\",\"pAug\",\"pSep\",\"pOct\",\"pNov\",\"pDec\")\n",
    "names(Temp_Stack) <- c(\"tJan\",\"tFeb\",\"tMar\",\"tApr\",\"tMay\",\"tJun\",\"tJul\",\"tAug\",\"tSep\",\"tOct\",\"tNov\",\"tDec\")\n",
    "\n",
    "combined_stack <- c(Prec_Stack, Temp_Stack)\n",
    "writeRaster(combined_stack, \"03_outputs/module1/climate_vars.tif\", overwrite = TRUE)\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 6.4  Prepare Site Descriptors (DEM, AWC, Lon/Lat)\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# DEM\n",
    "# Read raster\n",
    "elev <- rast(paste0(rasters_dir,\"/DEM_KANSAS.tif\"))\n",
    "# Align CRS to the common EPSG\n",
    "if(crs(elev) != epsg) elev <- project(elev, epsg, method = \"near\")\n",
    "\n",
    "# AWC (0–200 cm, in mm)\n",
    "# Read raster\n",
    "awc <- rast(paste0(rasters_dir,\"awc_KANSAS.tif\"))\n",
    "# Align CRS to the common EPSG\n",
    "if(crs(awc) != epsg) awc <- project(awc, epsg, method = \"near\")\n",
    "# Align geometries with the DEM raster\n",
    "awc <- resample(awc, elev)\n",
    "# Define name\n",
    "names(awc) <- \"awc\"\n",
    "\n",
    "# Temperature and Precipitation\n",
    "# Read climate stack\n",
    "newhall <- rast(\"03_outputs/module1/climate_vars.tif\")\n",
    "# Align CRS to the common EPSG\n",
    "if(crs(newhall) != epsg) newhall <- project(newhall, epsg, method = \"near\")\n",
    "# Align geometries with the DEM raster\n",
    "newhall <- resample(newhall, elev)\n",
    "\n",
    "# Longitude and Latitude \n",
    "# calculate lon/lat from the climate stack\n",
    "newhall$lonDD <- init(newhall[[1]], \"x\")\n",
    "newhall$latDD <- init(newhall[[1]], \"y\")\n",
    "\n",
    "# Mask climate stack to DEM valid area and add DEM/AWC to the stack\n",
    "newhall <- mask(newhall, elev)\n",
    "newhall$awc  <- awc\n",
    "newhall$elev <- elev\n",
    "\n",
    "# Save full NSM input stack\n",
    "writeRaster(newhall, \"03_outputs/module1/climate_newhall_vars.tif\", overwrite = TRUE)\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 6.5  Calculate NSM Outputs\n",
    "# -----------------------------------------------------------------------------\n",
    "# NOTE: This calculation can take hours at high resolution.\n",
    "# Optional aggregation (agg.factor > 1) reduces runtime.\n",
    "\n",
    "# Read the full Newhall input stack\n",
    "newhall <- rast(\"03_outputs/module1/climate_newhall_vars.tif\")\n",
    "\n",
    "# Optional: aggregate to speed up (e.g., fact = 2, 3, ...)\n",
    "# newhall <- aggregate(newhall, fact = agg.factor)\n",
    "\n",
    "system.time({\n",
    "  newhall_results <- jNSMR::newhall_batch(newhall, cores = 4)\n",
    "})\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# 6.6  Inspect and Export NSM Results\n",
    "# -----------------------------------------------------------------------------\n",
    "\n",
    "# Build a validity mask from the precipitation layer\n",
    "mask_valid <- !is.na(newhall_results$annualRainfall)\n",
    "# Apply the mask to all NSM layers\n",
    "results <- mask(newhall_results, mask_valid)\n",
    "# Quick inspection\n",
    "# -------------------------------------------------------------------\n",
    "names(results)\n",
    "plot(results)\n",
    "\n",
    "newhall_layers <- tribble(\n",
    "  ~layer, ~units, ~description,\n",
    "  \"annualRainfall\", \"mm/yr\", \"Total precipitation accumulated over the year.\",\n",
    "  \"waterHoldingCapacity\", \"mm\", \"Soil water-holding capacity used by the model (plant-available storage over the modeled profile/control section).\",\n",
    "  \"annualWaterBalance\", \"mm/yr\", \"Net annual water balance (precipitation − potential evapotranspiration) summed over the year; positive = surplus, negative = deficit.\",\n",
    "  \"annualPotentialEvapotranspiration\", \"mm/yr\", \"Total annual potential evapotranspiration (atmospheric evaporative demand) summed over the year.\",\n",
    "  \"summerWaterBalance\", \"mm (summer period)\", \"Net water balance during the model’s summer period (precipitation − potential evapotranspiration) for that season.\",\n",
    "  \"dryDaysAfterSummerSolstice\", \"days\", \"Number of days classified as dry after the summer solstice (indicator of summer dryness duration).\",\n",
    "  \"moistDaysAfterWinterSolstice\", \"days\", \"Number of days classified as moist after the winter solstice (indicator of post-winter moisture duration).\",\n",
    "  \"numCumulativeDaysDry\", \"days\", \"Cumulative number of dry days across the year (sum of days meeting the model’s dry condition).\",\n",
    "  \"numCumulativeDaysMoistDry\", \"days\", \"Cumulative number of moist-dry (intermediate) days across the year.\",\n",
    "  \"numCumulativeDaysMoist\", \"days\", \"Cumulative number of moist days across the year.\",\n",
    "  \"numCumulativeDaysDryOver5C\", \"days\", \"Cumulative number of dry days when temperature is > 5 °C.\",\n",
    "  \"numCumulativeDaysMoistDryOver5C\", \"days\", \"Cumulative number of moist-dry days when temperature is > 5 °C.\",\n",
    "  \"numCumulativeDaysMoistOver5C\", \"days\", \"Cumulative number of moist days when temperature is > 5 °C.\",\n",
    "  \"numConsecutiveDaysMoistInSomeParts\", \"days\", \"Consecutive days when the soil is moist in some part of the control section (persistence of partial-profile moisture).\",\n",
    "  \"numConsecutiveDaysMoistInSomePartsOver8C\", \"days\", \"Consecutive days when the soil is moist in some part of the control section and temperature is > 8 °C.\",\n",
    "  \"temperatureRegime\", \"class\", \"Categorical soil temperature regime output (e.g., mesic, thermic), derived from the model’s temperature criteria.\",\n",
    "  \"moistureRegime\", \"class\", \"Categorical soil moisture regime output (e.g., udic, ustic, xeric), derived from the model’s moisture criteria.\",\n",
    "  \"regimeSubdivision1\", \"class\", \"Van Wambeke modifier (Part 1): a simple qualifier such as Wet/Dry/Typic/Weak/Extreme that refines the moisture regime.\",\n",
    "  \"regimeSubdivision2\", \"class\", \"Van Wambeke base term (Part 2): the main regime term (often Temp* or Trop* forms like Tempustic/Tropustic or Tempudic/Tropudic). Combine with Part 1 to get the full label (e.g., 'Wet Tempustic').\"\n",
    ")\n",
    "newhall_layers\n",
    "\n",
    "# Export to GeoTIFF\n",
    "writeRaster(results, \"03_outputs/module1/newhall.tif\", overwrite = TRUE)\n",
    "\n",
    "# Optional: compact integer exports\n",
    "#    Multiply by 10 to preserve 1 decimal place, then store as integer.\n",
    "#    (Useful to reduce file size when decimals are not critical.)\n",
    "results_intx10 <- app(results, function(x) as.integer(round(x * 10)))\n",
    "writeRaster(results_intx10, \"03_outputs/module1/newhall_intx10.tif\", overwrite = TRUE)\n",
    "\n",
    "\n",
    "# Optional: \n",
    "## NSM calculation on a small area\n",
    "admin_3857 <- st_transform(admin, 3857)\n",
    "admin_3857 <- vect(admin_3857)            # convert sf -> terra vector\n",
    "riley <- admin_3857[admin_3857$NAME==\"Riley\"] # Select the extent of Riley county\n",
    "\n",
    "# Make sure CRS match (project polygon to raster CRS if needed)\n",
    "if (!same.crs(newhall, riley)) {\n",
    "  admin_3857 <- project(riley, crs(r))\n",
    "}\n",
    "r_crop <- crop(newhall, riley)       # trims extent\n",
    "system.time({\n",
    "  newhall_results <- jNSMR::newhall_batch(r_crop, cores = 4)\n",
    "})\n",
    "\n",
    "names(newhall_results)\n",
    "plot(newhall_results)\n",
    "plot(newhall_results[[16:19]])\n",
    "\n",
    "\n",
    "\n",
    "###############################################################################\n",
    "# END OF SESSION 4 — END OF MODULE 1\n",
    "#\n",
    "# The raster stack 'newhall_results' contains NSM outputs (soil moisture and\n",
    "# temperature regime classes and indices) ready to use as covariates in:\n",
    "#   - Module 2: Sampling Design (covariate space coverage)\n",
    "#   - Module 4: Digital Soil Mapping (environmental predictors)\n",
    "#\n",
    "# Summary of Module 1 sessions:\n",
    "#   Session 1 — R fundamentals (objects, structures, functions, tidyverse)\n",
    "#   Session 2 — KSSL data import; site/coord/depth cleaning (Part 1)\n",
    "#   Session 3 — Lab validation; duplicates; depth harmonization (Part 2)\n",
    "#   Session 4 — Spatial analysis with sf/terra; covariate preparation; NSM\n",
    "###############################################################################\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"
   ]
  }
 ]
}