{ "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 2 · Three-stage sampling design with covariate space coverage\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/2_sampling_design_csc_soilfer.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(\"RColorBrewer\", \"Rfast\", \"data.table\", \"doSNOW\", \"dplyr\", \"entropy\", \"fields\", \"ggplot2\", \"manipulate\", \"maps\", \"rassta\", \"raster\", \"sf\", \"sgsR\", \"snowfall\", \"sp\", \"stringr\", \"terra\", \"tibble\", \"tidyterra\", \"tripack\")\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(\"lemuscanovas/synoptReg\"))\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": [ "### MULTI-STAGE SAMPLING DESIGN FOR SOIL SAMPLING - Covariate Space Coverage" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "################################################################################\n", "# MULTI-STAGE SAMPLING DESIGN FOR SOIL SAMPLING - Covariate Space Coverage\n", "# \n", "# This script creates a three-stage sampling design:\n", "# - Primary Sampling Units (PSUs): 2x2 km grids\n", "# - Secondary Sampling Units (SSUs): 100x100 m cells within PSUs\n", "# - Tertiary Sampling Units (TSUs): Point locations within SSUs\n", "#\n", "# Author: Luis Rodriguez-Lado, PhD \n", "# Wanderson de Sousa Mendes, PhD\n", "# Date: 24 April 2026\n", "# Version: 3.0\n", "################################################################################\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1 - SET ENVIRONMENT AND LOAD LIBRARIES" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 1 - SET ENVIRONMENT AND LOAD LIBRARIES ======================================\n", "# Purpose: Load all required packages and set working directory\n", "\n", "# Verify working directory\n", "getwd()\n", "\n", "# List of required packages\n", "packages <- c(\"sp\",\"terra\",\"raster\",\"sf\", \"sgsR\",\"entropy\", \"tripack\",\"tibble\",\n", " \"manipulate\",\"dplyr\",\"synoptReg\", \"doSNOW\",\"Rfast\",\"fields\", \n", " \"ggplot2\", \"rassta\", \"snowfall\", \"maps\", \"tidyterra\",\n", " \"data.table\", \"RColorBrewer\", \"stringr\")\n", "\n", "# Load all packages silently\n", "invisible(lapply(packages, library, character.only = TRUE))\n", "rm(packages)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2 - DEFINE VARIABLES AND PARAMETERS" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 2 - DEFINE VARIABLES AND PARAMETERS =========================================\n", "# Purpose: Set all study-specific parameters in one place for easy modification\n", "\n", "# Country identification\n", "ISO.code <- \"KANSAS\" # 3-letter country code or full-name for file naming\n", "\n", "# Land use types to analyze (run script separately for each)\n", "landuse <- \"cropland\" # Options: \"cropland\", \"grassland\", \"forest\"\n", "\n", "# File paths\n", "raster.path <- \"01_data/module2/rasters/\"\n", "shp.path <- \"01_data/module2/shapes/\"\n", "other.path <- \"01_data/module2/other/\"\n", "results_dir <- paste0(\"03_outputs/module2/\")\n", "landuse_dir <- paste0(\"03_outputs/module2/\", landuse, \"/\")\n", "\n", "# Create directories if they don't exist\n", "dir.create(other.path, showWarnings = FALSE, recursive = TRUE)\n", "dir.create(\"03_outputs/module2\", showWarnings = FALSE, recursive = TRUE)\n", "dir.create(\"03_outputs/module2/img\", showWarnings = FALSE, recursive = TRUE)\n", "if (!file.exists(landuse_dir)) dir.create(landuse_dir)\n", "results.path <- landuse_dir\n", "\n", "# Coordinate reference system\n", "epsg <- \"EPSG:26714\" # NAD27 UTM Zone 14N - ADJUST FOR YOUR STUDY AREA\n", "\n", "# Sampling unit sizes (in meters)\n", "psu_size <- 2000 # PSU: 2 km x 2 km\n", "ssu_size <- 100 # SSU: 100 m x 100 m (1 hectare)\n", "\n", "# Number of SSUs per PSU\n", "num_primary_ssus <- 4 # Target SSUs (one per cluster)\n", "num_alternative_ssus <- 4 # Replacement SSUs (one per cluster)\n", "\n", "# Number of TSUs per SSU\n", "number_TSUs <- 3 # Point samples within each SSU\n", "\n", "# Optimization parameters\n", "iterations <- 10 # K-means clustering iterations\n", "\n", "# Minimum % of target land use required within a PSU for it to be selected\n", "percent_min <- 20 # for any of the target landuses selected PSUs must have >20% coverage on that target\n", "\n", "# Allocate number of sampling sites by land use type\n", "# ADJUST THESE PROPORTIONS based on your study objectives\n", "crop_prop <- 0.85\n", "grass_prop <- 0.10\n", "forest_prop <- 0.05\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3 - DEFINE CUSTOM FUNCTIONS" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 3 - DEFINE CUSTOM FUNCTIONS =================================================\n", "# Purpose: Create reusable functions for sampling design\n", "\n", "# Covariate Space Coverage with Legacy Data\n", "## This function performs constrained k-means clustering\n", "CSIS <- function(fixed, nsup, nstarts, mygrd) {\n", " # Args:\n", " # fixed: Data frame of fixed legacy points\n", " # nsup: Number of supplementary points to select\n", " # nstarts: Number of random starts for optimization\n", " # mygrd: Grid of all potential sampling locations\n", " \n", " n_fix <- nrow(fixed)\n", " p <- ncol(mygrd)\n", " units <- fixed$units\n", " mygrd_minfx <- mygrd[-units, ]\n", " MSSSD_cur <- NA\n", " \n", " for (s in 1:nstarts) {\n", " units <- sample(nrow(mygrd_minfx), nsup)\n", " centers_sup <- mygrd_minfx[units, ]\n", " centers <- rbind(fixed[, names(mygrd)], centers_sup)\n", " \n", " repeat {\n", " D <- rdist(x1 = centers, x2 = mygrd)\n", " cluster <- apply(X = D, MARGIN = 2, FUN = which.min) %>% as.factor(.)\n", " centers_cur <- centers\n", " \n", " for (i in 1:p) {\n", " centers[, i] <- tapply(mygrd[, i], INDEX = cluster, FUN = mean)\n", " }\n", " \n", " # Restore fixed centers (legacy data points don't move)\n", " centers[1:n_fix, ] <- centers_cur[1:n_fix, ]\n", " \n", " # Check convergence\n", " sumd <- diag(rdist(x1 = centers, x2 = centers_cur)) %>% sum(.)\n", " if (sumd < 1E-12) {\n", " D <- rdist(x1 = centers, x2 = mygrd)\n", " Dmin <- apply(X = D, MARGIN = 2, FUN = min)\n", " MSSSD <- mean(Dmin^2)\n", " \n", " if (s == 1 | MSSSD < MSSSD_cur) {\n", " centers_best <- centers\n", " clusters_best <- cluster\n", " MSSSD_cur <- MSSSD\n", " }\n", " break\n", " }\n", " }\n", " print(paste0(s,\" out of \",nstarts))\n", " }\n", " list(centers = centers_best, cluster = clusters_best)\n", "}\n", "\n", "# K-means with Progress Reporting\n", "## This function performs k-means clustering if legacy data is not available (progress reporting)\n", "kmeans_with_progress <- function(data, centers, iter.max = 10000, nstart = 100) {\n", " # Provides visual feedback during long k-means operations\n", " best_result <- NULL\n", " best_totss <- Inf\n", " \n", " cat(\"Running k-means clustering with\", nstart, \"random starts...\\n\")\n", " \n", " for (s in 1:nstart) {\n", " result <- kmeans(data, centers = centers, iter.max = iter.max, nstart = 1)\n", " \n", " if (result$tot.withinss < best_totss) {\n", " best_result <- result\n", " best_totss <- result$tot.withinss\n", " }\n", " \n", " print(paste0(s, \" out of \", nstart))\n", " }\n", " \n", " cat(\"Best tot.withinss:\", best_totss, \"\\n\")\n", " return(best_result)\n", "}\n", "\n", "# Generate TSU Points Within SSU\n", "## This function creates random point samples within each SSU polygon\n", "generate_tsu_points_within_ssu <- function(ssu, number_TSUs, index, ssu_type, crops) {\n", " # Args:\n", " # ssu: Single SSU polygon (sf object)\n", " # number_TSUs: Number of points to generate\n", " # index: SSU identifier\n", " # ssu_type: \"Target\" or \"Replacement\"\n", " # crops: Crop mask raster (20m resolution)\n", " \n", " ssu_vect <- vect(ssu)\n", " \n", " # Validate geometry\n", " if (is.null(ssu_vect) || nrow(ssu_vect) == 0 || is.na(ext(ssu_vect))) {\n", " warning(paste(\"SSU\", index, \"has invalid geometry. Skipping TSU generation.\"))\n", " return(NULL)\n", " }\n", " \n", " # Clip crop raster to SSU\n", " clipped_lu <- try(crop(crops, ssu_vect), silent = TRUE)\n", " if (inherits(clipped_lu, \"try-error\") || is.null(clipped_lu)) {\n", " warning(paste(\"SSU\", index, \"could not crop land use raster. Skipping.\"))\n", " return(NULL)\n", " }\n", " \n", " # Sample points (tries two methods)\n", " sampled_points <- try(sample_srs(clipped_lu, nSamp = number_TSUs), silent = TRUE)\n", " if (inherits(sampled_points, \"try-error\") || is.null(sampled_points) || nrow(sampled_points) == 0) {\n", " sampled_points <- try(spatSample(clipped_lu, size = number_TSUs, na.rm = TRUE, method = \"random\"), silent = TRUE)\n", " }\n", " \n", " if (inherits(sampled_points, \"try-error\") || is.null(sampled_points) || nrow(sampled_points) == 0) {\n", " warning(paste(\"SSU\", index, \"failed to generate TSUs. Skipping.\"))\n", " return(NULL)\n", " }\n", " \n", " # Add metadata\n", " sampled_points$PSU_ID <- selected_psu$ID\n", " sampled_points$SSU_ID <- index\n", " sampled_points$TSU_ID <- seq_len(nrow(sampled_points))\n", " sampled_points$SSU_Type <- ssu_type\n", " sampled_points$TSU_Name <- paste0(sampled_points$PSU_ID, \".\", index, \".\", seq_len(nrow(sampled_points)))\n", " \n", " return(sampled_points)\n", "}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4 - LOAD COUNTRY BOUNDARIES AND LEGACY DATA" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 4 - LOAD COUNTRY BOUNDARIES AND LEGACY DATA =================================\n", "# Purpose: Import study area boundaries and existing soil sample locations\n", "\n", "# Load country/region boundaries\n", "country_boundaries <- file.path(paste0(shp.path,\"roi_kansas_us_epsg_4326.shp\"))\n", "country_boundaries <- sf::st_read(country_boundaries, quiet=TRUE)\n", "\n", "# Reproject if necessary\n", "if(!same.crs(country_boundaries, epsg)){\n", " country_boundaries <- country_boundaries %>%\n", " st_as_sf() %>% sf::st_transform(crs=epsg)\n", "}\n", "\n", "# Load legacy soil data (optional - existing sample points)\n", "legacy <- file.path(paste0(shp.path,\"soil_legacy_data_kansas_epsg_4326.shp\"))\n", "\n", "if(file.exists(legacy)){\n", " legacy <- sf::st_read(legacy, quiet=TRUE)\n", " if(!same.crs(legacy, epsg)){\n", " legacy <- legacy %>% sf::st_transform(crs=epsg)\n", " }\n", "} else {\n", " rm(legacy) # Remove if doesn't exist\n", "} \n", "\n", "# Clean legacy data\n", "if(exists(\"legacy\")){\n", " legacy <- dplyr::select(legacy, geometry)\n", " legacy <- legacy[!duplicated(st_geometry(legacy)), ] # Remove duplicates\n", "}\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/legacy_map.png\")),\n", " width = 4000, # pixels\n", " height = 2800, # pixels\n", " res = 400 # DPI\n", ")\n", "# Visualize boundaries and legacy points\n", "ggplot() +\n", " geom_spatvector(data = country_boundaries, fill = NA, color = \"black\") +\n", " geom_spatvector(data = legacy, aes(geometry = geometry), size = 0.7, color = \"red\") +\n", " theme_minimal()\n", "\n", "dev.off()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5 - LOAD OPTIONAL EXCLUSION LAYERS" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 5 - LOAD OPTIONAL EXCLUSION LAYERS ==========================================\n", "# Purpose: Load masks to exclude protected areas, steep slopes, etc.\n", "\n", "# Protected areas (areas to EXCLUDE from sampling)\n", "npa <- file.path(paste0(shp.path,\"protected_areas_epsg_4326.shp\"))\n", "\n", "if(file.exists(npa)){\n", " npa <- sf::st_read(npa, quiet = FALSE)\n", " if(!same.crs(npa, epsg)){\n", " npa <- npa %>% sf::st_transform(crs = epsg)\n", " }\n", " npa <- sf::st_union(npa)\n", " npa <- sf::st_difference(country_boundaries, npa) # Create \"non-protected\" mask\n", "} else {\n", " rm(npa)\n", "}\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/non_protected_map.png\")),\n", " width = 4000, # pixels\n", " height = 2800, # pixels\n", " res = 400 # DPI\n", ")\n", "\n", "ggplot() +\n", " geom_sf(data = npa, fill = \"grey90\", color = \"black\", linewidth = 0.3) +\n", " coord_sf(datum = NA) +\n", " labs(title = \"Accessible area mask (non-protected in grey)\") +\n", " theme_minimal(base_size = 11)\n", "\n", "dev.off()\n", "\n", "# Slope mask (exclude areas with slope > threshold)\n", "# This should be a BINARY RASTER where:\n", " # - Value = 1 : Accessible slopes (≤ threshold, e.g., ≤30%)\n", " # - Value = NA: Inaccessible slopes (> threshold, excluded from sampling)\n", " #\n", " # HOW TO CREATE IN QGIS:\n", " # Step 1: Analyze your slope raster to determine appropriate threshold\n", " # (Raster → Raster Calculator or Slope tool from DEM)\n", " # Step 2: Open Raster Calculator (Raster → Raster Calculator)\n", " # Step 3: Use this formula to create binary mask:\n", " # ( \"Slope@1\" <= 30 ) * 1\n", " # Replace 30 with your chosen threshold (%)\n", " # Step 4: This outputs: 1 where slope ≤30%, 0 where slope >30%\n", " # Step 5: (Optional) Convert 0 to NA: ( \"Slope@1\" <= 30 ) * 1 + ( \"Slope@1\" > 30 ) * -9999\n", " # Then use \"Set Null\" tool to convert -9999 to NA\n", " #\n", " # THRESHOLD GUIDELINES:\n", " # - Gentle terrain: 0-15%\n", " # - Moderate terrain: 15-30% (typical threshold for field work)\n", " # - Steep terrain: >30% (usually excluded)\n", "slope <- file.path(paste0(raster.path,\"slope_mask_epsg_4326.tif\"))\n", "\n", "if(file.exists(slope)){\n", " slope <- rast(slope)\n", " if(!same.crs(slope, epsg)){\n", " slope <- project(slope, epsg, method=\"near\")\n", " }\n", " slope <- slope/slope # Convert to binary mask\n", " slope <- terra::mask(slope, country_boundaries)\n", "} else {\n", " rm(slope)\n", "}\n", "\n", "png(file.path(paste0(results_dir,\"img/slope_mask.png\")),\n", " width = 4000,\n", " height = 2800,\n", " res = 400,\n", " type = \"cairo\")\n", "\n", "plot(slope,\n", " col = \"#66c2a5\",\n", " legend = FALSE,\n", " main = \"Slope accessibility mask in green\")\n", "\n", "dev.off()\n", "\n", "# Geology data (for stratification) - If available\n", "geo <- file.path(paste0(shp.path,\"ecoregions_kansas_epsg_4326.shp\"))\n", "geo.classes <- \"US_L3NAME\" # Field name for geology classes\n", "\n", "if(file.exists(geo)){\n", " geo <- sf::st_read(geo, quiet=TRUE)\n", " if(!same.crs(geo, epsg)){\n", " geo <- geo %>% sf::st_transform(crs=epsg)\n", " }\n", " geo$GEO <- as.numeric(as.factor(geo[[geo.classes]]))\n", "} else {\n", " rm(geo)\n", "}\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/geology_ecoregion_map.png\")),\n", " width = 4000, # pixels\n", " height = 2800, # pixels\n", " res = 400 # DPI\n", ")\n", "\n", "ggplot() +\n", " geom_sf(data = geo, aes(fill = US_L3NAME), color = \"black\", linewidth = 0.2) +\n", " labs(title = \"Level III Ecoregions\",\n", " fill = \"Ecoregion (US_L3NAME)\") +\n", " theme_minimal(base_size = 11) +\n", " theme(legend.position = \"right\")\n", "\n", "dev.off()\n", "\n", "# Geomorphology data\n", "geomorph <- file.path(paste0(raster.path,\"/Geomorphon_Landforms_KANSAS.tif\"))\n", "geomorph.classes <- \"Class\"\n", "\n", "if (file.exists(geomorph)) {\n", " file_extension <- tools::file_ext(geomorph)\n", " \n", " if (file_extension == \"tif\") {\n", " geomorph <- rast(geomorph)\n", " names(geomorph) <- 'GEOMORPH'\n", " if(!same.crs(geomorph, epsg)){\n", " geomorph <- project(geomorph, epsg, method=\"near\")\n", " }\n", " } else if (file_extension == \"shp\") {\n", " geomorph <- sf::st_read(geomorph, quiet = TRUE)\n", " if (sf::st_crs(geomorph)$epsg != epsg) {\n", " geomorph <- sf::st_transform(geomorph, crs = epsg)\n", " }\n", " geomorph$GEOMORPH <- as.numeric(as.factor(geomorph[[geomorph.classes]]))\n", " }\n", " \n", " geomorph <- terra::mask(geomorph, country_boundaries)\n", "} else {\n", " rm(geomorph)\n", "}\n", "\n", "# Visualisation\n", "# Create a class lookup table (Geomorpho90m / geomorphons)\n", "geomorph_lut <- data.frame(\n", " value = 1:10,\n", " class = c(\n", " \"Flat\",\n", " \"Peak / summit\",\n", " \"Ridge\",\n", " \"Shoulder\",\n", " \"Spur\",\n", " \"Slope\",\n", " \"Hollow\",\n", " \"Footslope\",\n", " \"Valley\",\n", " \"Pit / depression\"\n", " )\n", ")\n", "\n", "# Semantically meaningful colors: warm tones for highs, cool tones for lows\n", "geomorph_colors <- c(\n", " \"Flat\" = \"#F5F5DC\", # beige – level ground\n", " \"Peak / summit\" = \"#8B0000\", # dark red – highest points\n", " \"Ridge\" = \"#CD5C5C\", # indian red – elongated highs\n", " \"Shoulder\" = \"#D2691E\", # chocolate – convex upper slopes\n", " \"Spur\" = \"#DAA520\", # goldenrod – diverging slopes\n", " \"Slope\" = \"#6B8E23\", # olive – planar slopes\n", " \"Hollow\" = \"#4682B4\", # steel blue – converging slopes\n", " \"Footslope\" = \"#5F9EA0\", # cadet blue – concave lower slopes\n", " \"Valley\" = \"#00008B\", # dark blue – linear lows\n", " \"Pit / depression\" = \"#191970\" # midnight – lowest points\n", ")\n", "\n", "# Attach labels to the raster as categories\n", "geomorph_cat <- as.factor(geomorph)\n", "levels(geomorph_cat) <- list(data.frame(\n", " ID = geomorph_lut$value,\n", " class = geomorph_lut$class\n", "))\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/geomorph_classes.png\")),\n", " width = 4000,\n", " height = 2800,\n", " res = 400\n", ")\n", "\n", "par(mar = c(3, 3, 4, 3), bg = \"white\")\n", "\n", "plot(\n", " geomorph_cat,\n", " col = geomorph_colors[levels(geomorph_cat)[[1]]$class],\n", " main = \"Geomorphology\",\n", " axes = FALSE,\n", " box = FALSE,\n", " plg = list(\n", " inset = c(0.02, 0.03),\n", " cex = 0.75,\n", " title = \"Geomorphons\",\n", " bty = \"n\"\n", " )\n", ")\n", "\n", "dev.off()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6 - LOAD AND PROCESS ENVIRONMENTAL COVARIATES (PSU LEVEL)" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 6 - LOAD AND PROCESS ENVIRONMENTAL COVARIATES (PSU LEVEL) ===================\n", "# Purpose: Load environmental data at 2km resolution for PSU selection from GEE code\n", "# Note: These covariates determine WHERE PSUs are placed\n", "\n", "cov.dat <- list.files(raster.path,\"Environmental_Covariates_250m_KANSAS.tif\", \n", " recursive = TRUE, full.names = TRUE)\n", "cov.dat <- terra::rast(cov.dat)\n", "\n", "# Reproject if necessary\n", "if(!same.crs(cov_dat, epsg)){\n", " cov.dat <- terra::project(cov.dat, epsg, method=\"near\")\n", "}\n", "\n", "# Resample to PSU resolution (2km)\n", "# IMPORTANT: This aggregation must match your PSU size\n", "psu_size_template <- rast(ext(cov.dat), resolution = psu_size, crs = crs(cov.dat))\n", "cov.dat <- resample(cov.dat, psu_size_template, method = \"bilinear\")\n", "cov.dat <- terra::mask(cov.dat, country_boundaries)\n", "names(cov.dat)\n", "\n", "# Load and process soil climate data\n", "newhall <- list.files(raster.path, pattern = \"newhall.tif$\", recursive = TRUE, full.names = TRUE)\n", "newhall <- terra::rast(newhall)\n", "\n", "if(crs(newhall)!=epsg){\n", " newhall <- terra::project(newhall, epsg, method=\"near\")\n", "}\n", "\n", "# Remove unnecessary layers\n", "newhall$regimeSubdivision1 <- NULL\n", "newhall$regimeSubdivision2 <- NULL\n", "\n", "# Process categorical variables (preserve factor levels)\n", "temperatureRegime <- project(newhall$temperatureRegime, cov.dat, method = \"near\")\n", "moistureRegime <- project(newhall$moistureRegime, cov.dat, method = \"near\")\n", "\n", "newhall$temperatureRegime <- NULL\n", "newhall$moistureRegime <- NULL\n", "newhall <- terra::resample(newhall, cov.dat)\n", "\n", "# Convert to dummy variables (one column per category)\n", "temperatureRegime <- as.factor(temperatureRegime)\n", "temperatureRegime <- dummies(ca.rast = temperatureRegime, preval = 1, absval = 0)\n", "moistureRegime <- as.factor(moistureRegime)\n", "moistureRegime <- dummies(ca.rast = moistureRegime, preval = 1, absval = 0)\n", "\n", "# Merge all covariates\n", "cov.dat <- c(cov.dat, newhall, temperatureRegime, moistureRegime)\n", "\n", "# Add geology if available\n", "if(exists(\"geo\")){\n", " geo <- rasterize(as(geo,\"SpatVector\"), cov.dat, field=\"GEO\")\n", " geo <- dummies(ca.rast = geo$GEO, preval = 1, absval = 0)\n", " cov.dat <- c(cov.dat, geo)\n", "}\n", "\n", "# Add geomorphology if available\n", "if (exists(\"geomorph\")) {\n", " if (!inherits(geomorph, \"SpatRaster\")) {\n", " geomorph <- rasterize(as(geomorph, \"SpatVector\"), cov.dat, field = \"GEOMORPH\")\n", " }\n", " geomorph <- dummies(ca.rast = geomorph$GEOMORPH, preval = 1, absval = 0)\n", " \n", " # Ensure matching extents\n", " if (!identical(ext(cov.dat), ext(geomorph))) {\n", " geomorph <- extend(geomorph, cov.dat)\n", " }\n", " if (!all(res(cov.dat) == res(geomorph))) {\n", " geomorph <- resample(geomorph, cov.dat, method = \"near\")\n", " }\n", " \n", " cov.dat <- c(cov.dat, geomorph)\n", "}\n", "\n", "# Clean up\n", "rm(newhall, geomorph)\n", "gc()\n", "\n", "# Crop to study area\n", "cov.dat <- crop(cov.dat, country_boundaries, mask=TRUE, overwrite=TRUE)\n", "writeRaster(cov.dat, file.path(paste0(results_dir,\"cov_dat_stack_psus.tif\")), overwrite=TRUE)\n", "\n", "names(cov.dat)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 7 - DIMENSIONALITY REDUCTION WITH PCA" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 7 - DIMENSIONALITY REDUCTION WITH PCA =======================================\n", "# Purpose: Reduce many covariates to fewer principal components\n", "# Why: Makes clustering faster and reduces multicollinearity\n", "\n", "pca <- scale(cov.dat)\n", "pca <- synoptReg::raster_pca(pca) # Fast PCA for rasters\n", "\n", "cov.dat <- pca$PCA\n", "\n", "# Keep only components explaining 99% of variance\n", "n_comps <- first(which(pca$summaryPCA[3,] > 0.99))\n", "cov.dat <- pca$PCA[[1:n_comps]]\n", "\n", "cat(sprintf(\"Using %d principal components (explaining >99%% variance)\\n\", n_comps))\n", "\n", "# Save PCA results\n", "writeRaster(cov.dat, file.path(paste0(results_dir,\"PCA_projected.tif\")), overwrite=TRUE)\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/cum_var_pcas_map.png\")),\n", " width = 4000, # pixels\n", " height = 2800, # pixels\n", " res = 400 # DPI\n", ")\n", "\n", "cum_var <- pca$summaryPCA[\"Cumulative\", ]\n", "\n", "plot(cum_var,\n", " type = \"l\",\n", " lwd = 2,\n", " xlab = \"Principal Component\",\n", " ylab = \"Cumulative Variance Explained\",\n", " main = \"Cumulative Variance\")\n", "\n", "abline(h = 0.99, col = \"red\", lty = 2)\n", "\n", "dev.off()\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/pcas_map.png\")),\n", " width = 4000, # pixels\n", " height = 2800, # pixels\n", " res = 400 # DPI\n", ")\n", "\n", "# Consistent color scale across PC1–PC3\n", "minmax_vals <- minmax(pca$PCA[[1:3]])\n", "zlim_vals <- range(minmax_vals)\n", "\n", "cols <- hcl.colors(100, \"Blue-Red 3\", rev = TRUE)\n", "\n", "par(mfrow = c(1,3), mar = c(3,3,3,6)) # extra space for legend\n", "\n", "plot(pca$PCA[[1]],\n", " col = cols,\n", " zlim = zlim_vals,\n", " main = \"PC1\",\n", " axes = FALSE,\n", " box = FALSE,\n", " plg = list(title = \"PC value\"))\n", "\n", "plot(pca$PCA[[2]],\n", " col = cols,\n", " zlim = zlim_vals,\n", " axes = FALSE,\n", " box = FALSE,\n", " main = \"PC2\",\n", " plg = list(title = \"PC value\"))\n", "\n", "plot(pca$PCA[[3]],\n", " col = cols,\n", " zlim = zlim_vals,\n", " main = \"PC3\",\n", " axes = FALSE,\n", " box = FALSE,\n", " plg = list(title = \"PC value\"))\n", "\n", "dev.off()\n", "\n", "par(mfrow = c(1,1))\n", "\n", "# Remove pca\n", "rm(pca)\n", "\n", "# Reload for further processing\n", "# Why this is useful:\n", "# - Previous sections (data loading, PCA) can take 30+ minutes\n", "# - If script crashes or needs modification, you can load here\n", "# - If the file exists in the folder from previous processing\n", "cov.dat <- rast(file.path(paste0(results_dir,\"PCA_projected.tif\")))\n", "# plot(cov.dat[[1]])\n", "\n", "# Reproject if necessary\n", "if(!same.crs(cov.dat, epsg)){\n", " cov.dat <- terra::project(cov.dat, epsg, method=\"near\")\n", "}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 8 - LOAD AND PREPARE LAND USE DATA" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 8 - LOAD AND PREPARE LAND USE DATA ==========================================\n", "# Purpose: Define sampling universe (where samples CAN be taken)\n", "# Uses TWO resolutions: 20m for TSU placement, 100m for PSU filtering\n", "\n", "# Load cropland mask\n", "landuse_file <- file.path(paste0(raster.path,\"Cropland_Mask_KANSAS_merged.tif\"))\n", "crops <- rast(landuse_file)\n", "crops <- crops/crops # Convert to binary (1=crop, NA=other)\n", "names(crops) <- \"lu\"\n", "\n", "# Reproject if necessary\n", "if(!same.crs(crops, epsg)){\n", " crops <- terra::project(crops, epsg, method=\"near\")\n", "}\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/crop_20m_map.png\")),\n", " width = 4000, # pixels\n", " height = 2800, # pixels\n", " res = 400 # DPI\n", ")\n", "# Visualize (it takes a few minutes)\n", "ggplot() +\n", " geom_spatraster(data = as.factor(crops)) +\n", " scale_fill_viridis_d(na.value = \"transparent\") +\n", " geom_spatvector(data = country_boundaries, fill = NA, color = \"black\") +\n", "# geom_spatvector(data = legacy, size = 0.7, color = \"red\") +\n", " theme_minimal()\n", "\n", "dev.off()\n", "\n", "# Resample to 20m resolution (for TSU placement)\n", "# CRITICAL: This resolution determines precision of TSU point placement\n", "lulc_size_template <- rast(ext(crops), resolution = 20, crs = crs(crops))\n", "crops <- as.factor(crops)\n", "crops <- resample(crops, lulc_size_template, method = \"near\")\n", "\n", "# Apply exclusion masks\n", "if(exists(\"npa\")){\n", " crops <- mask(crops, npa)\n", "}\n", "\n", "if(exists(\"slope\")){\n", " slope <- resample(slope, crops, method=\"near\")\n", " crops <- crops * slope\n", "}\n", "\n", "rm(npa, slope)\n", "\n", "# Save 20m resolution crop mask\n", "writeRaster(crops, file.path(paste0(results_dir,\"crop_mask_20m_clean.tif\")), overwrite=TRUE)\n", "# Same situation as loading \"cov.dat\"\n", "crops <- rast(file.path(paste0(results_dir,\"crop_mask_20m_clean.tif\")))\n", "\n", "# Create 100m resolution version for PSU filtering\n", "# Aggregate: (20m × 5) = 100m pixels\n", "lu <- aggregate(crops, 5, fun=modal, cores = 2, na.rm=T)\n", "names(lu) <- \"lu\"\n", "writeRaster(lu, file.path(paste0(results_dir,\"crop_mask_100m.tif\")), overwrite=TRUE)\n", "# Same situation as loading \"cov.dat\"\n", "lu <- rast(paste0(raster.path,\"crop_mask_100m.tif\"))\n", "\n", "# Filter legacy data to crop areas\n", "if (exists(\"legacy\")){\n", " legacy$INSIDE <- terra::extract(crops, legacy) %>% dplyr::select(lu)\n", " legacy <- legacy[!is.na(legacy$INSIDE),] %>% dplyr::select(-\"INSIDE\")\n", "}\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/crop_100m_legacy_map.png\")),\n", " width = 4000, # pixels\n", " height = 2800, # pixels\n", " res = 400 # DPI\n", ")\n", "\n", "# Visualize (it takes a few minutes)\n", "ggplot() +\n", " geom_spatraster(data = as.factor(lu)) +\n", " scale_fill_viridis_d(na.value = \"transparent\") +\n", " geom_spatvector(data = country_boundaries, fill = NA, color = \"black\") +\n", " geom_spatvector(data = legacy, size = 0.7, color = \"red\") +\n", " theme_minimal()\n", "\n", "dev.off()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 9 - GENERATE PSU GRID" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 9 - GENERATE PSU GRID ========================================================\n", "# Purpose: Create 2×2 km grid covering the study area\n", "\n", "psu_grid <- st_make_grid(country_boundaries, cellsize = c(psu_size, psu_size), square = TRUE)\n", "psu_grid <- st_sf(geometry = psu_grid)\n", "psu_grid$ID <- 1:nrow(psu_grid)\n", "\n", "# Clip to country boundary\n", "psu_grid <- psu_grid[country_boundaries[1],] # TIME CONSUMING!\n", "write_sf(psu_grid, file.path(paste0(results_dir,\"grid2k.shp\")), overwrite=TRUE)\n", "\n", "# Or load pre-saved grid (much faster)\n", "# psu_grid <- sf::st_read(file.path(paste0(results_dir,\"grid2k.shp\")))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 10 - FILTER PSUs BY CROP COVERAGE" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 10 - FILTER PSUs BY CROP COVERAGE ===========================================\n", "# Purpose: Keep only PSUs with sufficient cropland\n", "# Why: No point sampling in PSUs with <20% crops\n", "\n", "# Extract crop percentage for each PSU\n", "extracted_values <- terra::extract(lu, psu_grid)\n", "\n", "crop_perc <- extracted_values %>%\n", " group_by(ID) %>%\n", " summarize(crop_perc = sum(lu, na.rm = TRUE)*100/400) # 400 = total of 100m pixels in 2km PSU\n", "\n", "rm(extracted_values)\n", "\n", "# Join back to PSU grid\n", "psu_grid$crop_perc <- crop_perc$crop_perc\n", "write_sf(psu_grid, file.path(paste0(results_dir,\"psu_grid_counts.shp\")), overwrite=TRUE)\n", "\n", "# Reload and ensure correct projection\n", "## Same as cov.dat and land use data\n", "psu_grid <- sf::st_read(file.path(paste0(results_dir,\"psu_grid_counts.shp\")))\n", "if(!same.crs(psu_grid, epsg)){\n", " psu_grid <- psu_grid %>% sf::st_transform(crs=epsg)\n", "}\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/psu_grid_lu_coverage.png\")),\n", " width = 4000, # pixels\n", " height = 2800, # pixels\n", " res = 400 # DPI\n", ")\n", "\n", "# Visualize crop coverage\n", "ggplot() +\n", " geom_sf(data = psu_grid, aes(fill = crop_perc)) +\n", " scale_fill_distiller(palette = \"Spectral\") +\n", " labs(title = \"Crop Coverage by PSU\", fill = \"% Cropland\") +\n", " theme_minimal()\n", "\n", "dev.off()\n", "\n", "# Filter: Keep only PSUs with n > percent_min coverage\n", "psu_grid <- psu_grid[psu_grid$crop_perc > percent_min, \"ID\"]\n", "\n", "cat(sprintf(\"Retained %d PSUs with > %d%% crop coverage\\n\", \n", " nrow(psu_grid), percent_min))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 11 - RASTERIZE PSU GRID" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 11 - RASTERIZE PSU GRID =====================================================\n", "# Purpose: Convert vector PSUs to raster format for covariate extraction\n", "\n", "template <- rast(vect(psu_grid), res = psu_size)\n", "template <- rasterize(vect(psu_grid), template, field = \"ID\")\n", "\n", "# Crop covariates to eligible PSUs only\n", "cov.dat <- crop(cov.dat, psu_grid, mask=TRUE)\n", "PSU.r <- resample(cov.dat, template)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 12 - CALCULATE OPTIMAL SAMPLE SIZE" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 12 - CALCULATE OPTIMAL SAMPLE SIZE ==========================================\n", "# Purpose: Determine how many PSUs needed for representative coverage\n", "# What this does:\n", "# - Tests different sample sizes (50, 75, 100, ... up to 3000 PSUs)\n", "# - Measures how well each sample size represents environmental variability\n", "# - Identifies the \"sweet spot\" where adding more PSUs gives diminishing returns\n", "# - Uses Kullback-Leibler divergence to compare sample vs population distributions\n", "#\n", "# Method: Feature Space Coverage (FCS) algorithm\n", "# - Iteratively samples PSUs and compares to full covariate space\n", "# - Repeats 4 times per sample size to ensure stability\n", "# - Selects optimal N where coverage reaches 95% of maximum\n", "#\n", "# COMPUTATION TIME: \n", "# - Small areas (<1000 PSUs): 2-6 hours\n", "# - Medium areas (1000-5000 PSUs): 6-24 hours \n", "# - Large areas (>5000 PSUs): 1-3 days\n", "# - Depends on: # of PSUs, # of covariates, CPU cores available\n", "#\n", "# SKIP THIS SECTION IF:\n", "# - You already ran it and saved the result (optimal_N_KLD.RDS exists)\n", "# - You have a predetermined sample size (e.g., budget constraints)\n", "# - Running quick tests (use arbitrary N like 50-100 PSUs)\n", "#\n", "\n", "source(\"02_scripts/module2/opt_sample.R\") # Load optimization functions\n", "\n", "psu.r.df <- data.frame(PSU.r)\n", "\n", "# Optimization parameters\n", "initial.n <- 50\n", "final.n <- 3000\n", "by.n <- 25\n", "iters <- 4\n", "\n", "# Run optimization (can take several minutes)\n", "opt_N_fcs <- opt_sample(alg=\"fcs\",\n", " s_min=initial.n,\n", " s_max=final.n,\n", " s_step=by.n,\n", " s_reps=iters,\n", " covs = psu.r.df,\n", " cpus=4,\n", " conf=0.95)\n", "\n", "optimal_N_KLD <- opt_N_fcs$optimal_sites[1,2]\n", "cat(sprintf(\"Optimal sample size: %d PSUs\\n\", optimal_N_KLD))\n", "\n", "saveRDS(optimal_N_KLD, paste0(results.path,\"../optimal_N_KLD.RDS\"))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 13 - SELECT PSUs USING COVARIATE SPACE COVERAGE" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 13 - SELECT PSUs USING COVARIATE SPACE COVERAGE =============================\n", "# Purpose: Choose PSU locations that maximize environmental diversity\n", "# Method: Constrained k-means clustering (respects legacy data)\n", "\n", "optimal_N_KLD <- readRDS(\"01_data/module2/optimal_nos/optimal_N_KLD.RDS\")\n", "\n", "# Calculate number of PSUs for this land use\n", "n.psu <- round(optimal_N_KLD * crop_prop, 0) # For cropland\n", "\n", "# OPTIONAL: Add buffer to account for skipped PSUs\n", "n.psu <- round(optimal_N_KLD * crop_prop * 1.15, 0) # 15% buffer\n", "\n", "cat(sprintf(\"Targeting %d PSUs for %s\\n\", n.psu, landuse))\n", "\n", "# Prepare data\n", "PSU.df <- as.data.frame(PSU.r, xy=T)\n", "covs <- names(cov.dat)\n", "mygrd <- data.frame(scale(PSU.df[, covs]))\n", "\n", "# If legacy data exists, use constrained clustering\n", "if (exists(\"legacy\")){\n", " legacy <- st_filter(legacy, psu_grid)\n", " legacy_df <- st_coordinates(legacy)\n", " \n", " # Find nearest PSU for each legacy point\n", " units <- numeric(nrow(legacy_df))\n", " for (i in 1:nrow(legacy_df)) {\n", " distances <- sqrt((PSU.df$x - legacy_df[i, \"X\"])^2 + (PSU.df$y - legacy_df[i, \"Y\"])^2)\n", " units[i] <- which.min(distances)\n", " }\n", " \n", " fixed <- unique(data.frame(units, scale(PSU.df[, covs])[units, ]))\n", " \n", " # Run constrained clustering (REVIEW)\n", " res <- CSIS(fixed = fixed, nsup = n.psu, nstarts = iterations, mygrd = mygrd)\n", " \n", "} else {\n", " # No legacy data: standard k-means\n", " res <- kmeans_with_progress(mygrd, centers = n.psu, iter.max = 10000, nstart = iterations)\n", "}\n", "\n", "# Assign cluster IDs\n", "PSU.df$cluster <- res$cluster\n", "\n", "# Find PSU closest to each cluster center (these become sample locations)\n", "D <- rdist(x1 = res$centers, x2 = scale(PSU.df[, covs]))\n", "units <- apply(D, MARGIN = 1, FUN = which.min)\n", "\n", "myCSCsample <- PSU.df[units, c(\"x\", \"y\", covs)]\n", "\n", "# Label legacy vs new PSUs\n", "if (exists(\"legacy\")){\n", " myCSCsample$type <- c(rep(\"legacy\", nrow(fixed)), rep(\"new\", length(units)-nrow(fixed)))\n", "} else {\n", " myCSCsample$type <- \"new\"\n", "}\n", "\n", "# Convert to spatial object\n", "myCSCsample <- myCSCsample %>%\n", " st_as_sf(coords = c(\"x\", \"y\"), crs = epsg)\n", "\n", "# Separate legacy and new\n", "if (exists(\"legacy\")){\n", " legacy <- myCSCsample[myCSCsample$type==\"legacy\",]\n", "}\n", "new <- myCSCsample[myCSCsample$type==\"new\",]\n", "\n", "# Extract target PSU IDs\n", "PSUs <- sf::st_intersection(psu_grid, new) %>% dplyr::select(ID)\n", "target.PSUs <- psu_grid[psu_grid$ID %in% PSUs$ID,] %>% dplyr::select(ID)\n", "\n", "cat(sprintf(\"Selected %d target PSUs\\n\", nrow(target.PSUs)))\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/csc_psu_distribution.png\")),\n", " width = 4000, # pixels\n", " height = 2800, # pixels\n", " res = 400 # DPI\n", ")\n", "\n", "# Visualize\n", "ggplot() +\n", " geom_raster(data = as.data.frame(PSU.r$PC1, xy = TRUE), aes(x = x, y = y, fill = PC1)) +\n", " scale_fill_viridis_c() +\n", " geom_sf(data = target.PSUs, color = \"#101010\", fill = NA, lwd = 0.8) +\n", " geom_sf(data = new[1], color = \"#D81B60\", size = 0.5, shape = 19) +\n", " labs(title = \"Target Primary Sampling Units\") +\n", " theme_minimal()\n", "\n", "dev.off()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 14 - LOAD HIGH-RESOLUTION COVARIATES (SSU LEVEL)" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 14 - LOAD HIGH-RESOLUTION COVARIATES (SSU LEVEL) ============================\n", "# Purpose: Load 100m resolution data for SSU clustering WITHIN each PSU\n", "# Note: Different from PSU covariates - these determine SSU placement\n", "\n", "cov.dat.ssu <- terra::rast(paste0(raster.path, \"HighRes_Covariates_100m_KANSAS_merged.tif\"))\n", "names(cov.dat.ssu) <- gsub(\"(^\\\\d+_?S2_|^\\\\d+_|^S2_)\", \"\", names(cov.dat.ssu))\n", "\n", "if(!same.crs(cov.dat.ssu, epsg)){\n", " cov.dat.ssu <- terra::project(cov.dat.ssu, epsg, method=\"near\")\n", "}\n", "\n", "# Check for and replace NA values\n", "# Why: NA values cause complete.cases() to remove SSUs unnecessarily\n", "# Solution: Replace NA with 0 (assumes missing data = no feature present)\n", "if (any(is.na(values(cov.dat.ssu)))) {\n", " cat(\"Found NA values in SSU covariates. Replacing with 0...\\n\")\n", " na_count <- sum(is.na(values(cov.dat.ssu)))\n", " total_count <- ncell(cov.dat.ssu) * nlyr(cov.dat.ssu)\n", " cat(sprintf(\"Replacing %d NA values (%.2f%% of data)\\n\", \n", " na_count, (na_count/total_count)*100))\n", " \n", " cov.dat.ssu[is.na(cov.dat.ssu)] <- 0\n", " \n", " cat(\"NA values replaced with 0\\n\\n\")\n", "} else {\n", " cat(\"No NA values found in SSU covariates\\n\\n\")\n", "}\n", "\n", "names(cov.dat.ssu)\n", "\n", "writeRaster(cov.dat.ssu, file.path(paste0(results_dir,\"cov_dat_ssu_100m_clean.tif\")), overwrite=TRUE)\n", "# Load if needed. Similar process to cov.dat and land use vars.\n", "cov.dat.ssu <- rast(file.path(paste0(results_dir,\"cov_dat_ssu_100m_clean.tif\"))) # don't forget to check CRS\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 15 - GENERATE SSUs AND TSUs FOR TARGET PSUs" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 15 - GENERATE SSUs AND TSUs FOR TARGET PSUs =================================\n", "# Purpose: Within each target PSU, create SSUs and TSUs\n", "# This is the CORE of the three-stage sampling design\n", "\n", "# Initialize storage\n", "all_psus_tsus <- list()\n", "selected_ssus <- list()\n", "skipped_psus <- c()\n", "\n", "# MAIN LOOP: Process each target PSU\n", "# MAIN LOOP: Process each target PSU to generate SSUs and TSUs\n", "for (psu_id in 1:nrow(target.PSUs)) {\n", " \n", " # Select current PSU\n", " selected_psu <- target.PSUs[psu_id, ]\n", " \n", " # STEP 1: Generate 100m×100m SSU grid within the PSU boundary\n", " ssu_grid <- st_make_grid(selected_psu, cellsize = c(ssu_size, ssu_size), square = TRUE)\n", " ssu_grid_sf <- st_sf(geometry = ssu_grid)\n", " \n", " # Clip SSU grid to exact PSU boundary (removes partial cells outside PSU)\n", " ssu_grid_sf <- suppressWarnings(st_intersection(ssu_grid_sf, st_geometry(selected_psu)))\n", " \n", " # Convert to terra vector format for raster extraction\n", " ssu_grid_vect <- vect(ssu_grid_sf)\n", " \n", " # STEP 2: Calculate crop percentage for each SSU\n", " # Extracts from 20m resolution crop raster and calculates % coverage\n", " extracted_values <- extract(crops, ssu_grid_vect, fun = function(x) {\n", " sum(x > 0, na.rm = TRUE) / length(x) * 100\n", " })\n", " \n", " # Debug output: show range of crop coverage in this PSU\n", " cat(sprintf(\"\\nPSU %d \\n lu values range: %.2f to %.2f\\n\", \n", " psu_id, min(extracted_values[,2]), max(extracted_values[,2])))\n", " \n", " # STEP 3: Verify extraction succeeded and assign crop values\n", " if (ncol(extracted_values) >= 2) {\n", " ssu_grid_sf$lu <- extracted_values[, 2]\n", " \n", " # STEP 4: Handle split geometries from st_intersection\n", " # st_intersection sometimes creates MULTIPOLYGON - we need to merge them back\n", " ssu_grid_sf$ssu_temp_id <- 1:nrow(ssu_grid_sf) # Track original SSU IDs\n", " ssu_grid_sf <- st_cast(ssu_grid_sf, \"POLYGON\") # Split any MULTIPOLYGONs\n", " \n", " # Merge split parts back to single polygons per SSU\n", " ssu_grid_sf <- ssu_grid_sf %>%\n", " group_by(ssu_temp_id, lu) %>%\n", " summarise(geometry = st_union(geometry), .groups = \"drop\") %>%\n", " select(-ssu_temp_id) # Remove temporary ID\n", " \n", " # Store count before filtering for reporting\n", " total_ssus_before <- nrow(ssu_grid_sf)\n", " \n", " # STEP 5: Filter SSUs by actual 20m crop pixel count\n", " # This is more accurate than percentage and prevents TSU generation failures\n", " cat(\"Checking crop pixel availability for TSU generation...\\n\")\n", " ssu_grid_sf$crop_pixel_count <- sapply(1:nrow(ssu_grid_sf), function(i) {\n", " ssu_geom <- ssu_grid_sf[i, ]\n", " ssu_vect <- vect(ssu_geom)\n", " \n", " # Crop the 20m raster to this SSU and count valid pixels\n", " ssu_crop <- try(crop(crops, ssu_vect, mask = TRUE), silent = TRUE)\n", " if (inherits(ssu_crop, \"try-error\") || is.null(ssu_crop)) return(0)\n", " \n", " crop_vals <- values(ssu_crop, mat = FALSE)\n", " sum(crop_vals > 0, na.rm = TRUE)\n", " })\n", " \n", " # Minimum pixels needed: number of TSUs + safety buffer\n", " min_crop_pixels <- number_TSUs + 5\n", " ssu_grid_sf <- ssu_grid_sf[ssu_grid_sf$crop_pixel_count >= min_crop_pixels, ]\n", " \n", " # Report how many SSUs were removed\n", " cat(sprintf(\"SSUs after crop pixel filter: %d (removed %d)\\n\",\n", " nrow(ssu_grid_sf), total_ssus_before - nrow(ssu_grid_sf)))\n", " \n", " } else {\n", " # Extraction failed - skip this PSU\n", " warning(paste(\"PSU\", psu_id, \"returned insufficient extracted values. Skipping.\"))\n", " skipped_psus <- c(skipped_psus, psu_id)\n", " next\n", " }\n", " \n", " # Progress indicator\n", " cat(sprintf(\"\\rProgress: %.2f%% (%d out of %d)\\n\", \n", " (psu_id / nrow(target.PSUs)) * 100, psu_id, nrow(target.PSUs)))\n", " flush.console()\n", " \n", " # STEP 6: Check if enough SSUs remain for clustering\n", " total_ssus <- nrow(ssu_grid_sf)\n", " min_required_ssus <- max(num_primary_ssus + num_alternative_ssus, 8) # Default: 4+4=8\n", " \n", " if (total_ssus < min_required_ssus) {\n", " warning(paste(\"PSU\", psu_id, \"has only\", total_ssus, \n", " \"usable SSUs (min required:\", min_required_ssus, \"). Skipping.\"))\n", " skipped_psus <- c(skipped_psus, psu_id)\n", " next\n", " }\n", " \n", " # STEP 7: Extract 100m environmental covariates for each SSU\n", " ssu_grid_vect_filtered <- vect(ssu_grid_sf)\n", " ssu_covariates <- terra::extract(cov.dat.ssu, ssu_grid_vect_filtered, df = TRUE)\n", " \n", " # Aggregate covariates if any SSUs were split (collapse to one row per SSU)\n", " ssu_covariates <- ssu_covariates %>%\n", " group_by(ID) %>%\n", " summarise(across(everything(), ~mean(.x, na.rm = TRUE)))\n", " \n", " # Verify SSU count matches covariate count (critical for alignment)\n", " if (nrow(ssu_covariates) != nrow(ssu_grid_sf)) {\n", " warning(sprintf(\"PSU %d: Mismatch in SSU and covariate rows (%d vs %d). Skipping.\", \n", " psu_id, nrow(ssu_grid_sf), nrow(ssu_covariates)))\n", " skipped_psus <- c(skipped_psus, psu_id)\n", " next\n", " }\n", " \n", " # STEP 8: Prepare covariates for k-means clustering\n", " # Combine spatial data with environmental covariates\n", " ssu_data <- cbind(ssu_grid_sf, ssu_covariates[, -1])\n", " ssu_data_values <- st_drop_geometry(ssu_data)\n", " \n", " # Separate categorical variables (don't scale these)\n", " exclude <- grep(\"^geomorph_|^lu$\", names(ssu_data_values), value = TRUE)\n", " to_scale <- ssu_data_values[, !names(ssu_data_values) %in% exclude]\n", " to_keep <- ssu_data_values[, names(ssu_data_values) %in% exclude, drop = FALSE]\n", " \n", " # Remove columns with all NA values\n", " to_scale <- to_scale[, colSums(!is.na(to_scale)) > 0, drop = FALSE]\n", " \n", " # Identify zero-variance columns (would cause scaling errors)\n", " zero_variance_cols <- sapply(to_scale, function(x) sd(x, na.rm = TRUE) == 0)\n", " zero_variance_cols[is.na(zero_variance_cols)] <- TRUE\n", " \n", " # Scale only non-zero-variance columns\n", " scaled_part <- to_scale\n", " if (any(!zero_variance_cols)) {\n", " scaled_part[, !zero_variance_cols] <- scale(to_scale[, !zero_variance_cols])\n", " }\n", " \n", " # STEP 9: Remove incomplete cases and keep ssu_data aligned\n", " mygrd_ssu <- cbind(to_keep, scaled_part)\n", " complete_rows <- complete.cases(mygrd_ssu)\n", " mygrd_ssu <- mygrd_ssu[complete_rows, ]\n", " ssu_data <- ssu_data[complete_rows, ] # CRITICAL: Keep aligned with mygrd_ssu!\n", " \n", " # Verify enough SSUs remain after removing incomplete cases\n", " if (nrow(mygrd_ssu) < 4) {\n", " warning(paste(\"PSU\", psu_id, \"has too few SSUs (\", nrow(mygrd_ssu), \") to form 4 clusters. Skipping.\"))\n", " skipped_psus <- c(skipped_psus, psu_id)\n", " next\n", " }\n", " \n", " # STEP 10: Cluster SSUs using k-means\n", " # Number of clusters = number of target SSUs needed\n", " optimal_k <- num_primary_ssus\n", " \n", " kmeans_result <- kmeans(mygrd_ssu[, -1, drop = FALSE], centers = optimal_k, iter.max = 10000, nstart = 10)\n", " ssu_data$cluster <- as.factor(kmeans_result$cluster)\n", " \n", " # STEP 11: Select SSUs closest to cluster centers\n", " # These become target and replacement SSUs\n", " D <- rdist(x1 = kmeans_result$centers, x2 = mygrd_ssu[, -1, drop = FALSE])\n", " target_units <- apply(D, 1, function(x) order(x)[1]) # 1st closest = target\n", " replacement_units <- apply(D, 1, function(x) order(x)[2]) # 2nd closest = replacement\n", " \n", " # Verify we got enough SSUs\n", " if (length(target_units) < num_primary_ssus || length(replacement_units) < num_alternative_ssus) {\n", " warning(sprintf(\"PSU %d did not yield\", num_primary_ssus ,\"target and\", num_alternative_ssus, \"replacement SSUs. Skipping.\", psu_id))\n", " skipped_psus <- c(skipped_psus, psu_id)\n", " next\n", " }\n", " \n", " # STEP 12: Extract selected SSUs and add metadata\n", " target_ssus <- ssu_data[target_units, ]\n", " replacement_ssus <- ssu_data[replacement_units, ]\n", " \n", " target_ssus$SSU_Type <- \"Target\"\n", " replacement_ssus$SSU_Type <- \"Replacement\"\n", " target_ssus$SSU_ID <- 1:nrow(target_ssus)\n", " replacement_ssus$SSU_ID <- (nrow(target_ssus) + 1):(2 * nrow(target_ssus))\n", " \n", " # Link each replacement SSU to its corresponding target SSU (same cluster)\n", " replacement_ssus$replacement_for <- sapply(replacement_ssus$cluster, function(cl) {\n", " matched <- target_ssus$SSU_ID[target_ssus$cluster == cl]\n", " if (length(matched) > 0) return(matched[1]) else return(NA)\n", " })\n", " \n", " target_ssus$replacement_for <- NA\n", " \n", " # Add PSU identifier\n", " psu_actual_id <- selected_psu$ID\n", " target_ssus$PSU_ID <- psu_actual_id\n", " replacement_ssus$PSU_ID <- psu_actual_id\n", " \n", " # Combine target and replacement SSUs\n", " combined_ssus <- rbind(target_ssus, replacement_ssus)\n", " \n", " # STEP 13: Final check - must have exactly expected number of SSUs\n", " if (nrow(combined_ssus) != (num_primary_ssus + num_alternative_ssus)) {\n", " warning(sprintf(\"PSU %d has %d SSUs instead of\", num_primary_ssus + num_alternative_ssus, \". Skipping.\", psu_id, nrow(combined_ssus)))\n", " skipped_psus <- c(skipped_psus, psu_id)\n", " next\n", " }\n", " \n", " # Store SSUs for this PSU\n", " selected_ssus[[psu_id]] <- combined_ssus\n", " \n", " # STEP 14: Generate TSU point samples within each SSU\n", " primary_tsus <- lapply(1:nrow(target_ssus), function(i) {\n", " generate_tsu_points_within_ssu(target_ssus[i, ], number_TSUs, target_ssus$SSU_ID[i], \"Target\", crops)\n", " })\n", " \n", " alternative_tsus <- lapply(1:nrow(replacement_ssus), function(i) {\n", " generate_tsu_points_within_ssu(replacement_ssus[i, ], number_TSUs, replacement_ssus$SSU_ID[i], \"Replacement\", crops)\n", " })\n", " \n", " # STEP 15: Verify all TSUs were generated successfully\n", " all_tsus <- c(primary_tsus, alternative_tsus)\n", " tsu_counts <- sapply(all_tsus, function(x) if(is.null(x)) 0 else nrow(x))\n", " \n", " # Check: each SSU must have exactly number_TSUs points\n", " if (any(tsu_counts != number_TSUs)) {\n", " warning(sprintf(\"PSU %d: TSU generation failed for some SSUs. Expected %d TSUs per SSU, got: %s. Skipping.\", \n", " psu_id, number_TSUs, paste(tsu_counts, collapse=\",\")))\n", " skipped_psus <- c(skipped_psus, psu_id)\n", " # Remove this PSU from selected SSUs (failed at final step)\n", " selected_ssus[[psu_id]] <- NULL\n", " next\n", " }\n", " \n", " # Store all TSUs for this PSU\n", " all_psus_tsus[[psu_id]] <- do.call(rbind, Filter(Negate(is.null), all_tsus))\n", "}\n", "\n", "# Summary\n", "cat(sprintf(\"Successfully processed: %d PSUs\\n\", length(selected_ssus)))\n", "cat(sprintf(\"Skipped: %d PSUs\\n\", length(skipped_psus)))\n", "if (length(skipped_psus) > 0) {\n", " cat(sprintf(\"Skipped PSU IDs: %s\\n\", paste(skipped_psus, collapse = \", \")))\n", "}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 16 - COMBINE AND STRUCTURE TSU DATA" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 16 - COMBINE AND STRUCTURE TSU DATA =========================================\n", "# Purpose: Merge all TSUs into single spatial dataset with proper structure\n", "\n", "all_ssus <- do.call(rbind, selected_ssus)\n", "all_ssus <- all_ssus %>% mutate_at(vars(PSU_ID, SSU_ID), as.numeric)\n", "\n", "all_tsus <- do.call(rbind, all_psus_tsus)\n", "all_tsus <- all_tsus %>% mutate_at(vars(PSU_ID, SSU_ID), as.numeric)\n", "\n", "# Join TSU and SSU metadata\n", "all_tsus <- st_join(all_tsus, all_ssus[c(\"PSU_ID\", \"SSU_ID\", \"SSU_Type\", \"replacement_for\")])\n", "\n", "all_tsus <- all_tsus %>%\n", " select(PSU_ID = PSU_ID.x, SSU_ID = SSU_ID.y, SSU_Type = SSU_Type.y,\n", " Replacement_for = replacement_for, TSU_ID, geometry)\n", "\n", "# Label TSU types (Target = primary sample, Alternative = backup)\n", "all_tsus <- all_tsus %>%\n", " group_by(PSU_ID) %>%\n", " filter(n() == (num_primary_ssus + num_alternative_ssus) * number_TSUs) %>% # 24 for default\n", " group_by(PSU_ID, SSU_ID) %>%\n", " mutate(TSU_Type = ifelse(row_number() == 1, \"Target\", \"Alternative\")) %>%\n", " ungroup()\n", "\n", "all_tsus$PSU_Type <- \"Target\"\n", "\n", "all_tsus <- all_tsus %>%\n", " dplyr::select(\"PSU_ID\", \"SSU_ID\", \"SSU_Type\", \"Replacement_for\", \n", " \"TSU_ID\", \"TSU_Type\", \"geometry\")\n", "\n", "# Count target samples\n", "n_target_tsus <- all_tsus %>%\n", " filter(SSU_Type == \"Target\" & TSU_Type == \"Target\") %>%\n", " distinct(PSU_ID, SSU_ID) %>%\n", " nrow()\n", "\n", "cat(sprintf(\"\\nTotal target sampling locations: %d\\n\", n_target_tsus))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 17 - VISUALIZE SAMPLING DESIGN" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 17 - VISUALIZE SAMPLING DESIGN ==============================================\n", "# Purpose: Create map showing PSU, SSUs, and TSUs\n", "\n", "# Select one PSU for detailed visualization\n", "viz_psu_id <- selected_ssus[[1]]$PSU_ID[1]\n", "selected_psu_viz <- target.PSUs[target.PSUs$ID == viz_psu_id, ]\n", "\n", "bbox_psu <- st_bbox(selected_psu_viz)\n", "lu_bbox <- terra::crop(crops, selected_psu_viz, mask = TRUE)\n", "\n", "# Filter data for this PSU\n", "tsus_plot <- all_tsus[all_tsus$PSU_ID == viz_psu_id, ]\n", "ssus_plot <- all_ssus[all_ssus$PSU_ID == viz_psu_id, ]\n", "\n", "# Create detailed map\n", "labels <- c(\"Target PSU\", \"Target SSUs\", \"Replacement SSUs\", \"TSUs\")\n", "\n", "png(\n", " filename = file.path(paste0(results_dir,\"img/target_psu_ssu_tsu.png\")),\n", " width = 4000, # pixels\n", " height = 2800, # pixels\n", " res = 400 # DPI\n", ")\n", "\n", "ggplot() +\n", " geom_raster(data = as.data.frame(lu_bbox, xy = TRUE), aes(x = x, y = y, fill = lu)) +\n", " guides(fill = \"none\") +\n", " geom_sf(data = selected_psu_viz, fill = NA, \n", " aes(color = labels[1]), lwd = 0.8, show.legend = TRUE) +\n", " geom_sf(data = ssus_plot[ssus_plot$SSU_Type == \"Target\", ], fill = NA, \n", " aes(color = labels[2]), lwd = 0.6, show.legend = TRUE) +\n", " geom_sf(data = ssus_plot[ssus_plot$SSU_Type == \"Replacement\", ], fill = NA, \n", " aes(color = labels[3]), lwd = 0.6, show.legend = TRUE) +\n", " geom_sf(data = tsus_plot, aes(geometry = geometry, color = labels[4]), \n", " size = 0.5, shape = 19, show.legend = TRUE) +\n", " coord_sf(xlim = c(bbox_psu[\"xmin\"], bbox_psu[\"xmax\"]),\n", " ylim = c(bbox_psu[\"ymin\"], bbox_psu[\"ymax\"])) +\n", " labs(title = \"Example: three-stage sampling design\",\n", " subtitle = sprintf(\"PSU %d\", viz_psu_id),\n", " x = \"Longitude\", y = \"Latitude\", color = \"Legend\") +\n", " scale_color_manual(values = c(\"Target PSU\" = \"blue\",\n", " \"Target SSUs\" = \"blue\",\n", " \"Replacement SSUs\" = \"red\",\n", " \"TSUs\" = \"black\")) +\n", " theme_minimal()\n", "\n", "dev.off()\n", "\n", "ggsave(file.path(paste0(results_dir,\"img/sampling_design_example.png\")), \n", " width = 10, height = 10, dpi = 300)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 18 - EXPORT TARGET SAMPLING UNITS" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 18 - EXPORT TARGET SAMPLING UNITS ===========================================\n", "# Purpose: Save all shapefiles for field work\n", "\n", "# Add cluster information\n", "dfr <- PSU.df[,c(\"x\",\"y\",\"cluster\")]\n", "dfr$cluster <- as.numeric(dfr$cluster)\n", "dfr <- rasterFromXYZ(dfr)\n", "crs(dfr) <- epsg\n", "\n", "valid.PSU_clusters <- target.PSUs %>% \n", " mutate(cluster = extract(dfr, target.PSUs, fun = mean, na.rm = TRUE))\n", "\n", "all.PSU_clusters <- psu_grid %>% \n", " mutate(cluster = extract(dfr, psu_grid, fun = mean, na.rm = TRUE))\n", "\n", "all.PSU_clusters <- na.omit(all.PSU_clusters)\n", "\n", "valid.PSU_clusters <- valid.PSU_clusters %>% rename(Replace_ID = cluster)\n", "\n", "# Join cluster info to TSUs\n", "all_tsus <- st_join(all_tsus, valid.PSU_clusters)\n", "\n", "# Add sampling order\n", "all_tsus <- all_tsus %>%\n", " group_by(PSU_ID) %>%\n", " mutate(order = match(SSU_ID, unique(SSU_ID))) %>%\n", " ungroup()\n", "\n", "# Create unique site IDs\n", "# Format: [COUNTRY CODE][PSU ID]-[SSU ID]-[TSU ID][LAND USE]\n", "# Example: KANSAS0001-1-1C (PSU 1, SSU 1, TSU 1, Cropland)\n", "all_tsus$site_id <- paste0(ISO.code, sprintf(\"%04d\", all_tsus$PSU_ID), \n", " \"-\", all_tsus$SSU_ID, \n", " \"-\", all_tsus$TSU_ID, \"C\") # C for Cropland\n", "\n", "# Filter valid PSUs (those with complete TSU sets)\n", "psus_with_tsus <- unique(all_tsus$PSU_ID)\n", "valid.PSU_clusters_filtered <- valid.PSU_clusters %>%\n", " filter(ID %in% psus_with_tsus)\n", "\n", "# Export shapefiles\n", "write_sf(valid.PSU_clusters_filtered, \n", " paste0(results.path,\"PSUs_target.shp\"), overwrite=TRUE)\n", "write_sf(all_tsus, \n", " paste0(results.path,\"/TSUs_target.shp\"), overwrite=TRUE)\n", "write_sf(all.PSU_clusters, \n", " paste0(results.path,\"/PSU_pattern_cl.shp\"), overwrite=TRUE)\n", "writeRaster(dfr, paste0(results.path,\"/clusters.tif\"), overwrite=TRUE)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 19 - CALCULATE REPLACEMENT PSUs" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 19 - CALCULATE REPLACEMENT PSUs ============================================\n", "# Purpose: For each target PSU, find a similar replacement in same cluster\n", "# Why: Field teams need backups if target PSU becomes inaccessible\n", "\n", "# Find PSUs NOT selected as targets\n", "remaining.PSU_clusters <- all.PSU_clusters %>%\n", " filter(!(ID %in% valid.PSU_clusters$ID))\n", "\n", "# Get unique cluster IDs from targets\n", "unique_cluster <- distinct(valid.PSU_clusters, Replace_ID)$Replace_ID\n", "\n", "# Sample one replacement per cluster\n", "sampled_indices <- integer(0)\n", "\n", "for (clust in unique_cluster) {\n", " candidates_indices <- which(remaining.PSU_clusters$cluster == clust)\n", " \n", " if (length(candidates_indices) > 0) {\n", " sampled_index <- sample(candidates_indices, size = 1)\n", " sampled_indices <- c(sampled_indices, sampled_index)\n", " }\n", "}\n", "\n", "replacements <- remaining.PSU_clusters[sampled_indices, ]\n", "\n", "cat(sprintf(\"Selected %d replacement PSUs\\n\", nrow(replacements)))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 20 - GENERATE SSUs AND TSUs FOR REPLACEMENT PSUs" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 20 - GENERATE SSUs AND TSUs FOR REPLACEMENT PSUs ============================\n", "# Purpose: Create same sampling structure in replacement PSUs\n", "# Note: Code is identical to target PSU loop (Section 15)\n", "\n", "alt_psus_tsus_sf <- list()\n", "selected_ssus_sf <- list()\n", "skipped_psus_sf <- c()\n", "\n", "# DUPLICATE OF MAIN LOOP FOR REPLACEMENTS\n", "# (See Section 15 for detailed comments)\n", "for (psu_id in 1:nrow(replacements)) {\n", " selected_psu <- replacements[psu_id, ]\n", " \n", " # Generate SSUs within the selected PSU\n", " ssu_grid <- st_make_grid(selected_psu, cellsize = c(ssu_size, ssu_size), square = TRUE)\n", " ssu_grid_sf <- st_sf(geometry = ssu_grid)\n", " ssu_grid_sf <- suppressWarnings(st_intersection(ssu_grid_sf, st_geometry(selected_psu)))\n", " ssu_grid_vect <- vect(ssu_grid_sf)\n", " \n", " # Extract land use values (LU)\n", " extracted_values <- extract(crops, ssu_grid_vect, fun = function(x) {\n", " sum(x > 0, na.rm = TRUE) / length(x) * 100\n", " })\n", " \n", " # ADD THIS DEBUG LINE:\n", " cat(sprintf(\"\\nPSU %d\\n lu values range: %.2f to %.2f\\n\", \n", " psu_id, min(extracted_values[,2]), max(extracted_values[,2])))\n", " \n", " # Defensive check: ensure extracted_values has at least 2 columns\n", " if (ncol(extracted_values) >= 2) {\n", " ssu_grid_sf$lu <- extracted_values[, 2]\n", " \n", " # FIX: Handle MULTIPOLYGON geometries that cause row mismatches\n", " ssu_grid_sf$ssu_temp_id <- 1:nrow(ssu_grid_sf)\n", " ssu_grid_sf <- st_cast(ssu_grid_sf, \"POLYGON\")\n", " \n", " # Group split geometries back together\n", " ssu_grid_sf <- ssu_grid_sf %>%\n", " group_by(ssu_temp_id, lu) %>%\n", " summarise(geometry = st_union(geometry), .groups = \"drop\") %>%\n", " select(-ssu_temp_id)\n", " \n", " # Store count before filtering\n", " total_ssus_before <- nrow(ssu_grid_sf)\n", " \n", " # ONLY USE THIS FILTER: Direct pixel count (accurate!)\n", " cat(\"Checking crop pixel availability for TSU generation...\\n\")\n", " ssu_grid_sf$crop_pixel_count <- sapply(1:nrow(ssu_grid_sf), function(i) {\n", " ssu_geom <- ssu_grid_sf[i, ]\n", " ssu_vect <- vect(ssu_geom)\n", " ssu_crop <- try(crop(crops, ssu_vect, mask = TRUE), silent = TRUE)\n", " if (inherits(ssu_crop, \"try-error\") || is.null(ssu_crop)) return(0)\n", " crop_vals <- values(ssu_crop, mat = FALSE)\n", " sum(crop_vals > 0, na.rm = TRUE)\n", " })\n", " \n", " min_crop_pixels <- number_TSUs + 5\n", " ssu_grid_sf <- ssu_grid_sf[ssu_grid_sf$crop_pixel_count >= min_crop_pixels, ]\n", " \n", " cat(sprintf(\"SSUs after crop pixel filter: %d (removed %d)\\n\",\n", " nrow(ssu_grid_sf), total_ssus_before - nrow(ssu_grid_sf)))\n", " \n", " } else {\n", " warning(paste(\"PSU\", psu_id, \"returned insufficient extracted values. Skipping.\"))\n", " skipped_psus_sf <- c(skipped_psus_sf, psu_id)\n", " next\n", " }\n", " \n", " cat(sprintf(\"\\rProgress: %.2f%% (%d out of %d)\\n\", \n", " (psu_id / nrow(replacements)) * 100, psu_id, nrow(replacements)))\n", " flush.console()\n", " \n", " # Count available SSUs\n", " total_ssus <- nrow(ssu_grid_sf)\n", " min_required_ssus <- max(num_primary_ssus + num_alternative_ssus, 8) # 4 clusters x 2 SSUs each\n", " \n", " if (total_ssus < min_required_ssus) {\n", " warning(paste(\"PSU\", psu_id, \"has only\", total_ssus, \n", " \"usable SSUs (min required:\", min_required_ssus, \"). Skipping.\"))\n", " skipped_psus_sf <- c(skipped_psus_sf, psu_id)\n", " next\n", " }\n", " \n", " # Extract covariates\n", " # Convert filtered SSU grid to vector and extract covariates *after filtering*\n", " ssu_grid_vect_filtered <- vect(ssu_grid_sf)\n", " ssu_covariates <- terra::extract(cov.dat.ssu, ssu_grid_vect_filtered, df = TRUE)\n", " \n", " ssu_covariates <- ssu_covariates %>%\n", " group_by(ID) %>%\n", " summarise(across(everything(), ~mean(.x, na.rm = TRUE)))\n", " \n", " # Ensure alignment by checking row counts before binding\n", " if (nrow(ssu_covariates) != nrow(ssu_grid_sf)) {\n", " warning(sprintf(\"PSU %d: Mismatch in SSU and covariate rows (%d vs %d). Skipping.\", \n", " psu_id, nrow(ssu_grid_sf), nrow(ssu_covariates)))\n", " skipped_psus_sf <- c(skipped_psus_sf, psu_id)\n", " next\n", " }\n", " \n", " ssu_data <- cbind(ssu_grid_sf, ssu_covariates[, -1])\n", " ssu_data_values <- st_drop_geometry(ssu_data)\n", " \n", " exclude <- grep(\"^geomorph_|^lu$\", names(ssu_data_values), value = TRUE)\n", " to_scale <- ssu_data_values[, !names(ssu_data_values) %in% exclude]\n", " to_keep <- ssu_data_values[, names(ssu_data_values) %in% exclude, drop = FALSE]\n", " \n", " to_scale <- to_scale[, colSums(!is.na(to_scale)) > 0, drop = FALSE]\n", " zero_variance_cols <- sapply(to_scale, function(x) sd(x, na.rm = TRUE) == 0)\n", " zero_variance_cols[is.na(zero_variance_cols)] <- TRUE\n", " \n", " scaled_part <- to_scale\n", " if (any(!zero_variance_cols)) {\n", " scaled_part[, !zero_variance_cols] <- scale(to_scale[, !zero_variance_cols])\n", " }\n", " \n", " mygrd_ssu <- cbind(to_keep, scaled_part)\n", " complete_rows <- complete.cases(mygrd_ssu)\n", " mygrd_ssu <- mygrd_ssu[complete_rows, ]\n", " # Also filter ssu_data to match\n", " ssu_data <- ssu_data[complete_rows, ]\n", " \n", " if (nrow(mygrd_ssu) < 4) {\n", " warning(paste(\"PSU\", psu_id, \"has too few SSUs (\", nrow(mygrd_ssu), \") to form 4 clusters. Skipping.\"))\n", " skipped_psus_sf <- c(skipped_psus_sf, psu_id)\n", " next\n", " }\n", " \n", " # Fixed number of clusters\n", " optimal_k <- num_primary_ssus\n", " \n", " kmeans_result <- kmeans(mygrd_ssu[, -1, drop = FALSE], centers = optimal_k, iter.max = 10000, nstart = 10)\n", " ssu_data$cluster <- as.factor(kmeans_result$cluster)\n", " \n", " # Compute distances and pick SSUs closest to centers\n", " D <- rdist(x1 = kmeans_result$centers, x2 = mygrd_ssu[, -1, drop = FALSE])\n", " target_units <- apply(D, 1, function(x) order(x)[1])\n", " replacement_units <- apply(D, 1, function(x) order(x)[2])\n", " \n", " if (length(target_units) < num_primary_ssus || length(replacement_units) < num_alternative_ssus) {\n", " warning(sprintf(\"PSU %d did not yield\", num_primary_ssus ,\"target and\", num_alternative_ssus, \"replacement SSUs. Skipping.\", psu_id))\n", " skipped_psus_sf <- c(skipped_psus_sf, psu_id)\n", " next\n", " }\n", " \n", " target_ssus <- ssu_data[target_units, ]\n", " replacement_ssus <- ssu_data[replacement_units, ]\n", " \n", " target_ssus$SSU_Type <- \"Target\"\n", " replacement_ssus$SSU_Type <- \"Replacement\"\n", " target_ssus$SSU_ID <- 1:nrow(target_ssus)\n", " replacement_ssus$SSU_ID <- (nrow(target_ssus) + 1):(2 * nrow(target_ssus))\n", " \n", " replacement_ssus$replacement_for <- sapply(replacement_ssus$cluster, function(cl) {\n", " matched <- target_ssus$SSU_ID[target_ssus$cluster == cl]\n", " if (length(matched) > 0) return(matched[1]) else return(NA)\n", " })\n", " \n", " target_ssus$replacement_for <- NA\n", " \n", " # Add PSU_ID and combine\n", " psu_actual_id <- selected_psu$ID\n", " target_ssus$PSU_ID <- psu_actual_id\n", " replacement_ssus$PSU_ID <- psu_actual_id\n", " # Combine SSUs\n", " combined_ssus <- rbind(target_ssus, replacement_ssus)\n", " \n", " # CHECK: Must have exactly 8 SSUs (4 target + 4 replacement)\n", " if (nrow(combined_ssus) != (num_primary_ssus + num_alternative_ssus)) {\n", " warning(sprintf(\"PSU %d has %d SSUs instead of\", num_primary_ssus + num_alternative_ssus, \". Skipping.\", psu_id, nrow(combined_ssus)))\n", " skipped_psus_sf <- c(skipped_psus_sf, psu_id)\n", " next\n", " }\n", " \n", " # Store only if we have exactly 8 SSUs\n", " selected_ssus_sf[[psu_id]] <- combined_ssus\n", " \n", " ### Generate TSUs ###\n", " primary_tsus <- lapply(1:nrow(target_ssus), function(i) {\n", " generate_tsu_points_within_ssu(target_ssus[i, ], number_TSUs, target_ssus$SSU_ID[i], \"Target\", crops)\n", " })\n", " \n", " alternative_tsus <- lapply(1:nrow(replacement_ssus), function(i) {\n", " generate_tsu_points_within_ssu(replacement_ssus[i, ], number_TSUs, replacement_ssus$SSU_ID[i], \"Replacement\", crops)\n", " })\n", " \n", " # CHECK: Verify all TSUs were generated successfully\n", " all_tsus <- c(primary_tsus, alternative_tsus)\n", " tsu_counts <- sapply(all_tsus, function(x) if(is.null(x)) 0 else nrow(x))\n", " \n", " # Each SSU should have exactly number_TSUs (3) TSUs\n", " if (any(tsu_counts != number_TSUs)) {\n", " warning(sprintf(\"PSU %d: TSU generation failed for some SSUs. Expected %d TSUs per SSU, got: %s. Skipping.\", \n", " psu_id, number_TSUs, paste(tsu_counts, collapse=\",\")))\n", " skipped_psus_sf <- c(skipped_psus_sf, psu_id)\n", " # Remove this PSU from selected_ssus\n", " selected_ssus_sf[[psu_id]] <- NULL\n", " next\n", " }\n", " \n", " alt_psus_tsus_sf[[psu_id]] <- do.call(rbind, Filter(Negate(is.null), all_tsus))\n", "}\n", "\n", "cat(sprintf(\"Replacement PSUs: %d successful, %d skipped\\n\", \n", " length(selected_ssus_sf), length(skipped_psus_sf)))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 21 - STRUCTURE REPLACEMENT TSU DATA" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 21 - STRUCTURE REPLACEMENT TSU DATA =========================================\n", "\n", "all_ssus_combined_sf <- do.call(rbind, selected_ssus_sf)\n", "all_ssus_combined_sf <- all_ssus_combined_sf %>%\n", " mutate_at(vars(PSU_ID, SSU_ID), as.numeric)\n", "\n", "alt_tsus_combined_sf <- do.call(rbind, alt_psus_tsus_sf)\n", "alt_tsus_combined_sf <- alt_tsus_combined_sf %>%\n", " mutate_at(vars(PSU_ID, SSU_ID), as.numeric)\n", "\n", "alt_tsus_combined_sf <- st_join(alt_tsus_combined_sf, \n", " all_ssus_combined_sf[c(\"PSU_ID\", \"SSU_ID\", \"SSU_Type\", \"replacement_for\")])\n", "\n", "alt_tsus_combined_sf <- alt_tsus_combined_sf %>%\n", " select(PSU_ID = PSU_ID.x, SSU_ID = SSU_ID.y, SSU_Type = SSU_Type.y,\n", " Replacement_for = replacement_for, TSU_ID, geometry)\n", "\n", "alt_tsus_combined_sf <- alt_tsus_combined_sf %>%\n", " group_by(PSU_ID) %>%\n", " filter(n() == (num_primary_ssus + num_alternative_ssus) * number_TSUs) %>%\n", " group_by(PSU_ID, SSU_ID) %>%\n", " mutate(TSU_Type = ifelse(row_number() == 1, \"Target\", \"Alternative\")) %>%\n", " ungroup()\n", "\n", "alt_tsus_combined_sf$PSU_Type <- \"Replacement\"\n", "\n", "alt_tsus_combined_sf <- alt_tsus_combined_sf %>%\n", " dplyr::select(\"PSU_ID\", \"SSU_ID\", \"SSU_Type\", \"Replacement_for\", \n", " \"TSU_ID\", \"TSU_Type\", \"geometry\")\n", "\n", "n_replacement_tsus <- alt_tsus_combined_sf %>%\n", " filter(SSU_Type == \"Target\" & TSU_Type == \"Target\") %>%\n", " distinct(PSU_ID, SSU_ID) %>%\n", " nrow()\n", "\n", "cat(sprintf(\"Total replacement sampling locations: %d\\n\", n_replacement_tsus))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 22 - EXPORT REPLACEMENT SAMPLING UNITS" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 22 - EXPORT REPLACEMENT SAMPLING UNITS ======================================\n", "\n", "replacements <- replacements %>% rename(Replace_ID = cluster)\n", "\n", "alt_tsus_combined_sf <- st_join(alt_tsus_combined_sf, replacements)\n", "\n", "alt_tsus_combined_sf <- alt_tsus_combined_sf %>%\n", " group_by(PSU_ID) %>%\n", " mutate(order = match(SSU_ID, unique(SSU_ID))) %>%\n", " ungroup()\n", "\n", "# Create site IDs\n", "alt_tsus_combined_sf$site_id <- paste0(ISO.code, sprintf(\"%04d\", alt_tsus_combined_sf$PSU_ID), \n", " \"-\", alt_tsus_combined_sf$SSU_ID, \n", " \"-\", alt_tsus_combined_sf$TSU_ID, \"C\")\n", "\n", "# Filter and export\n", "psus_with_tsus_sf <- unique(alt_tsus_combined_sf$PSU_ID)\n", "replacements_filtered <- replacements %>%\n", " filter(ID %in% psus_with_tsus_sf)\n", "\n", "write_sf(replacements_filtered, \n", " paste0(results.path,\"/PSUs_replacements.shp\"), overwrite=TRUE)\n", "write_sf(alt_tsus_combined_sf, \n", " paste0(results.path,\"/TSUs_replacements.shp\"), overwrite=TRUE)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 23 - AVAILABILITY ANALYSIS" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 23 - AVAILABILITY ANALYSIS ==================================================\n", "# Purpose: Count how many backup PSUs are available per cluster\n", "# Useful for risk assessment and field planning\n", "\n", "valid_counts <- valid.PSU_clusters %>%\n", " group_by(Replace_ID) %>%\n", " summarise(Count = n())\n", "\n", "remaining_counts <- remaining.PSU_clusters %>%\n", " group_by(cluster) %>%\n", " summarise(Count = n())\n", "\n", "availability <- st_join(valid_counts, remaining_counts, \n", " by = \"cluster\", suffix = c(\"_valid\", \"_remaining\"))\n", "\n", "write_sf(availability, paste0(results.path,\"/availability.shp\"), overwrite=TRUE)\n", "\n", "################################################################################\n", "## PART 2: MERGE MULTIPLE LAND USES AND CREATE UNIFIED IDs\n", "################################################################################\n", "# Purpose: Combine cropland, grassland, and forest sampling designs\n", "# Creates country-wide unique identifiers for all sites\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 24 - SETUP FOR MERGING" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 24 - SETUP FOR MERGING ======================================================\n", "\n", "# Define folder for merged outputs\n", "folder_all <- paste0(results_dir, \"all/\")\n", "if (!file.exists(folder_all)){\n", " dir.create(folder_all)\n", "}\n", "\n", "# Load country boundaries (for provincial statistics)\n", "country_boundaries <- sf::st_read(paste0(shp.path,\"roi_kansas_adm2_us_epsg_4326.shp\"), quiet=TRUE)\n", "country_boundaries$country <- ISO.code\n", "head(country_boundaries, 5)\n", "country_boundaries$province <- country_boundaries$NAM_2\n", "\n", "if(!same.crs(country_boundaries, epsg)){\n", " country_boundaries <- country_boundaries %>%\n", " st_as_sf() %>% sf::st_transform(crs=epsg)\n", "}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 25 - IMPORT ALL LAND USE SHAPEFILES" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 25 - IMPORT ALL LAND USE SHAPEFILES =========================================\n", "# Purpose: Load PSUs and TSUs from all three land use types\n", "\n", "# Initialize empty lists\n", "psus_target_list <- list()\n", "tsus_target_list <- list()\n", "psus_repl_list <- list()\n", "tsus_repl_list <- list()\n", "\n", "# Define land use types to check\n", "landuse_types <- c(\"cropland\", \"grassland\", \"forest\")\n", "landuse_codes <- c(\"C\", \"G\", \"F\")\n", "\n", "# TARGET PSUs - Load only if files exist\n", "for (i in seq_along(landuse_types)) {\n", " file_path <- paste0(results_dir, landuse_types[i], \"/PSUs_target.shp\")\n", " \n", " if (file.exists(file_path)) {\n", " temp_psu <- sf::st_read(file_path, quiet = TRUE)\n", " temp_psu$lulc <- landuse_codes[i]\n", " psus_target_list[[landuse_types[i]]] <- temp_psu\n", " cat(sprintf(\" ✓ Loaded %s (%d PSUs)\\n\", landuse_types[i], nrow(temp_psu)))\n", " } else {\n", " cat(sprintf(\" ✗ Skipped %s (file not found)\\n\", landuse_types[i]))\n", " }\n", "}\n", "\n", "# TARGET TSUs - Load only if files exist\n", "for (i in seq_along(landuse_types)) {\n", " file_path <- paste0(results_dir, landuse_types[i], \"/TSUs_target.shp\")\n", " \n", " if (file.exists(file_path)) {\n", " temp_tsu <- sf::st_read(file_path, quiet = TRUE)\n", " temp_tsu$lulc <- landuse_codes[i]\n", " tsus_target_list[[landuse_types[i]]] <- temp_tsu\n", " cat(sprintf(\" ✓ Loaded %s (%d TSUs)\\n\", landuse_types[i], nrow(temp_tsu)))\n", " } else {\n", " cat(sprintf(\" ✗ Skipped %s (file not found)\\n\", landuse_types[i]))\n", " }\n", "}\n", "\n", "# REPLACEMENT PSUs - Load only if files exist\n", "for (i in seq_along(landuse_types)) {\n", " file_path <- paste0(results_dir, landuse_types[i], \"/PSUs_replacements.shp\")\n", " \n", " if (file.exists(file_path)) {\n", " temp_psu <- sf::st_read(file_path, quiet = TRUE)\n", " temp_psu$lulc <- landuse_codes[i]\n", " psus_repl_list[[landuse_types[i]]] <- temp_psu\n", " cat(sprintf(\" ✓ Loaded %s (%d PSUs)\\n\", landuse_types[i], nrow(temp_psu)))\n", " } else {\n", " cat(sprintf(\" ✗ Skipped %s (file not found)\\n\", landuse_types[i]))\n", " }\n", "}\n", "\n", "# REPLACEMENT TSUs - Load only if files exist\n", "for (i in seq_along(landuse_types)) {\n", " file_path <- paste0(results_dir, landuse_types[i], \"/TSUs_replacements.shp\")\n", " \n", " if (file.exists(file_path)) {\n", " temp_tsu <- sf::st_read(file_path, quiet = TRUE)\n", " temp_tsu$lulc <- landuse_codes[i]\n", " tsus_repl_list[[landuse_types[i]]] <- temp_tsu\n", " cat(sprintf(\" ✓ Loaded %s (%d TSUs)\\n\", landuse_types[i], nrow(temp_tsu)))\n", " } else {\n", " cat(sprintf(\" ✗ Skipped %s (file not found)\\n\", landuse_types[i]))\n", " }\n", "}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 26 - MERGE LAND USE TYPES" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 26 - MERGE LAND USE TYPES ===================================================\n", "# Purpose: Combine all available land uses into single datasets\n", "\n", "# Merge TARGET PSUs\n", "psus_target <- sf::st_as_sf(data.table::rbindlist(psus_target_list))\n", "cat(sprintf(\"✓ Merged %d target PSUs from %d land use type(s)\\n\", \n", " nrow(psus_target), length(psus_target_list)))\n", "\n", "# Merge TARGET TSUs\n", "tsus_target <- sf::st_as_sf(data.table::rbindlist(tsus_target_list))\n", "cat(sprintf(\"✓ Merged %d target TSUs from %d land use type(s)\\n\", \n", " nrow(tsus_target), length(tsus_target_list)))\n", "\n", "# Merge REPLACEMENT PSUs (if any exist)\n", "if (length(psus_repl_list) > 0) {\n", " psus_repl <- sf::st_as_sf(data.table::rbindlist(psus_repl_list))\n", " cat(sprintf(\"✓ Merged %d replacement PSUs from %d land use type(s)\\n\", \n", " nrow(psus_repl), length(psus_repl_list)))\n", "} else {\n", " warning(\"No replacement PSU files found - skipping replacement PSU merge\")\n", "}\n", "\n", "# Merge REPLACEMENT TSUs (if any exist)\n", "if (length(tsus_repl_list) > 0) {\n", " tsus_repl <- sf::st_as_sf(data.table::rbindlist(tsus_repl_list))\n", " cat(sprintf(\"✓ Merged %d replacement TSUs from %d land use type(s)\\n\", \n", " nrow(tsus_repl), length(tsus_repl_list)))\n", "} else {\n", " warning(\"No replacement TSU files found - skipping replacement TSU merge\")\n", "}\n", "\n", "# Summary by land use\n", "for (lu in unique(psus_target$lulc)) {\n", " lu_name <- switch(lu,\n", " \"C\" = \"Cropland\",\n", " \"G\" = \"Grassland\", \n", " \"F\" = \"Forest\",\n", " lu)\n", " n_psus <- sum(psus_target$lulc == lu)\n", " n_tsus <- sum(tsus_target$lulc == lu)\n", " cat(sprintf(\" %s: %d PSUs, %d TSUs\\n\", lu_name, n_psus, n_tsus))\n", "}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 27 - CREATE UNIQUE IDs FOR TARGET PSUs" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 27 - CREATE UNIQUE IDs FOR TARGET PSUs ======================================\n", "# Purpose: Assign sequential IDs across all land uses\n", "# ID Structure: Country-wide sequential number\n", "\n", "# Create composite ID (cluster + land use)\n", "psus_target$PSU_R_LULC_ID <- paste0(psus_target$Replace_ID, \"-\", psus_target$lulc)\n", "\n", "# Assign sequential country-wide IDs\n", "psus_target[[paste0(ISO.code,\"_PSU_ID\")]] <- 1:nrow(psus_target)\n", "\n", "psus_target <- psus_target %>%\n", " select(ID, Replace_ID, lulc, PSU_R_LULC_ID, \n", " all_of(paste0(ISO.code, \"_PSU_ID\")), everything())\n", "\n", "cat(sprintf(\"Assigned %d unique PSU IDs\\n\", nrow(psus_target)))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 28 - CREATE UNIQUE IDs FOR REPLACEMENT PSUs" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 28 - CREATE UNIQUE IDs FOR REPLACEMENT PSUs =================================\n", "\n", "psus_repl$PSU_R_LULC_ID <- paste0(psus_repl$Replace_ID, \"-\", psus_repl$lulc)\n", "\n", "# Continue numbering after target PSUs\n", "start_id <- nrow(psus_target) + 1\n", "psus_repl[[paste0(ISO.code,\"_PSU_ID\")]] <- start_id:(start_id + nrow(psus_repl) - 1)\n", "\n", "psus_repl <- psus_repl %>%\n", " select(ID, Replace_ID, lulc, PSU_R_LULC_ID, \n", " all_of(paste0(ISO.code, \"_PSU_ID\")), everything())\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 29 - LINK TARGETS AND REPLACEMENTS" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 29 - LINK TARGETS AND REPLACEMENTS ==========================================\n", "# Purpose: Each target PSU knows its replacement PSU ID\n", "\n", "# Match replacement IDs to targets\n", "index <- match(psus_target$PSU_R_LULC_ID, psus_repl$PSU_R_LULC_ID)\n", "psus_target[[\"PSU_R_ID\"]] <- psus_repl[[paste0(ISO.code,\"_PSU_ID\")]][index]\n", "\n", "psus_target <- psus_target %>%\n", " select(ID, Replace_ID, lulc, PSU_R_LULC_ID, PSU_R_ID, \n", " all_of(paste0(ISO.code, \"_PSU_ID\")), everything())\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 30 - ASSIGN UNIQUE IDs TO TARGET TSUs" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 30 - ASSIGN UNIQUE IDs TO TARGET TSUs =======================================\n", "# Purpose: Create site-level unique identifiers\n", "# ID Format: COUNTRY[####]-[#]-[#]L\n", "# Example: TUN0001-1-1C = Tunisia, PSU 1, SSU 1, TSU 1, Cropland\n", "\n", "# Create matching key\n", "psus_target$PSU_T_LULC_ID <- paste0(psus_target$ID, \"-\", psus_target$lulc)\n", "tsus_target$PSU_T_LULC_ID <- paste0(tsus_target$PSU_ID, \"-\", tsus_target$lulc)\n", "\n", "# Transfer PSU IDs to TSUs\n", "index <- match(tsus_target$PSU_T_LULC_ID, psus_target$PSU_T_LULC_ID)\n", "tsus_target[[paste0(ISO.code,\"_PSU_ID\")]] <- psus_target[[paste0(ISO.code,\"_PSU_ID\")]][index]\n", "tsus_target[[\"PSU_R_ID\"]] <- psus_target[[\"PSU_R_ID\"]][index]\n", "\n", "# Clean up column names\n", "tsus_target <- tsus_target %>%\n", " rename_with(~ str_replace_all(., c(\"Typ\" = \"Type\", \"^Rplcmn_\" = \"SSU_Repl\")))\n", "\n", "tsus_target$PSU_Type <- \"Target\"\n", "\n", "tsus_target <- tsus_target %>%\n", " select(all_of(paste0(ISO.code, \"_PSU_ID\")), PSU_Type, order, \n", " SSU_Type, SSU_Repl, TSU_ID, TSU_Type, site_id, lulc, PSU_R_ID)\n", "\n", "# Rename for clarity\n", "names(tsus_target) <- gsub(paste0(ISO.code, \"_PSU_ID\"), \"PSU_ID\", names(tsus_target))\n", "names(tsus_target) <- gsub(\"order\", \"SSU_ID\", names(tsus_target))\n", "\n", "# FINAL SITE ID GENERATION\n", "# This is the ID field teams will use in the field\n", "tsus_target$site_id <- paste0(ISO.code, sprintf(\"%04d\", tsus_target$PSU_ID), \n", " \"-\", tsus_target$SSU_ID, \n", " \"-\", tsus_target$TSU_ID, tsus_target$lulc)\n", "\n", "cat(sprintf(\"Created %d unique site IDs\\n\", nrow(tsus_target)))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 31 - ASSIGN UNIQUE IDs TO REPLACEMENT TSUs" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 31 - ASSIGN UNIQUE IDs TO REPLACEMENT TSUs ==================================\n", "\n", "psus_repl$PSU_T_LULC_ID <- paste0(psus_repl$ID, \"-\", psus_repl$lulc)\n", "tsus_repl$PSU_T_LULC_ID <- paste0(tsus_repl$PSU_ID, \"-\", tsus_repl$lulc)\n", "\n", "# Transfer IDs\n", "index <- match(tsus_repl$PSU_T_LULC_ID, psus_repl$PSU_T_LULC_ID)\n", "tsus_repl[[paste0(ISO.code,\"_PSU_ID\")]] <- psus_repl[[paste0(ISO.code,\"_PSU_ID\")]][index]\n", "\n", "# Link to target PSUs\n", "tsus_repl[[\"PSU_T_ID\"]] <- psus_target[[paste0(ISO.code, \"_PSU_ID\")]][\n", " match(tsus_repl[[paste0(ISO.code,\"_PSU_ID\")]], psus_target$PSU_R_ID)]\n", "\n", "tsus_repl <- tsus_repl %>%\n", " rename_with(~ str_replace_all(., c(\"Typ\" = \"Type\", \"^Rplcmn_\" = \"SSU_Repl\")))\n", "\n", "tsus_repl$PSU_Type <- \"Replacement\"\n", "\n", "tsus_repl <- tsus_repl %>%\n", " select(all_of(paste0(ISO.code, \"_PSU_ID\")), PSU_Type, order, \n", " SSU_Type, SSU_Repl, TSU_ID, TSU_Type, site_id, lulc, PSU_T_ID)\n", "\n", "names(tsus_repl) <- gsub(paste0(ISO.code, \"_PSU_ID\"), \"PSU_ID\", names(tsus_repl))\n", "names(tsus_repl) <- gsub(\"order\", \"SSU_ID\", names(tsus_repl))\n", "\n", "tsus_repl$site_id <- paste0(ISO.code, sprintf(\"%04d\", tsus_repl$PSU_ID), \n", " \"-\", tsus_repl$SSU_ID, \n", " \"-\", tsus_repl$TSU_ID, tsus_repl$lulc)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 32 - EXPORT MERGED SHAPEFILES" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 32 - EXPORT MERGED SHAPEFILES ===============================================\n", "# Purpose: Save final, field-ready shapefiles\n", "\n", "# TARGET PSUs (simplified)\n", "psus_target_export <- psus_target %>%\n", " select(PSU_ID = all_of(paste0(ISO.code, \"_PSU_ID\")), \n", " Replace_ID = PSU_R_ID, lulc)\n", "\n", "write_sf(psus_target_export, paste0(folder_all, \"all_psus_target.shp\"), overwrite = T)\n", "\n", "# REPLACEMENT PSUs\n", "psus_repl_export <- psus_repl %>%\n", " select(Replace_ID = all_of(paste0(ISO.code, \"_PSU_ID\")), lulc)\n", "\n", "write_sf(psus_repl_export, paste0(folder_all, \"all_psus_replacements.shp\"), overwrite = T)\n", "\n", "# TARGET TSUs\n", "names(tsus_target) <- gsub(\"PSU_R_ID\", \"Replace_ID\", names(tsus_target))\n", "\n", "write_sf(tsus_target, paste0(folder_all, \"all_tsus_target.shp\"), overwrite = T)\n", "\n", "# REPLACEMENT TSUs\n", "names(tsus_repl) <- gsub(\"PSU_T_ID\", \"PSU_Target_ID\", names(tsus_repl))\n", "\n", "write_sf(tsus_repl, paste0(folder_all, \"all_tsus_replacements.shp\"), overwrite = T)\n", "\n", "cat(\"\\n✓ All merged shapefiles exported to:\", folder_all, \"\\n\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 33 - GENERATE SUMMARY STATISTICS" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 33 - GENERATE SUMMARY STATISTICS ============================================\n", "# Purpose: Create distribution tables and graphs\n", "\n", "# Extract unique target sites only (1 per SSU)\n", "tsus_uniq_sites <- dplyr::filter(tsus_target, \n", " SSU_Type == \"Target\" & TSU_Type == \"Target\")\n", "\n", "# Join with provinces\n", "tsus_uniq_sites <- sf::st_join(tsus_uniq_sites, country_boundaries)\n", "tsus_uniq_sites <- tsus_uniq_sites %>%\n", " dplyr::select(c(site_id, country, province, lulc, geometry))\n", "\n", "# Country-level statistics\n", "sites_distribution <- as.data.frame(tsus_uniq_sites) %>%\n", " group_by(country, lulc) %>%\n", " summarise(Sites = n(), .groups = \"drop\")\n", "\n", "print(sites_distribution)\n", "\n", "# Visualize country distribution\n", "custom_colors <- RColorBrewer::brewer.pal(3, \"BrBG\")\n", "\n", "ggplot(sites_distribution, aes(x = lulc, y = Sites, fill = lulc)) +\n", " geom_bar(stat = \"identity\", width = 0.7) +\n", " geom_text(aes(label = Sites), vjust = -0.5, size = 4) +\n", " labs(title = \"Sampling Site Distribution by Land Use\",\n", " x = \"Land Use\", y = \"Number of Sites\", fill = \"Land Use\") +\n", " scale_x_discrete(labels = c(\"C\" = \"Cropland\", \"F\" = \"Forest\", \"G\" = \"Grassland\")) +\n", " scale_fill_manual(values = custom_colors, \n", " labels = c(\"C\" = \"Cropland\", \"F\" = \"Forest\", \"G\" = \"Grassland\")) +\n", " ylim(0, max(sites_distribution$Sites) * 1.2) +\n", " theme_minimal() +\n", " theme(plot.title = element_text(hjust = 0.5, size = 16, face = \"bold\"),\n", " axis.text = element_text(size = 12),\n", " legend.position = \"top\")\n", "\n", "ggsave(paste0(results_dir, \"img//final_site_distribution.png\"), \n", " width = 10, height = 6, dpi = 300)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 34 - PROVINCIAL STATISTICS" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "## 34 - PROVINCIAL STATISTICS ==================================================\n", "\n", "# Province-level statistics\n", "sites_distribution_prov <- as.data.frame(tsus_uniq_sites) %>%\n", " group_by(province, lulc) %>%\n", " summarise(Sites = n(), .groups = \"drop\")\n", "\n", "print(sites_distribution_prov)\n", "\n", "# Visualize provincial distribution\n", "ncolors <- length(unique(sites_distribution_prov$province))\n", "custom_colors_prov <- colorRampPalette(brewer.pal(11, \"Spectral\"))(ncolors)\n", "\n", "ggplot(sites_distribution_prov, aes(x = lulc, y = Sites, fill = province)) +\n", " geom_bar(stat = \"identity\", position = \"dodge\", width = 0.7) +\n", " geom_text(aes(label = Sites), \n", " position = position_dodge(width = 0.7), \n", " vjust = -0.5, size = 3) +\n", " labs(title = \"Site Distribution by Province and Land Use\",\n", " x = \"Land Use\", y = \"Number of Sites\", fill = \"Province\") +\n", " scale_x_discrete(labels = c(\"C\" = \"Cropland\", \"F\" = \"Forest\", \"G\" = \"Grassland\")) +\n", " scale_fill_manual(values = custom_colors_prov) +\n", " ylim(0, max(sites_distribution_prov$Sites) * 1.2) +\n", " theme_minimal() +\n", " theme(plot.title = element_text(hjust = 0.5, size = 16, face = \"bold\"),\n", " legend.position = \"top\")\n", "\n", "# The user can adjust the width and height for better visualisation\n", "ggsave(paste0(results_dir, \"img/final_site_distribution_province.png\"), \n", " width = 12, height = 8, dpi = 300)\n", "\n", "################################################################################\n", "## END\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" ] } ] }