H-HSM: Hierarchical Habitat Suitability Modeling Workflow

Source Code Supplement to:

Broad, E., C. J. Brown, N. Shackell, and K. Robert. 2026. “Hierarchical Habitat-Suitability Models Highlight the Influence of Benthic Seascapes on Snow Crab Distribution.” Landscape Ecology.

OVERVIEW

This code implements a two-tier Hierarchical Habitat Suitability Model
(H-HSM) using a covariate strategy to evaluate species distribution:

Tier 1 (Regional Level): Models presence/absence at ~11 km
scale using BNAM (Wang et al., 2018) model derrived ocean climate predictors.

Downscaling Step: Interpolates the Tier 1 prediction raster to a
0.2 km grid using bilinear interpolation via the 'terra' package.

Tier 2 (Landscape Level): Models distribution at a 0.2 km scale
using localized seabed terrain derivatives and the downscaled
climate suitability raster as an added covariate.

PREREQUISITES & R PACKAGES

Run the following command in R to install all required dependencies:

install.packages(c("terra", "dplyr", "usdm", "randomForest",
"caret", "sf", "blockCV", "pROC", "RColorBrewer"))

Data Coordination:

The test data coordinates (small male snow crab) are projected using the Albers
Equal Area coordinate system reference (AEAcrs).

DATA DIRECTORY SETUP

The raw data files are too large for GitHub and are hosted on the 
Federated Research Data Repository (FRDR).

1. Download the dataset from FRDR: [https://doi.org/10.20383/103.01702]
2. In your local repository directory, create a folder named "raw".
3. Extract/place the downloaded files into that folder:

Your_Project_Folder/
│
├── raw/
│   ├── snowcrab_regional.csv    
│   ├── snowcrab_landscape.csv  
│   ├── modelsm.tif         
│   └── 200m_AE_Depth.tif 

STEP-BY-STEP CODE WORKFLOW

Step 1: Custom RF Optimization
The script defines a custom 'caret' learner ("customRF") to
simultaneously tune both 'mtry' and 'ntree' variables.

Step 2: Tier 1 Regional Modeling

Filters occurrence records for the target demographic class.

Uses 'blockCV' to partition data into 10 spatially independent
folds using 11 km hexagonal blocks to mitigate autocorrelation.

Trains a Random Forest with downsampled, balanced presence-absence
classes across an iterative hyperparameter tuning grid (150-600 trees).

Step 3: Bilinear Resampling

Loads the generated 11.6 km Tier 1 regional prediction raster.

Uses terra::disagg(..., fact=58, method='bilinear') followed by
terra::resample() to scale it precisely to the 0.2 km landscape grid.

Step 4: Tier 2 Landscape Modeling

Extracts the downscaled climate suitability values at your
high-resolution sampling positions.

Appends this raster layer to the terrain data matrix as a covariate.

Splits data (75% train / 25% test), runs spatial block CV, and
trains the final micro-scale Random Forest model.

*Note* to create habitat suitability maps you will need to download and combine
rasters as single terra raster objects for each level of the H-HSM,
projected in the AEAcrs. Please see data sources linked in the manuscript
for BNAM and other pedictors. The 0.2 km bathyemtric surface is publically avalible 
from FRDR repositories, as used in the paper https://doi.org/10.20383/103.01702, and in its original format in https://doi.org/10.20383/103.01557.

CITATION for H-HSM code: 

Broad, E., C. J. Brown, N. Shackell, and K. Robert. 2026. “Hierarchical Habitat-Suitability Models Highlight the Influence of Benthic Seascapes on Snow Crab Distribution.” Landscape Ecology.