{
 "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": [
    "# Step 1 · Climate data and soil-climate covariates (Newhall model)\n",
    "**Module 2 · Soil sampling design** · SoilFER Training\n",
    "\n",
    "[Course page](https://training.yigini.net/modules/02-sampling-design/) · [Manual chapter](https://training.yigini.net/manual/sampling-design-for-soil-surveys.html) · [Original script](https://github.com/SoilFER/SoilFER-Training-Resources/blob/main/02_scripts/module2/1_climate_soil_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 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\", \"sf\", \"sp\", \"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",
    "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": [
    "### Google Earth Engine exports\n",
    "Module 2 starts from rasters exported with the Earth Engine script `02_scripts/module2/gee_soilfer_env_covariates.txt` (run it in the [GEE Code Editor](https://code.earthengine.google.com)). The exports land in your Google Drive. Share that Drive folder as *Anyone with the link → Viewer*, paste the link below and run the cell to bring the rasters into this notebook."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "gee_folder <- \"\"   # <- paste your Google Drive folder link here\n",
    "if (nzchar(gee_folder)) {\n",
    "  dest <- file.path(root, \"02_scripts/rasters/GEE_Exports\")   # where 1_climate_soil_covariates.R looks\n",
    "  dir.create(dest, recursive = TRUE, showWarnings = FALSE)\n",
    "  system(\"python3 -m pip -q install gdown\")\n",
    "  system(paste(\"python3 -m gdown --folder --quiet\", shQuote(gee_folder), \"-O\", shQuote(dest)))\n",
    "  dir.create(file.path(root, \"01_data/module2/rasters\"), recursive = TRUE, showWarnings = FALSE)\n",
    "  file.copy(list.files(dest, full.names = TRUE, recursive = TRUE), file.path(root, \"01_data/module2/rasters\"), overwrite = FALSE)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Run time:** the complete sampling design for a whole country runs for many hours, longer than a Colab session allows. In Colab use the Kansas example or a single region, or run the full design on a desktop computer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Climate Data Processing and Soil-Climate Environmental Covariates"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################################################################################\n",
    "# Climate Data Processing and Soil-Climate Environmental Covariates\n",
    "################################################################################\n",
    "# \n",
    "# PURPOSE:\n",
    "# This script processes raw temperature and precipitation data from Google Earth\n",
    "# Engine (GEE), calculates monthly averages, and generates soil-climate \n",
    "# environmental covariates using the Newhall model.\n",
    "#\n",
    "# WORKFLOW:\n",
    "# 1. Process raw temperature and precipitation data from GEE\n",
    "# 2. Calculate monthly averages across all years\n",
    "# 3. Combine climate data with elevation and soil water capacity\n",
    "# 4. Run Newhall soil-climate model\n",
    "# 5. Generate and export results\n",
    "#\n",
    "# REQUIREMENTS:\n",
    "# - Temperature and precipitation data from GEE (see soilfer_env_covariates script)\n",
    "# - Digital elevation model (DEM)\n",
    "# - Available water capacity (AWC) data (0-200 cm depth)\n",
    "# - Region of interest (ROI) shapefile\n",
    "#\n",
    "# AUTHORS: Luis Rodriguez-Lado, PhD \n",
    "#          Wanderson de Sousa Mendes, PhD\n",
    "# DATE: 29 October 2025\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0 - ENVIRONMENT SETUP"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 0 - ENVIRONMENT SETUP =======================================================\n",
    "\n",
    "# Set working directory to source file location\n",
    "# Note: This requires running the script in RStudio\n",
    "setwd(file.path(Sys.getenv(\"SOILFER_ROOT\"), \"02_scripts/module2\"))  # folder of this script (was rstudioapi::getActiveDocumentContext, RStudio only)\n",
    "setwd(\"../\")  # Move up to main project folder\n",
    "getwd()       # Display current working directory\n",
    "\n",
    "# List of required packages\n",
    "packages <- c(\"sp\",         # Spatial data classes and methods\n",
    "              \"terra\",      # Raster data manipulation\n",
    "              \"sf\",         # Simple features for vector data\n",
    "              \"jNSMR\",      # Newhall Simulation Model in R\n",
    "              \"dplyr\")      # Data manipulation\n",
    "\n",
    "# Load all packages\n",
    "invisible(lapply(packages, library, character.only = TRUE))\n",
    "rm(packages)  # Clean up environment\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1 - USER-DEFINED VARIABLES"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 1 - USER-DEFINED VARIABLES ==================================================\n",
    "\n",
    "# Define file paths (modify these according to your directory structure)\n",
    "raster.path <- \"rasters/\"           # Path to raster data\n",
    "shp.path <- \"shapes/\"               # Path to shapefiles\n",
    "\n",
    "\n",
    "# Define coordinate reference system\n",
    "# EPSG:3857 is Web Mercator projection (units in meters)\n",
    "# Verify on https://epsg.io/\n",
    "epsg <- \"EPSG:3857\"\n",
    "\n",
    "# Aggregation factor for upscaling rasters (if needed)\n",
    "# Set to 1 for no aggregation, >1 to reduce resolution\n",
    "agg.factor <- 1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 - PROCESS TEMPERATURE DATA"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 2 - PROCESS TEMPERATURE DATA ================================================\n",
    "# Temperature data should be downloaded from GEE script\n",
    "# File format: TerraClimate_AvgTemp_1981_2023_KANSAS.tif (monthly data for multiple years)\n",
    "\n",
    "# Load raw temperature data (scaled by 10, needs conversion to Celsius)\n",
    "tmp <- terra::rast(paste0(raster.path, \"GEE_Exports/TerraClimate_AvgTemp_1981_2023_KANSAS.tif\"))\n",
    "\n",
    "# Get reference layer (first month) for initialization\n",
    "tmp_Jan_1 <- tmp[[1]]\n",
    "cat(\"Temperature data dimensions:\", dim(tmp), \"\\n\")\n",
    "cat(\"Total months in dataset:\", dim(tmp)[3], \"\\n\")\n",
    "cat(\"Number of years:\", dim(tmp)[3]/12, \"\\n\")\n",
    "\n",
    "# Create empty list to store monthly averages\n",
    "temp_list <- list()\n",
    "\n",
    "# Calculate average temperature for each month across all years\n",
    "# Loop through 12 months\n",
    "for (i in 1:12) { \n",
    "  \n",
    "  # Initialize sum raster with zeros\n",
    "  var_sum <- tmp_Jan_1 * 0\n",
    "  \n",
    "  # Start with month i\n",
    "  k <- i\n",
    "  \n",
    "  # Loop through all years\n",
    "  for (j in 1:(dim(tmp)[3]/12)) {\n",
    "    cat(\"Processing temperature month:\", k, \"\\n\")\n",
    "    \n",
    "    # Add current month to sum\n",
    "    var_sum <- var_sum + tmp[[k]]\n",
    "    \n",
    "    # Move to same month next year\n",
    "    k <- k + 12\n",
    "  }\n",
    "  \n",
    "  # Calculate average by dividing sum by number of years\n",
    "  var_avg <- var_sum / (dim(tmp)[3]/12)\n",
    "  \n",
    "  # Store in list\n",
    "  temp_list[[i]] <- var_avg\n",
    "}\n",
    "\n",
    "# Create raster stack from list\n",
    "Temp_Stack <- terra::rast(temp_list)\n",
    "\n",
    "# Convert from scaled values to degrees Celsius (data is scaled by 10)\n",
    "Temp_Stack <- Temp_Stack * 0.1\n",
    "\n",
    "# Save intermediate temperature result\n",
    "terra::writeRaster(Temp_Stack, \n",
    "                   filename = paste0(raster.path, 'Temp_Stack_01-22_TC.tif'),\n",
    "                   overwrite = TRUE)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3 - PROCESS PRECIPITATION DATA"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 3 - PROCESS PRECIPITATION DATA ==============================================\n",
    "# Precipitation data should be downloaded from GEE \n",
    "# File format: TerraClimate_Precip_1981_2023_KANSAS.tif (monthly data for multiple years)\n",
    "\n",
    "# Load raw precipitation data (mm)\n",
    "prec_mm <- terra::rast(paste0(raster.path, \"/GEE_Exports/TerraClimate_Precip_1981_2023_KANSAS.tif\"))\n",
    "\n",
    "# Get reference layer (first month) for initialization\n",
    "pre_Jan_1 <- prec_mm[[1]]\n",
    "cat(\"Precipitation data dimensions:\", dim(prec_mm), \"\\n\")\n",
    "\n",
    "# Create empty list to store monthly averages\n",
    "prec_list <- list()\n",
    "\n",
    "# Calculate average precipitation for each month across all years\n",
    "# Loop through 12 months\n",
    "for (i in 1:12) { \n",
    "  \n",
    "  # Initialize sum raster with zeros\n",
    "  var_sum <- pre_Jan_1 * 0\n",
    "  \n",
    "  # Start with month i\n",
    "  k <- i\n",
    "  \n",
    "  # Loop through all years\n",
    "  for (j in 1:(dim(prec_mm)[3]/12)) {\n",
    "    cat(\"Processing precipitation month:\", k, \"\\n\")\n",
    "    \n",
    "    # Add current month to sum\n",
    "    var_sum <- var_sum + prec_mm[[k]]\n",
    "    \n",
    "    # Move to same month next year\n",
    "    k <- k + 12\n",
    "  }\n",
    "  \n",
    "  # Calculate average by dividing sum by number of years\n",
    "  var_avg <- var_sum / (dim(prec_mm)[3]/12)\n",
    "  \n",
    "  # Store in list\n",
    "  prec_list[[i]] <- var_avg\n",
    "}\n",
    "\n",
    "# Create raster stack from list\n",
    "Prec_Stack <- terra::rast(prec_list)\n",
    "\n",
    "# Save intermediate precipitation result\n",
    "terra::writeRaster(Prec_Stack, \n",
    "                   filename = paste0(raster.path, 'Prec_Stack_01-22_TC.tif'),\n",
    "                   overwrite = TRUE)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4 - MERGE CLIMATE DATA"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 4 - MERGE CLIMATE DATA ======================================================\n",
    "# Combine temperature and precipitation into single climate variables file\n",
    "\n",
    "# Define layer names following Newhall model conventions\n",
    "# p = precipitation (mm), t = temperature (°C)\n",
    "prec_colnames <- c(\"pJan\", \"pFeb\", \"pMar\", \"pApr\", \"pMay\", \"pJun\", \n",
    "                   \"pJul\", \"pAug\", \"pSep\", \"pOct\", \"pNov\", \"pDec\")\n",
    "temp_colnames <- c(\"tJan\", \"tFeb\", \"tMar\", \"tApr\", \"tMay\", \"tJun\", \n",
    "                   \"tJul\", \"tAug\", \"tSep\", \"tOct\", \"tNov\", \"tDec\")\n",
    "\n",
    "# Rename precipitation layers\n",
    "names(Prec_Stack) <- prec_colnames\n",
    "\n",
    "# Rename temperature layers\n",
    "names(Temp_Stack) <- temp_colnames\n",
    "\n",
    "# Combine into single raster stack (24 layers total)\n",
    "combined_stack <- c(Prec_Stack, Temp_Stack)\n",
    "\n",
    "print(names(combined_stack))\n",
    "\n",
    "# Export merged climate data (primary input for Newhall model)\n",
    "terra::writeRaster(combined_stack, \n",
    "                   filename = paste0(raster.path, 'climate_vars.tif'),\n",
    "                   overwrite = TRUE)\n",
    "\n",
    "# Clean up temporary variables\n",
    "rm(tmp, prec_mm, temp_list, prec_list, tmp_Jan_1, pre_Jan_1, \n",
    "   Temp_Stack, Prec_Stack, prec_colnames, temp_colnames)\n",
    "gc()  # Garbage collection to free memory\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5 - IMPORT ADDITIONAL SPATIAL DATA"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 5 - IMPORT ADDITIONAL SPATIAL DATA ==========================================\n",
    "\n",
    "# Load region of interest (ROI) boundary shapefile\n",
    "# Make sure the shapefile is uploaded in the shapes folder\n",
    "country_boundaries <- sf::st_read(file.path(paste0(shp.path, \"roi_kansas_us_epsg_4326.shp\")), \n",
    "                                  quiet = TRUE)\n",
    "\n",
    "# Check and reproject if necessary to match target CRS\n",
    "if(sf::st_crs(country_boundaries) != epsg) {\n",
    "  cat(\"Reprojecting country boundaries to\", epsg, \"\\n\")\n",
    "  country_boundaries <- country_boundaries %>%\n",
    "    sf::st_as_sf() %>% \n",
    "    sf::st_transform(crs = epsg)\n",
    "}\n",
    "\n",
    "# Load digital elevation model (DEM)\n",
    "# This file should contain elevation in meters\n",
    "# The DEM was downloaded using GEE script\n",
    "# Move this file to the raster folder\n",
    "elev <- list.files(raster.path, \n",
    "                   pattern = \"DEM_KANSAS.tif$\",  \n",
    "                   recursive = TRUE, \n",
    "                   full.names = TRUE)\n",
    "elev <- terra::rast(elev)\n",
    "#elev <- elev$dtm_elevation_250m\n",
    "\n",
    "# Reproject elevation to match target CRS if needed\n",
    "if(terra::crs(elev) != epsg) {\n",
    "  cat(\"Reprojecting DEM data to\", epsg, \"\\n\")\n",
    "  elev <- terra::project(elev, epsg, method = \"near\")\n",
    "}\n",
    "\n",
    "# Load available water capacity (AWC) data\n",
    "# AWC represents the soil's ability to hold water (0-200 cm depth)\n",
    "# Data source: Hengl and Gupta (https://zenodo.org/records/2629149)\n",
    "# This folder should be downloaded in the raster folder.\n",
    "awc <- terra::rast(paste0(raster.path, \"sol_available.water.capacity_usda.mm_m_250m_0..200cm_1950..2017_v0.1.tif\"))\n",
    "\n",
    "# Reproject AWC to match target CRS if needed (12 min)\n",
    "if(terra::crs(awc) != epsg) {\n",
    "  cat(\"Reprojecting AWC data to\", epsg, \"\\n\")\n",
    "  awc <- terra::project(awc, epsg, method = \"near\")\n",
    "}\n",
    "\n",
    "# Mask to the AOI geometry (sets outside cells to NA) (12 min)\n",
    "awc <- terra::crop(awc, country_boundaries)\n",
    "names(awc) <- \"awc\"\n",
    "plot(awc)\n",
    "\n",
    "# Resample AWC to match elevation grid\n",
    "awc <- terra::resample(awc, elev)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6 - PREPARE NEWHALL MODEL INPUTS"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 6 - PREPARE NEWHALL MODEL INPUTS ============================================\n",
    "\n",
    "# Load the climate variables created in previous steps\n",
    "newhall <- terra::rast(paste0(raster.path, \"climate_vars.tif\"))\n",
    "\n",
    "# Reproject climate data if needed\n",
    "if(terra::crs(newhall) != epsg) {\n",
    "  cat(\"Reprojecting climate data to\", epsg, \"\\n\")\n",
    "  newhall <- terra::project(newhall, epsg, method = \"near\")\n",
    "}\n",
    "\n",
    "# Resample to match elevation grid (ensures all layers align)\n",
    "newhall <- terra::resample(newhall, elev)\n",
    "\n",
    "# Add longitude and latitude as covariates\n",
    "# These are required for accurate Newhall calculations\n",
    "newhall$lonDD <- terra::init(newhall[[1]], 'x')  # Longitude in decimal degrees\n",
    "newhall$latDD <- terra::init(newhall[[1]], 'y')  # Latitude in decimal degrees\n",
    "names(newhall)\n",
    "\n",
    "# Mask climate data to elevation extent (remove NoData areas)\n",
    "newhall <- terra::mask(newhall, elev)\n",
    "\n",
    "# Add AWC and elevation to the stack\n",
    "newhall$awc <- awc\n",
    "newhall$elev <- elev\n",
    "\n",
    "# Visual check of some climate layers\n",
    "terra::plot(newhall[\"tJan\"], main = \"January Temperature (°C)\")\n",
    "terra::plot(newhall[\"pJan\"], main = \"January Precipitation (mm)\")\n",
    "\n",
    "# Display summary statistics\n",
    "print(summary(newhall))\n",
    "\n",
    "# Optional: Crop to administrative boundary\n",
    "# Uncomment the following line if you want to limit analysis to specific region\n",
    "# newhall <- terra::crop(newhall, country_boundaries, mask = TRUE)\n",
    "\n",
    "# Save complete Newhall input dataset\n",
    "terra::writeRaster(newhall, \n",
    "                   paste0(raster.path, \"climate_newhall_vars.tif\"), \n",
    "                   overwrite = TRUE)\n",
    "\n",
    "# Clean up\n",
    "rm(elev, awc)\n",
    "gc()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7 - COMPUTE NEWHALL SOIL-CLIMATE MODEL"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 7 - COMPUTE NEWHALL SOIL-CLIMATE MODEL ======================================\n",
    "\n",
    "# Optional: Aggregate raster to reduce computation time\n",
    "# Uncomment if processing large datasets\n",
    "# newhall <- terra::aggregate(newhall, fact = agg.factor)\n",
    "\n",
    "# Run Newhall batch processing\n",
    "# This calculates soil temperature and moisture regimes\n",
    "# cores = 4 enables parallel processing (adjust based on your CPU) - 6-7 hrs\n",
    "system.time({\n",
    "  newhall_results <- jNSMR::newhall_batch(newhall, cores = 4)\n",
    "})\n",
    "\n",
    "# Create mask from results (used to remove invalid areas)\n",
    "mask <- newhall_results$annualRainfall - newhall_results$annualRainfall + 1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8 - PLOT NEWHALL MODEL RESULTS"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 8 - PLOT NEWHALL MODEL RESULTS ==============================================\n",
    "\n",
    "# Water balance plots\n",
    "terra::plot(newhall_results$annualWaterBalance,\n",
    "            cex.main = 0.9, \n",
    "            main = \"Annual Water Balance (P-PET) [mm]\")\n",
    "\n",
    "terra::plot(newhall_results$summerWaterBalance,\n",
    "            cex.main = 0.9, \n",
    "            main = \"Summer Water Balance [mm]\")\n",
    "\n",
    "# Regime classifications\n",
    "terra::plot(newhall_results$temperatureRegime, \n",
    "            main = \"Temperature Regime\")\n",
    "\n",
    "terra::plot(newhall_results$moistureRegime, \n",
    "            main = \"Moisture Regime\")\n",
    "\n",
    "# Temporal metrics\n",
    "terra::plot(newhall_results$numCumulativeDaysDry, \n",
    "            col = grDevices::terrain.colors(20),\n",
    "            cex.main = 0.75, \n",
    "            main = \"Number of Cumulative Days Dry\")\n",
    "\n",
    "terra::plot(newhall_results$numCumulativeDaysDryOver5C, \n",
    "            col = grDevices::terrain.colors(20),\n",
    "            cex.main = 0.75, \n",
    "            main = \"Number of Cumulative Days Dry over 5°C\")\n",
    "\n",
    "terra::plot(newhall_results$numConsecutiveDaysMoistInSomePartsOver8C, \n",
    "            col = rev(grDevices::terrain.colors(20)),\n",
    "            cex.main = 0.75, \n",
    "            main = \"Number of Consecutive Days Moist\\nin some parts over 8°C\")\n",
    "\n",
    "terra::plot(newhall_results$dryDaysAfterSummerSolstice, \n",
    "            col = grDevices::terrain.colors(20),\n",
    "            cex.main = 0.75, \n",
    "            main = \"Number of Dry Days After Summer Solstice\")\n",
    "\n",
    "terra::plot(newhall_results$moistDaysAfterWinterSolstice,\n",
    "            col = rev(grDevices::terrain.colors(50)),\n",
    "            cex.main = 0.75, \n",
    "            main = \"Number of Moist Days After Winter Solstice\")\n",
    "\n",
    "# Regime subdivisions (if calculated)\n",
    "if(\"regimeSubdivision1\" %in% names(newhall_results)) {\n",
    "  terra::plot(newhall_results$regimeSubdivision1,\n",
    "              cex.main = 0.75, \n",
    "              main = \"Regime Subdivision 1\")\n",
    "}\n",
    "\n",
    "if(\"regimeSubdivision2\" %in% names(newhall_results)) {\n",
    "  terra::plot(newhall_results$regimeSubdivision2,\n",
    "              cex.main = 0.75, \n",
    "              main = \"Regime Subdivision 2\")\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9 - EXPORT RESULTS"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 9 - EXPORT RESULTS ==========================================================\n",
    "\n",
    "# Apply mask to results\n",
    "results <- terra::mask(newhall_results, mask)\n",
    "\n",
    "# Export main Newhall results\n",
    "terra::writeRaster(results, \n",
    "                   paste0(raster.path, \"newhall.tif\"), \n",
    "                   overwrite = TRUE)\n",
    "\n",
    "\n",
    "# Export climate input as integer (multiplied by 10 to preserve 1 decimal place)\n",
    "# This reduces file size while maintaining sufficient precision\n",
    "results_climate_int <- terra::app(newhall, function(x) { as.integer(x * 10) })\n",
    "names(results_climate_int) <- names(newhall)\n",
    "\n",
    "terra::writeRaster(results_climate_int,\n",
    "                   paste0(raster.path, \"climate_intx10.tif\"), \n",
    "                   overwrite = TRUE)\n",
    "\n",
    "# Export Newhall results as integer (multiplied by 10)\n",
    "results_newhall_int <- terra::app(results, function(x) { as.integer(x * 10) })\n",
    "names(results_newhall_int) <- names(results)\n",
    "\n",
    "terra::writeRaster(results_newhall_int,\n",
    "                   paste0(raster.path, \"newhall_intx10.tif\"), \n",
    "                   overwrite = TRUE)\n",
    "\n",
    "# Final cleanup\n",
    "rm(results_climate_int, results_newhall_int, mask, newhall, newhall_results)\n",
    "gc()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10 - MERGE TILED GEE EXPORTS"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## 10 - MERGE TILED GEE EXPORTS ================================================\n",
    "\n",
    "# Function to merge tiles for a given dataset (SPATIAL merge, not stacking)\n",
    "merge_gee_tiles <- function(base_name, input_path, output_path) {\n",
    "  \n",
    "  cat(\"\\n========================================\\n\")\n",
    "  cat(\"Processing:\", base_name, \"\\n\")\n",
    "  cat(\"========================================\\n\")\n",
    "  \n",
    "  # Find all tiles for this dataset\n",
    "  tile_files <- list.files(input_path,\n",
    "                           pattern = paste0(\"^\", base_name, \".*\\\\.tif$\"),\n",
    "                           full.names = TRUE)\n",
    "  \n",
    "  if(length(tile_files) == 0) {\n",
    "    cat(\"WARNING: No tiles found for\", base_name, \"\\n\")\n",
    "    return(NULL)\n",
    "  }\n",
    "  \n",
    "  cat(\"Found\", length(tile_files), \"tile(s)\\n\")\n",
    "  print(basename(tile_files))\n",
    "  \n",
    "  # If only one tile, just copy it\n",
    "  if(length(tile_files) == 1) {\n",
    "    cat(\"Single tile - copying directly\\n\")\n",
    "    merged <- terra::rast(tile_files[1])\n",
    "  } else {\n",
    "    # Load all tiles\n",
    "    cat(\"Loading tiles...\\n\")\n",
    "    tile_list <- lapply(tile_files, terra::rast)\n",
    "    \n",
    "    # SPATIAL MERGE using mosaic (stitches tiles together)\n",
    "    cat(\"Mosaicking tiles spatially...\\n\")\n",
    "    merged <- do.call(terra::mosaic, tile_list)\n",
    "  }\n",
    "  \n",
    "  # Create output filename\n",
    "  output_file <- paste0(output_path, base_name, \"_merged.tif\")\n",
    "  \n",
    "  # Export merged raster\n",
    "  cat(\"Exporting to:\", output_file, \"\\n\")\n",
    "  terra::writeRaster(merged, \n",
    "                     filename = output_file,\n",
    "                     overwrite = TRUE)\n",
    "  \n",
    "  cat(\"SUCCESS: Merged\", base_name, \"\\n\\n\")\n",
    "  return(merged)\n",
    "}\n",
    "\n",
    "# List of datasets to merge\n",
    "datasets_to_merge <- c(\n",
    "  \"ESA_WorldCover_2021_KANSAS\",\n",
    "  \"Cropland_Mask_KANSAS\",\n",
    "  \"HighRes_Covariates_100m_KANSAS\"\n",
    ")\n",
    "\n",
    "# Merge EACH dataset separately (not stacked together)\n",
    "for(dataset in datasets_to_merge) {\n",
    "  merge_gee_tiles(\n",
    "    base_name = dataset,\n",
    "    input_path = paste0(raster.path, \"/GEE_Exports/\"),\n",
    "    output_path = raster.path\n",
    "  )\n",
    "}\n",
    "\n",
    "\n",
    "################################################################################\n",
    "# END OF SCRIPT\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"
   ]
  }
 ]
}