Technical Guide & Modeling Protocol

This protocol describes the development and application of BhutanBioClims—a comprehensive, high-resolution (250 m) bioclimatic dataset for Bhutan derived from CMIP6 Global Climate Models (GCMs). The pipeline provides 121 sets of 19 standard BIOCLIM variables (Booth et al., 2014), spanning historical (1986–2015) and three future time slices (2021–2050, 2051–2080, 2071–2100) across four Shared Socioeconomic Pathways (SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) using 10 GCMs.

CMIP6 GCM outputs were downscaled against historical data developed with Bhutan's national weather station network (Stewart et al., 2017; Stewart et al., 2021) using the delta change method applied to anomalies interpolated via bivariate thin plate splines. Bioclimatic variables were computed in R using the terra package (Hijmans et al., 2023), with automatic unit conversion from Kelvin to °C and from kg m⁻² s⁻¹ to mm month⁻¹ where required. The pipeline — not the rasters — is what this repository provides: a reproducible, scripted workflow that any user can run on the CSIRO source data to obtain analysis-ready bioclimatic layers. The dataset is rigorously quality-controlled for spatial alignment and multicollinearity, providing a modeling-ready foundation for species distribution modeling (SDM), climate vulnerability assessments, and environmental applications across the Eastern Himalayas.

Comparative Feature Analysis

Preprocessing Time Comparison

Primary Beneficiaries

3.1 Source Data & Preprocessing

The primary inputs are monthly mean temperature, minimum temperature, maximum temperature, and precipitation grids obtained from Dorji et al. (2025). Source data were generated via downscaling to a 250 m resolution. Monthly files are organized by GCM, SSP, and time slice in the 01_raw_cmip6_data directory.

3.2 Bias Correction

3.3 Downscaling

The source data from Dorji et al. (2025) were downscaled to 250 m resolution using statistical downscaling methods. The source dataset provides monthly temperature (mean, minimum, maximum) and precipitation grids for Bhutan. Details of the downscaling technique are available in the source publication.

3.4 BIOCLIM Derivation

The 19 standard bioclimatic variables (BIO1–BIO19) were computed from monthly tasmin, tasmax, and pr rasters using the terra R package (Hijmans et al., 2023), following the standard definitions of O'Donnell & Ignizio (2012) and Nix (1986). These variables capture annual trends, seasonality, and extreme or limiting environmental factors.

Unit conversions are applied automatically by the pipeline: temperatures exceeding 150 (indicating Kelvin) are shifted by −273.15 to Celsius; precipitation values below 5 (indicating kg m⁻² s⁻¹) are multiplied by 2,628,000 to obtain mm month⁻¹. Quarterly variables (BIO8–BIO11, BIO16–BIO19) are computed using all 12 rolling 3-month windows via pixel-wise index operations — no approximation to fixed calendar quarters.

3.5 Ensemble Aggregation

Multi-model ensembles are calculated across the 10 GCMs. Aggregation uses equal (unweighted) weighting for all models, computing the arithmetic mean across all available GCMs for each time slice, scenario, and BIOCLIM variable.

3.6 Quality Control Acceptance Criteria

A file is only promoted to release status if it meets all criteria:

• Alignment: 100% pixel-perfect alignment with the Bhutan National Mask — verified by qc_validation.R
• CRS: Correct projection (EPSG:5266 — DRUKREF 03) embedded in metadata
• Extent: xmin: 125,694.5 · ymin: 2,954,438 · xmax: 460,694.5 · ymax: 3,125,938 (EPSG:5266)
• Completeness: No missing pixels within the Bhutan administrative boundary
• Range Check: Values within physically plausible bounds (e.g., BIO01 between −30°C and +40°C)
• Integrity: File integrity verified — see 06_quality_control/integrity_summary.json

The repository follows a structured 10-stage pipeline to ensure data provenance and reproducibility.

Pipeline Scripts

Running the Pipeline

The main entry point is bioclim_master.R. It auto-detects all GCMs, time periods, and SSP scenarios from the input directory structure — no configuration file needed.

# Minimal — auto-detect all models, periods, and scenarios
Rscript 10_scripts/r/bioclim_master.R \
  --input_root ./01_raw_cmip6_data/01_cmip6_monthly \
  --output_root ./03_bioclim_variables/01_bioclim_by_gcm

# Selective run — specify models and scenarios, tune memory
Rscript 10_scripts/r/bioclim_master.R \
  --input_root ./01_raw_cmip6_data/01_cmip6_monthly \
  --output_root ./03_bioclim_variables/01_bioclim_by_gcm \
  --models "ACCES_CM2,MIROC6,NorESM2-MM" \
  --scenarios "SSP245,SSP585" \
  --overwrite TRUE \
  --memfrac 0.5 \
  --tempdir "D:/terra_scratch"

# Run multicollinearity screening after BIO computation
Rscript 10_scripts/r/multicollinearity_screening.R \
  --project_root . \
  --time_slice "1986_2015" \
  --scenario "historical"

Output Filename Convention & GeoTIFF Specifications

All output rasters use a standardized naming pattern that encodes model, scenario, period, and variable:

Example: bhutan_cmip6_access_cm2_ssp245_2021_2050_bio01_v1_0.tif

6.1 Multicollinearity Screening

To avoid overfitting, screen predictors before training models:

• Consult 05_multicollinearity_analysis/vif_results.csv
• VIF < 5: Preferred for high-rigor models
• VIF < 10: Acceptable but requires justification
• Always prioritize ecological relevance over statistical correlation alone

6.2 Spatial Autocorrelation

High-resolution grids often exhibit strong spatial autocorrelation. We recommend:

• Spatial Thinning: Ensure occurrence points are not clustered within a single 250 m cell
• Block Cross-Validation: Use spatial blocks for model validation rather than random k-fold

6.3 Scenario Reporting Template

When reporting results, use a template such as:

6.4 Loading Pipeline Outputs in R

All pipeline outputs are standard compressed GeoTIFFs loadable with any R spatial library. Use terra (recommended, matches the pipeline) or sf/stars as needed.

library(terra)

# ── Single ensemble mean raster (historical BIO01) ───────────────────
# Filename pattern: bhutan_cmip6_ensemble_{stat}_{ssp}_{period}_{bio}_v1_0.tif
bio01_hist <- rast(
  "04_ensemble_products/ensemble_mean/1986_2015/historical/bhutan_cmip6_ensemble_mean_historical_1986_2015_bio01_v1_0.tif"
)

# ── Full 19-variable ensemble mean stack for a future scenario ────────
ssp245_2050 <- rast(
  list.files(
    "04_ensemble_products/ensemble_mean/2021_2050/ssp245/",
    pattern = "^bhutan_cmip6_ensemble_mean.*\\.tif$",
    full.names = TRUE
  )
)

# ── Ensemble standard deviation (uncertainty map) ─────────────────────
bio01_sd <- rast(
  "04_ensemble_products/ensemble_standard_deviation/2021_2050/ssp245/bhutan_cmip6_ensemble_standard_deviation_ssp245_2021_2050_bio01_v1_0.tif"
)

# ── VIF-screened predictor set for SDM ───────────────────────────────
# Consult 05_multicollinearity_analysis/selected_predictors.txt for retained variables.
# Load only those from the ensemble mean directory:
selected_vars <- readLines("05_multicollinearity_analysis/selected_predictors.txt")
mean_files <- list.files(
  "04_ensemble_products/ensemble_mean/1986_2015/historical/",
  pattern = "\\.tif$", full.names = TRUE
)
sdm_predictors <- rast(
  mean_files[sapply(mean_files, function(f)
    any(sapply(selected_vars, function(v) grepl(v, basename(f), fixed = TRUE)))
  )]
)

# ── Stack individual GCM outputs to quantify inter-model spread ───────
# Per-GCM file pattern: bhutan_cmip6_{model}_{scenario}_{period}_{bio}_v1_0.tif
gcm_bio01 <- rast(
  list.files(
    "03_bioclim_variables/01_bioclim_by_gcm/",
    pattern = "bio01_v1_0\\.tif$",
    recursive = TRUE, full.names = TRUE
  )
)
bio01_range <- max(gcm_bio01) - min(gcm_bio01)  # pixel-wise inter-model range

# ── Extract values at species occurrence points ───────────────────────
# occ_pts <- vect("occurrences.gpkg")
# climate_vals <- extract(sdm_predictors, occ_pts)

• File integrity verified — see 06_quality_control/integrity_summary.json
• Quality Control (QC) alignment report passed — see 06_quality_control/raster_alignment_report.txt
• GCM inventory (n=10) confirmed
• EPSG code confirmed (EPSG:5266 — DRUKREF 03 / Bhutan National Grid)
• Bias correction documented (Applied by Dorji et al. 2025 — see source publication)
• Nodata value confirmed (−9999)
• Spatial extent confirmed (xmin: 125,694.5 · ymin: 2,954,438 · xmax: 460,694.5 · ymax: 3,125,938)
• Ensemble weighting confirmed (Equal/unweighted mean across 10 GCMs)
• Version tag (v1.0.0) applied to all files
• Internal navigation links verified