{
 "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": [
    "# Part 4 · Additional exercises\n",
    "**Module 4 · Soil spectroscopy for digital soil mapping** · SoilFER Training\n",
    "\n",
    "[Course page](https://training.yigini.net/modules/04-soil-spectroscopy/) · [Manual chapter](https://training.yigini.net/manual/soil-spectroscopy-for-digital-soil-mapping.html) · [Original script](https://github.com/SoilFER/SoilFER-Training-Resources/blob/main/02_scripts/module4/004-Soil%20Spectroscopy%20for%20Digital%20Soil%20Mapping.R)\n",
    "\n",
    "**Before you start**\n",
    "1. Check that the runtime is **R**: *Runtime → Change runtime type → R* (this notebook should open in R automatically).\n",
    "2. Run the **Setup** cell below once per session (≈1–3 minutes). It downloads the training project, installs the R packages and downloads the course rasters and MIR data (≈1.2 GB) and sets the working folder.\n",
    "3. Then run the cells in order with **Shift + Enter**.\n",
    "\n",
    "> Colab resets when you close it or after ~90 minutes without activity. Save results you want to keep with *Files → Download* (left sidebar), or re-run the Setup cell after a reset.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# ==== SoilFER · Colab setup (run first, once per session) ==================\n",
    "options(repos = c(CRAN = \"https://cloud.r-project.org\"), timeout = 3600)\n",
    "# Works in Google Colab and in any other Jupyter (JupyterHub, JupyterLab on your computer)\n",
    "on_colab <- nzchar(Sys.getenv(\"COLAB_RELEASE_TAG\")) || dir.exists(\"/content/sample_data\")\n",
    "root <- if (on_colab) \"/content/SoilFER-Training-Resources\" else path.expand(\"~/SoilFER-Training-Resources\")\n",
    "Sys.setenv(SOILFER_ROOT = root)\n",
    "\n",
    "# 1. Training project (scripts, small data, outputs, assignments)\n",
    "if (!dir.exists(root))\n",
    "  system(paste(\"git clone --depth 1 https://github.com/SoilFER/SoilFER-Training-Resources\", root))\n",
    "\n",
    "# 2. R packages (in Colab the R runtime installs ready-made binaries, so this is fast)\n",
    "pkgs <- c(\"aqp\", \"dplyr\", \"ggplot2\", \"matrixStats\", \"prospectr\", \"ranger\", \"readxl\", \"terra\")\n",
    "need <- setdiff(pkgs, rownames(installed.packages()))\n",
    "if (length(need)) install.packages(need)\n",
    "\n",
    "# 3. Course rasters + MIR spectra (Google Drive folder of the SoilFER training, ≈1.2 GB)\n",
    "td <- file.path(root, \"01_data/module1/training_data\")\n",
    "if (!file.exists(file.path(td, \"MIR_KANSAS_data.xlsx\"))) {\n",
    "  system(\"python3 -m pip -q install gdown\")\n",
    "  drv <- file.path(dirname(root), \"soilfer_drive\")\n",
    "  system(paste(\"python3 -m gdown --folder --quiet https://drive.google.com/drive/folders/1K7tq9zX5HsqbqWcNoT27WtfPtehcKBCu -O\", shQuote(drv)))\n",
    "  f <- list.files(drv, recursive = TRUE, full.names = TRUE)\n",
    "  file.copy(f, td, overwrite = FALSE)\n",
    "}\n",
    "\n",
    "setwd(root)\n",
    "cat(\"Ready. Working folder:\", getwd(), \"\\n\")\n",
    "missing <- setdiff(pkgs, rownames(installed.packages()))\n",
    "if (length(missing)) message(\"Not installed: \", paste(missing, collapse = \", \"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fresh runtime? Re-create the objects from the previous session(s)\n",
    "This session uses objects created earlier. If you just opened this notebook (or Colab was reset), run this cell first. It takes a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "for (s in c(\"02_scripts/module4/001-Soil Spectroscopy for Digital Soil Mapping.R\", \"02_scripts/module4/002-Soil Spectroscopy for Digital Soil Mapping.R\", \"02_scripts/module4/003-Soil Spectroscopy for Digital Soil Mapping.R\")) { message(\"Running \", s); source(file.path(root, s), echo = FALSE) }\n",
    "setwd(root)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CODE FOR SOILFER MODULE 4:"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "## CODE FOR SOILFER MODULE 4: \n",
    "## SOIL SECTROSCOPY FOR DIGITAL SOIL MAPPING\n",
    "## Alex Wadoux and Leonardo Ramirez-Lopez\n",
    "\n",
    "## Additional excerisies part\n",
    "\n",
    "## Preceeding code: 003-Soil Spectroscopy for Digital Soil Mapping.R\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: PCA on processed spectra"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Step 1: PCA on processed spectra ##############################\n",
    "################################################################################\n",
    "\n",
    "library(prospectr)\n",
    "library(ggplot2)\n",
    "\n",
    "pca_out <- prcomp(dat$spc_processed, center = TRUE, scale. = FALSE)\n",
    "scores_raw <- pca_out$x[, 1:10]\n",
    "scores_sd <- apply(scores_raw, 2, sd)\n",
    "scores <- sweep(scores_raw, 2, scores_sd, \"/\")\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: KDE of the full population"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Step 2: KDE of the full population ############################\n",
    "################################################################################\n",
    "\n",
    "n_grid <- 512\n",
    "kde_full <- vector(\"list\", ncol(scores))\n",
    "\n",
    "for (j in seq_len(ncol(scores))) {\n",
    "  kde_full[[j]] <- density(\n",
    "    scores[, j],\n",
    "    bw = \"nrd0\", n = n_grid,\n",
    "    from = min(scores[, j]), to = max(scores[, j])\n",
    "  )\n",
    "}\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: MSD across set sizes"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Step 3: MSD across set sizes ##################################\n",
    "################################################################################\n",
    "\n",
    "# Range of calibration set sizes to evaluate\n",
    "set_sizes <- seq(10, 200, by = 10)\n",
    "\n",
    "# k-means is stochastic → repeat each set size several times\n",
    "repetitions <- 10\n",
    "\n",
    "# Storage: MSD matrix (rows = sizes, cols = reps), PC1 KDE per size, mean sample count\n",
    "msd_matrix <- matrix(NA, nrow = length(set_sizes), ncol = repetitions)\n",
    "kde_subsets_pc1 <- vector(\"list\", length(set_sizes))\n",
    "n_samples_vec <- numeric(length(set_sizes))\n",
    "\n",
    "# Loop over candidate calibration set sizes\n",
    "for (i in seq_along(set_sizes)) {\n",
    "  n_samp_rep <- numeric(repetitions)\n",
    "  \n",
    "  # Repeat to capture k-means initialisation variability\n",
    "  for (r in seq_len(repetitions)) {\n",
    "    set.seed(r * i)\n",
    "    \n",
    "    # Select k profiles by k-means in PC space\n",
    "    kms_i <- prospectr::naes(scores, k = set_sizes[i], iter.max = 1000)\n",
    "    \n",
    "    # Profile-complete: include all horizons of selected profiles\n",
    "    idx_i <- which(dat$ProfID %in% unique(dat$ProfID[kms_i$model]))\n",
    "    n_samp_rep[r] <- length(idx_i)\n",
    "    \n",
    "    # MSD per PC, then averaged across the 10 PCs\n",
    "    msd_pcs <- numeric(ncol(scores))\n",
    "    for (j in seq_len(ncol(scores))) {\n",
    "      # KDE on the same grid/bandwidth as the full-population KDE\n",
    "      kde_j <- density(\n",
    "        scores[idx_i, j],\n",
    "        bw = kde_full[[j]]$bw, n = n_grid,\n",
    "        from = min(scores[, j]), to = max(scores[, j])\n",
    "      )\n",
    "      # Mean squared distance between the two density estimates\n",
    "      msd_pcs[j] <- mean((kde_j$y - kde_full[[j]]$y)^2)\n",
    "    }\n",
    "    msd_matrix[i, r] <- mean(msd_pcs)\n",
    "    \n",
    "    # Keep the first PC1 KDE per set size for the convergence plot\n",
    "    if (r == 1) {\n",
    "      kde_subsets_pc1[[i]] <- density(\n",
    "        scores[idx_i, 1],\n",
    "        bw = kde_full[[1]]$bw, n = n_grid,\n",
    "        from = min(scores[, 1]), to = max(scores[, 1])\n",
    "      )\n",
    "    }\n",
    "  }\n",
    "  # Average number of horizons selected (used for the dual axis later)\n",
    "  n_samples_vec[i] <- round(mean(n_samp_rep))\n",
    "}\n",
    "\n",
    "# Summarise MSD across repetitions\n",
    "msd_mean <- rowMeans(msd_matrix)\n",
    "msd_sd <- apply(msd_matrix, 1, sd)\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create the DRS calibration subsets"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Create the DRS calibration subsets ############################\n",
    "################################################################################\n",
    "\n",
    "n_profiles <- c(20, 50, 85, 140, 200, 300)\n",
    "cal_indices <- vector(\"list\", length(n_profiles))\n",
    "names(cal_indices) <- as.character(n_profiles)\n",
    "\n",
    "for (i in seq_along(n_profiles)) {\n",
    "  set.seed(201909)\n",
    "  kms_i <- prospectr::naes(\n",
    "    dat$spc_processed,\n",
    "    k = n_profiles[i],\n",
    "    pc = 10,\n",
    "    iter.max = 1000\n",
    "  )\n",
    "  cal_indices[[i]] <- which(dat$ProfID %in% unique(dat$ProfID[kms_i$model]))\n",
    "}\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit DRS models for each size"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Fit DRS models for each size ##################################\n",
    "################################################################################\n",
    "\n",
    "drs_preds <- vector(\"list\", length(n_profiles))\n",
    "names(drs_preds) <- as.character(n_profiles)\n",
    "\n",
    "for (i in seq_along(n_profiles)) {\n",
    "  cal_idx <- cal_indices[[i]]\n",
    "  pred_idx <- setdiff(seq_len(nrow(dat)), cal_idx)\n",
    "  \n",
    "  drs_model_i <- ranger(\n",
    "    x = dat$spc_processed[cal_idx, ], y = dat$SOC[cal_idx],\n",
    "    quantreg = TRUE, num.trees = 3500,\n",
    "    sample.fraction = 1, replace = TRUE,\n",
    "    splitrule = \"maxstat\", min.node.size = 10, seed = 201909\n",
    "  )\n",
    "  \n",
    "  preds_i <- predict(drs_model_i,\n",
    "                     data = dat$spc_processed[pred_idx, ],\n",
    "                     type = \"quantiles\",\n",
    "                     what = function(x) sample(x, 100, replace = TRUE)\n",
    "  )$predictions\n",
    "  \n",
    "  drs_preds[[i]] <- data.frame(\n",
    "    row_idx = pred_idx,\n",
    "    SOC_predRF = rowMeans(preds_i),\n",
    "    SOC_sdRF = sqrt(rowVars(preds_i))\n",
    "  )\n",
    "}\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assemble the augmented datasets"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Assemble the augmented datasets ###############################\n",
    "################################################################################\n",
    "\n",
    "aug_datasets <- vector(\"list\", length(n_profiles))\n",
    "names(aug_datasets) <- as.character(n_profiles)\n",
    "\n",
    "my_col_names <- c(\"smp_id\", \"ProfID\", \"Long_Site.x\", \"Lat_Site.x\",\n",
    "                  \"Top_depth_cm.x\", \"Bottom_depth_cm.x\")\n",
    "\n",
    "for (i in seq_along(n_profiles)) {\n",
    "  cal_idx <- cal_indices[[i]]\n",
    "  pred_idx <- drs_preds[[i]]$row_idx\n",
    "  \n",
    "  aug_i <- as.data.frame(dat)[, my_col_names]\n",
    "  aug_i$SOC <- NA\n",
    "  aug_i$SOC_sd <- NA\n",
    "  \n",
    "  aug_i$SOC[cal_idx] <- dat$SOC[cal_idx]\n",
    "  aug_i$SOC_sd[cal_idx] <- 0.15\n",
    "  \n",
    "  aug_i$SOC[pred_idx] <- drs_preds[[i]]$SOC_predRF\n",
    "  aug_i$SOC_sd[pred_idx] <- drs_preds[[i]]$SOC_sdRF\n",
    "  \n",
    "  aug_datasets[[i]] <- aug_i\n",
    "}\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hold out a fixed validation set"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Hold out a fixed validation set ###############################\n",
    "################################################################################\n",
    "\n",
    "soil_df_ref$ProfID <- soc_ref_0_30_xy$SOC_0_30[!is.na(soc_ref_0_30_xy$SOC_0_30)]\n",
    "\n",
    "set.seed(202409)\n",
    "val_profiles <- sample(unique(soil_df_ref$ProfID), size = 50)\n",
    "val_data <- soil_df_ref[soil_df_ref$ProfID %in% val_profiles, ]\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Harmonise to 0-30 cm"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Harmonise to 0-30 cm ##########################################\n",
    "################################################################################\n",
    "\n",
    "depth_weighted_mean_0_30 <- function(p) {\n",
    "  h <- horizons(p)\n",
    "  w <- pmax(0, pmin(h$Bottom_depth_cm.x, 30) - pmax(h$Top_depth_cm.x, 0))\n",
    "  keep <- w > 0 & !is.na(h$SOC)\n",
    "  if (sum(w[keep]) < 30) return(NA)\n",
    "  sum(w[keep] * h$SOC[keep]) / sum(w[keep])\n",
    "}\n",
    "\n",
    "aug_0_30 <- vector(\"list\", length(n_profiles))\n",
    "names(aug_0_30) <- as.character(n_profiles)\n",
    "\n",
    "for (i in seq_along(n_profiles)) {\n",
    "  cal_idx <- cal_indices[[i]]\n",
    "  aug_i <- aug_datasets[[i]]\n",
    "  depths(aug_i) <- ProfID ~ Top_depth_cm.x + Bottom_depth_cm.x\n",
    "  \n",
    "  soc_mean_i <- profileApply(aug_i, depth_weighted_mean_0_30)\n",
    "  \n",
    "  soc_sd_i <- horizons(aug_i) |>\n",
    "    as.data.frame() |>\n",
    "    group_by(ProfID) |>\n",
    "    summarise(SOC_sd_0_30 = weighted_sd_0_30(\n",
    "      Top_depth_cm.x, Bottom_depth_cm.x, SOC_sd),\n",
    "      .groups = \"drop\")\n",
    "  \n",
    "  coords_i <- as.data.frame(aug_i)[, c(\"ProfID\", \"Lat_Site.x\", \"Long_Site.x\")]\n",
    "  coords_i <- coords_i[!duplicated(coords_i$ProfID), ]\n",
    "  \n",
    "  soc_0_30_i <- soc_sd_i |>\n",
    "    mutate(SOC_0_30 = soc_mean_i) |>\n",
    "    select(ProfID, SOC_0_30, SOC_sd_0_30) |>\n",
    "    merge(coords_i, by = \"ProfID\", all.x = TRUE)\n",
    "  \n",
    "  # Lab uncertainty fixed at profile level\n",
    "  cal_profiles <- unique(dat$ProfID[cal_idx])\n",
    "  soc_0_30_i$SOC_sd_0_30[soc_0_30_i$ProfID %in% cal_profiles] <- 0.15\n",
    "  \n",
    "  aug_0_30[[i]] <- soc_0_30_i[!soc_0_30_i$ProfID %in% val_profiles, ]\n",
    "}\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit DSM models for each scenario"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Fit DSM models for each scenario ##############################\n",
    "################################################################################\n",
    "\n",
    "dsm_models <- vector(\"list\", length(n_profiles))\n",
    "soc_maps <- vector(\"list\", length(n_profiles))\n",
    "names(dsm_models) <- as.character(n_profiles)\n",
    "names(soc_maps) <- as.character(n_profiles)\n",
    "\n",
    "my_formula <- as.formula(paste(\"SOC_0_30 ~\", paste(cov_names, collapse = \" + \")))\n",
    "nms <- c(\"x\", \"y\", \"SOC_0_30\", \"SOC_sd_0_30\", cov_names)\n",
    "\n",
    "for (i in seq_along(n_profiles)) {\n",
    "  pts_i <- vect(aug_0_30[[i]],\n",
    "                geom = c(\"Long_Site.x\", \"Lat_Site.x\"),\n",
    "                crs = \"EPSG:4326\")\n",
    "  pts_i <- project(pts_i, covs)\n",
    "  pts_i <- cbind(pts_i, terra::extract(covs, pts_i)[, -1])\n",
    "  \n",
    "  df_i <- as.data.frame(pts_i, geom = \"XY\")[, nms]\n",
    "  df_i <- df_i[!is.na(df_i$SOC_0_30) & !is.na(df_i$SOC_sd_0_30), ]\n",
    "  df_i$weight <- 1 / df_i$SOC_sd_0_30^2\n",
    "  df_i$weight <- df_i$weight / max(df_i$weight, na.rm = TRUE)\n",
    "  \n",
    "  dsm_models[[i]] <- ranger(\n",
    "    formula = my_formula, data = df_i,\n",
    "    case.weights = df_i$weight,\n",
    "    replace = FALSE, sample.fraction = 0.632,\n",
    "    num.trees = 500, seed = 201909\n",
    "  )\n",
    "  soc_maps[[i]] <- predict(covs, dsm_models[[i]], na.rm = TRUE)\n",
    "}\n",
    "\n",
    "\n",
    "################################################################################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate on the held-out set"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "################ Evaluate on the held-out set ##################################\n",
    "################################################################################\n",
    "\n",
    "dsm_metrics <- data.frame()\n",
    "\n",
    "for (i in seq_along(n_profiles)) {\n",
    "  preds_val <- predict(dsm_models[[i]], data = val_data)$predictions\n",
    "  \n",
    "  obs <- val_data$SOC_0_30\n",
    "  resid <- preds_val - obs\n",
    "  \n",
    "  me <- mean(resid, na.rm = TRUE)\n",
    "  rmse <- sqrt(mean(resid^2, na.rm = TRUE))\n",
    "  nse <- 1 - sum(resid^2, na.rm = TRUE) /\n",
    "    sum((obs - mean(obs, na.rm = TRUE))^2, na.rm = TRUE)\n",
    "  \n",
    "  dsm_metrics <- rbind(\n",
    "    dsm_metrics,\n",
    "    data.frame(\n",
    "      n_profiles = n_profiles[i],\n",
    "      n_samples = length(cal_indices[[i]]),\n",
    "      ME = round(me, 3),\n",
    "      RMSE = round(rmse, 3),\n",
    "      NSE = round(nse, 3)\n",
    "    )\n",
    "  )\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"
   ]
  }
 ]
}