{
 "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": [
    "# Modelling and mapping SOC with Quantile Regression Forest\n",
    "**Module 3 · Digital soil mapping** · SoilFER Training\n",
    "\n",
    "[Course page](https://training.yigini.net/modules/03-digital-soil-mapping/) · [Manual chapter](https://training.yigini.net/manual/digital-soil-mapping.html) · [Original script](https://github.com/SoilFER/SoilFER-Training-Resources/blob/main/02_scripts/module3/modelling_%26_mapping_v2.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(\"Boruta\", \"caret\", \"dplyr\", \"mapview\", \"ranger\", \"terra\", \"tidyverse\")\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": [
    "### Digital Soil Mapping: Modelling & Mapping"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "#\n",
    "# Digital Soil Mapping: Modelling & Mapping\n",
    "# SoilFER Training - Module 3\n",
    "#\n",
    "# Structure:\n",
    "#   SESSION 1: Data preparation, feature selection, \n",
    "#              modelling and accuracy assessment.\n",
    "#   SESSION 2: Prediction of conditional mean and standard deviation of SOC\n",
    "#\n",
    "\n",
    "rm(list = ls())\n",
    "gc()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SESSION 1"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# SESSION 1 ====================================================================\n",
    "\n",
    "# 1. Setup and packages --------------------------------------------------------\n",
    "\n",
    "library(tidyverse)\n",
    "library(caret)\n",
    "library(terra)\n",
    "library(Boruta)\n",
    "library(ranger)\n",
    "library(mapview)\n",
    "\n",
    "\n",
    "Sys.setenv(PROJ_LIB = \"\")  \n",
    "Sys.setenv(PROJ_DATA = system.file(\"proj\", package = \"terra\"))\n",
    "\n",
    "# Set working directory to script location (works in RStudio)\n",
    "setwd(file.path(Sys.getenv(\"SOILFER_ROOT\"), \"02_scripts/module3\"))  # folder of this script (was rstudioapi::getActiveDocumentContext, RStudio only)\n",
    "setwd(\"../../\")\n",
    "\n",
    "# Create output directories\n",
    "for (d in c(\"03_outputs/module3/models/\", \n",
    "            \"03_outputs/module3/validation/\", \n",
    "            \"03_outputs/module3/tiles/\",\n",
    "            \"03_outputs/module3/maps/\",\n",
    "            \"03_outputs/module3/maps/aoa/\",\n",
    "            \"terra_tmp\")) {\n",
    "  if (!dir.exists(d)) dir.create(d, recursive = TRUE)\n",
    "}\n",
    "\n",
    "# Terra settings: show progress bar, use 60% of RAM, use local temp folder\n",
    "terraOptions(progress = 1, memfrac = 0.6, tempdir = file.path(getwd(), \"terra_tmp\"))\n",
    "\n",
    "\n",
    "# 2. Load covariates -----------------------------------------------------------\n",
    "\n",
    "covs <- rast(\"01_data/module1/training_data/Environmental_Covariates_250m_KANSAS.tif\")  # EDIT THIS: your covariate stack\n",
    "\n",
    "cov_names <- names(covs)\n",
    "\n",
    "# Inspect the covariate raster\n",
    "covs[[1]]  # Shows first covariate resolution, extent, CRS\n",
    "nlyr(covs) # Number of layers\n",
    "plot(covs[[1]])\n",
    "\n",
    "# 3. Load and transform soil data ----------------------------------------------\n",
    "\n",
    "dat <- read_csv(\"03_outputs/module1/KSSL_DSM_0-30.csv\") \n",
    "\n",
    "# Inspect the data\n",
    "dat\n",
    "summary(dat)\n",
    "\n",
    "## 3.1 Create bulk density with a pedotransfer function (Saxton)----------------\n",
    "\n",
    "dat <- dat %>%\n",
    "  mutate(BD = 1.35 + 0.0045 * Sand + 0.0035 * Clay - 0.06 * 1.72 * SOC)\n",
    "\n",
    "# Bonus track: Estimate Available Water Capacity (Rawl and Saxto 1982)\n",
    "\n",
    "# dat <- dat %>% \n",
    "#   mutate(\n",
    "#     # 1. Convert percentage inputs to decimal fractions for model constants\n",
    "#     S = Sand / 100,\n",
    "#     C = Clay / 100,\n",
    "#     # Convert SOC to Organic Matter (OM) using the Van Bemmelen factor and scale to decimal\n",
    "#     OM = (SOC * 1.724) / 100,\n",
    "#     \n",
    "#     # 2. Estimate Wilting Point (WP) at 1500 kPa\n",
    "#     # First solution for 1500 kPa moisture\n",
    "#     theta_1500t = -0.024 * S + 0.487 * C + 0.006 * OM + \n",
    "#       0.005 * (S * OM) - 0.013 * (C * OM) + \n",
    "#       0.068 * (S * C) + 0.031,\n",
    "#     \n",
    "#     # Final WP adjustment using the lack-of-fit equation\n",
    "#     WP = theta_1500t + (0.14 * theta_1500t - 0.02),\n",
    "#     \n",
    "#     # 3. Estimate Field Capacity (FC) at 33 kPa\n",
    "#     # First solution for 33 kPa moisture\n",
    "#     theta_33t = -0.251 * S + 0.195 * C + 0.011 * OM + \n",
    "#       0.006 * (S * OM) - 0.027 * (C * OM) - \n",
    "#       0.452 * (S * C) + 0.299,\n",
    "#     \n",
    "#     # Final FC adjustment using the lack-of-fit equation\n",
    "#     FC = theta_33t + (1.283 * (theta_33t^2) - 0.374 * theta_33t - 0.015),\n",
    "#     \n",
    "#     # 4. Calculate Available Water Capacity (AWC)\n",
    "#     # Result is in volumetric water content (%v) as a decimal fraction\n",
    "#     AWC = FC - WP\n",
    "#   ) %>%\n",
    "#   # Clean up auxiliary variables\n",
    "#   select(-S, -C, -theta_1500t, -theta_33t, OM)\n",
    "\n",
    "\n",
    "# 4. Merge soil data and covariate values --------------------------------------\n",
    "# Covert the dataframe into a spatial object (points)\n",
    "dat_pts <- vect(dat, geom = c(\"lon\", \"lat\"), crs = \"epsg:4326\")\n",
    "dat_pts <- terra::project(dat_pts, covs)\n",
    "mapview(dat_pts, cex=1.5) #+ mapview(covs[[1]])\n",
    "\n",
    "# Extract covariates at sample locations \n",
    "extracted_covs <- terra::extract(x = covs, y = dat_pts, xy = TRUE, ID = FALSE)\n",
    "summary(extracted_covs)\n",
    "dat <- as.data.frame(dat_pts)\n",
    "dat_cov <- bind_cols(dat, extracted_covs)\n",
    "\n",
    "# 5. Feature selection with Boruta ---------------------------------------------\n",
    "# Define the target soil property\n",
    "target <- \"SOC\"  # EDIT THIS: your continuous target variable\n",
    "\n",
    "# Prepare training data (target + covariates, complete cases only)\n",
    "d <- dat_cov %>%\n",
    "  dplyr::select(all_of(target), all_of(cov_names)) %>%\n",
    "  na.omit() %>% \n",
    "  as.data.frame()\n",
    "# Observe d\n",
    "View(d)\n",
    "\n",
    "# Run Boruta feature selection\n",
    "set.seed(1)\n",
    "boruta_result <- Boruta(\n",
    "  y = d[,target],\n",
    "  x = d[,cov_names],\n",
    "  maxRuns = 25, # this should be >100\n",
    "  doTrace = 1  # Show progress\n",
    ")\n",
    "\n",
    "# Save Boruta plot\n",
    "png(paste0(\"03_outputs/module3/figures/boruta_\",target,\".png\"), \n",
    "    width = 15, height = 20, units = \"cm\", res = 150)\n",
    "par(las = 1, mar = c(4, 10, 4, 2) + 0.1)\n",
    "plot(boruta_result, horizontal = TRUE,las=1, \n",
    "          ylab = \"\", xlab = \"Importance\", cex.axis = 0.6)\n",
    "dev.off()\n",
    "\n",
    "# Observe the figure\n",
    "par(las = 1, mar = c(4, 10, 4, 2) + 0.1)\n",
    "plot(boruta_result, horizontal = TRUE,las=1, \n",
    "     ylab = \"\", xlab = \"Importance\", cex.axis = 0.6)\n",
    "\n",
    "# Get selected features (confirmed (green) + tentative (yellow))\n",
    "selected_features <- getSelectedAttributes(boruta_result, withTentative = TRUE)\n",
    "\n",
    "\n",
    "# 6. Train and validate the Quantile Regression Forest (qrf) model -------------\n",
    "\n",
    "# Set up Cross-validation process\n",
    "cv_control <- trainControl(\n",
    "  method = \"repeatedcv\",\n",
    "  number = 5,   # ideal number 10-20\n",
    "  repeats = 5,   # idel number 10-20\n",
    "  savePredictions = TRUE\n",
    ")\n",
    "\n",
    "# Set up hyperparameter grid for ranger\n",
    "mtry_base <- round(length(selected_features) / 3)\n",
    "tune_grid <- expand.grid(\n",
    "  mtry = c(mtry_base - round(mtry_base/2), \n",
    "           mtry_base, \n",
    "           mtry_base + round(mtry_base/2)\n",
    "           ),\n",
    "  min.node.size = 5,\n",
    "  splitrule = c(\"variance\", \"extratrees\")\n",
    ")\n",
    "\n",
    "# Train the model\n",
    "# quantreg=TRUE enables quantile predictions (for uncertainty maps)\n",
    "# num.threads uses multiple CPU cores for faster training\n",
    "model <- caret::train(\n",
    "  y = d[,target],\n",
    "  x = d[,selected_features],\n",
    "  method = \"ranger\",\n",
    "  quantreg = TRUE,\n",
    "  importance = \"permutation\",\n",
    "  trControl = cv_control,\n",
    "  tuneGrid = tune_grid,\n",
    "  num.threads = max(1, parallel::detectCores() - 1)  # EDIT THIS: number of CPU threads\n",
    ")\n",
    "\n",
    "model\n",
    "\n",
    "model$bestTune\n",
    "\n",
    "# Variance importance\n",
    "varImp(model)\n",
    "\n",
    "# Save variable importance plot\n",
    "png(paste0(\"03_outputs/module3/figures/varImp_\",target,\".png\"),\n",
    "    width = 15, height = 15, units = \"cm\", res = 150)\n",
    "plot(varImp(model), main = paste(\"Variable Importance -\", target))\n",
    "dev.off()\n",
    "\n",
    "# View the variable importance\n",
    "plot(varImp(model), main = paste(\"Variable Importance -\", target))\n",
    "\n",
    "\n",
    "# 7. Accuracy assessment -------------------------------------------------------\n",
    "\n",
    "# Extract CV predictions for the best hyperparameter combination\n",
    "cv_preds <- model$pred %>%\n",
    "  filter(mtry == model$bestTune$mtry,\n",
    "         splitrule == model$bestTune$splitrule,\n",
    "         min.node.size == model$bestTune$min.node.size)\n",
    "\n",
    "obs_vals <- cv_preds$obs\n",
    "pred_vals <- cv_preds$pred\n",
    "\n",
    "# Calculate validation metrics\n",
    "load(\"02_scripts/module3/eval.RData\")\n",
    "accuracy <- eval(pred_vals, obs_vals)[, 1:6]\n",
    "\n",
    "accuracy\n",
    "\n",
    "# Scatterplot of observed vs predicted\n",
    "residuals <- tibble(Observed = obs_vals, Predicted = pred_vals)\n",
    "g_scatter <- ggplot(residuals, aes(x = Observed, y = Predicted)) +\n",
    "  geom_point(alpha = 0.3) +\n",
    "  geom_abline(slope = 1, intercept = 0, color = \"red\", linewidth = 1) +\n",
    "  ylim(c(min(residuals$Observed), max(residuals$Observed))) + theme(aspect.ratio=1)+\n",
    "  coord_fixed() +\n",
    "  labs(title = paste(\"Observed vs Predicted -\", target)) +\n",
    "  theme_minimal()\n",
    "\n",
    "g_scatter\n",
    "\n",
    "ggsave(g_scatter, \n",
    "       filename = paste0(\"03_outputs/module3/figures/scatterplot_\",target,\".png\"),\n",
    "       width = 12, height = 12, units = \"cm\")\n",
    "\n",
    "\n",
    "write_csv(accuracy, paste0(\"03_outputs/module3/validation/\", target,\"_accuracy.csv\"))\n",
    "\n",
    "\n",
    "# Save the qrf model\n",
    "saveRDS(model, paste0(\"03_outputs/module3/models/ranger_model_\",target,\".rds\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SESSION 2"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# SESSION 2 ====================================================================\n",
    "\n",
    "# Set working directory to script location (works in RStudio)\n",
    "setwd(file.path(Sys.getenv(\"SOILFER_ROOT\"), \"02_scripts/module3\"))  # folder of this script (was rstudioapi::getActiveDocumentContext, RStudio only)\n",
    "setwd(\"../../\")\n",
    "\n",
    "# Load required libraries and objects\n",
    "\n",
    "library(tidyverse)\n",
    "library(caret)\n",
    "library(terra)\n",
    "library(Boruta)\n",
    "library(ranger)\n",
    "library(mapview)\n",
    "Sys.setenv(PROJ_LIB = \"\")  \n",
    "Sys.setenv(PROJ_DATA = system.file(\"proj\", package = \"terra\"))\n",
    "\n",
    "covs <- rast(\"01_data/module1/training_data/Environmental_Covariates_250m_KANSAS.tif\")\n",
    "cov_names <- names(covs)\n",
    "dat <- read_csv(\"03_outputs/module1/KSSL_DSM_0-30.csv\") \n",
    "target <- \"SOC\"\n",
    "model <- readRDS(paste0(\"03_outputs/module3/models/ranger_model_\",target,\".rds\"))\n",
    "\n",
    "\n",
    "# 8. Spatial prediction (tiled) ------------------------------------------------\n",
    "\n",
    "## 8.1 Create tile grid for memory-efficient prediction ------------------------\n",
    "\n",
    "r <-covs[[1]]\n",
    "t <- rast(nrows = 5, ncols = 10, extent = ext(r), crs = crs(r))\n",
    "tile <- makeTiles(r, t,overwrite=TRUE,filename=\"03_outputs/module3/tiles/tiles.tif\")\n",
    "\n",
    "# vew the first tile\n",
    "rast(tile[1]) %>% plot\n",
    "\n",
    "## 8.2 Prepare model for prediction (single-threaded per tile) -----------------\n",
    "ranger_model <- model$finalModel\n",
    "\n",
    "# Prediction function for quantile RF (returns transposed predictions)\n",
    "pfun <- function(...) {\n",
    "  predict(...)$predictions |> t()\n",
    "}\n",
    "\n",
    "# GDAL write options: compress output, use tiled format for faster reading\n",
    "gdal_opts <- c(\"COMPRESS=LZW\")\n",
    "\n",
    "\n",
    "## 8.3 Prediction for one tile -------------------------------------------------\n",
    "\n",
    "t1 <- rast(tile[1])\n",
    "\n",
    "# crop the selected covariates with the tile 1\n",
    "covs_t1 <- crop(covs, t1)\n",
    "\n",
    "plot(covs_t1[[1:9]])\n",
    "\n",
    "# Prediction of the conditional mean\n",
    "pred_mean_t1 <- terra::interpolate(covs_t1, \n",
    "                              model = ranger_model, \n",
    "                              fun=pfun, \n",
    "                              na.rm=TRUE,\n",
    "                              type = \"quantiles\",\n",
    "                              what=mean)\n",
    "\n",
    "# Prediction of the conditional standard deviation\n",
    "pred_sd_t1 <- terra::interpolate(covs_t1, \n",
    "                                 model = ranger_model, \n",
    "                                 fun=pfun, \n",
    "                                 na.rm=TRUE,\n",
    "                                 type = \"quantiles\",\n",
    "                                 what=sd)\n",
    "\n",
    "# Compute the coeficient of variation\n",
    "pred_cv_t1 <- pred_sd_t1/pred_mean_t1\n",
    "\n",
    "\n",
    "# Explore the results at one tile\n",
    "mapview(pred_mean_t1) \n",
    "\n",
    "\n",
    "\n",
    "## 8.4 Prediction at all tiles (it may take time) ------------------------------\n",
    "\n",
    "# Loop through tiles and predict mean + SD\n",
    "for (j in seq_along(tile)) {\n",
    "  print(paste(\"Processing tile\", j, \"of\", length(tile)))\n",
    "\n",
    "  # Crop covariates to this tile\n",
    "  t <- rast(tile[j])\n",
    "  covs_tile <- crop(covs, t)\n",
    "\n",
    "  # Predict MEAN\n",
    "  mean_file <- paste0(\"03_outputs/module3/tiles/\", target, \"_mean_\", j, \".tif\")\n",
    "  terra::interpolate(\n",
    "    covs_tile,\n",
    "    model = ranger_model,\n",
    "    fun = pfun,\n",
    "    na.rm = TRUE,\n",
    "    type = \"quantiles\",\n",
    "    what = mean,\n",
    "    filename = mean_file,\n",
    "    overwrite = TRUE,\n",
    "    wopt = list(datatype = \"FLT4S\", gdal = gdal_opts)\n",
    "  )\n",
    "\n",
    "  # Predict SD (uncertainty)- this will take some time to load\n",
    "  sd_file <- paste0(\"03_outputs/module3/tiles/\", target, \"_sd_\", j, \".tif\")\n",
    "  terra::interpolate(\n",
    "    covs_tile,\n",
    "    model = ranger_model,\n",
    "    fun = pfun,\n",
    "    na.rm = TRUE,\n",
    "    type = \"quantiles\",\n",
    "    what = sd,\n",
    "    filename = sd_file,\n",
    "    overwrite = TRUE,\n",
    "    wopt = list(datatype = \"FLT4S\", gdal = gdal_opts)\n",
    "  )\n",
    "}\n",
    "\n",
    "# 9. Mosaic tiles into final maps ----------------------------------------------\n",
    "\n",
    "## 9.1 Mosaic all mean tiles ---------------------------------------------------\n",
    "# Find all mean tiles and mosaic them\n",
    "mean_tiles <- list.files(\"03_outputs/module3/tiles/\",\n",
    "                         pattern = paste0(\"^\", target, \"_mean_.*\\\\.tif$\"),\n",
    "                         full.names = TRUE)\n",
    "\n",
    "# Create a SpatRasterCollection and mosaic\n",
    "mean_rasters <- list()\n",
    "for (i in 1:length(mean_tiles)) {\n",
    "  mean_rasters[[i]] <- rast(mean_tiles[i])\n",
    "}\n",
    "pred_mean <- mosaic(sprc(mean_rasters), fun=\"first\")\n",
    "names(pred_mean) <- paste(target, \"_mean\")\n",
    "\n",
    "# Save final mean map\n",
    "writeRaster(pred_mean, paste0(\"03_outputs/module3/maps/mean_\", target, \".tif\"),\n",
    "            overwrite = TRUE, wopt = list(datatype = \"FLT4S\", gdal = gdal_opts))\n",
    "\n",
    "## 9.2 Mosaic all sd tiles -----------------------------------------------------\n",
    "# Same for SD tiles\n",
    "sd_tiles <- list.files(\"03_outputs/module3/tiles/\",\n",
    "                       pattern = paste0(\"^\", target, \"_sd_.*\\\\.tif$\"),\n",
    "                       full.names = TRUE)\n",
    "\n",
    "sd_rasters <- list()\n",
    "for (i in 1:length(sd_tiles)) {\n",
    "  sd_rasters[[i]] <- rast(sd_tiles[i])\n",
    "}\n",
    "pred_sd <- mosaic(sprc(sd_rasters), fun=\"first\")\n",
    "names(pred_sd) <- paste(target, \"_sd\")\n",
    "\n",
    "writeRaster(pred_sd, paste0(\"03_outputs/module3/maps/sd_\", target, \".tif\"),\n",
    "            overwrite = TRUE, wopt = list(datatype = \"FLT4S\", gdal = gdal_opts))\n",
    "\n",
    "# Display maps\n",
    "plot(pred_mean, main = paste(\"Mean\", target), col = hcl.colors(100, \"Viridis\"))\n",
    "plot(pred_sd, main = paste(\"SD\", target), col = hcl.colors(100, \"Viridis\"))\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"
   ]
  }
 ]
}