7 Habitat Suitability Index models with ecorest
This module will teach you how to use the ecorest package in R to efficiently and reproducibly carry out HSI modeling and generate useful output like this:
Authors: Darixa Hernandez-Abrams (writing, code), Kyle McKay (writing, code), Ed Stowe (writing, code, editing), ; Background from the draft ecorest tech report, authors: Kyle McKay, Darixa Hernández-Abrams, Rachel Nifong, Todd Swannack
Last update: 2025-10-30
Acknowledgements:
7.1 Learning objectives
- Understand the benefits of conducting HSI modeling programmatically with
ecorest - Learn how to use
ecorestto implementblue bookHSI models - Build a custom HSI model using
ecorest
7.2 Background on HSI modeling with ecorest
Before ecorest, there was no standard platform for computing outcomes from habitat or other index models. Then and now, users often develop ad hoc spreadsheet models, but these can be highly prone to numerical errors (McKay 2009). The common quantity-quality structure of index models, however, provides an opportunity to develop a consistent, error-checked index modeling calculator adaptable to a variety of applications across USACE. That tool is ecorest, which can greatly increase the efficiency of performing habitat modeling as well as decrease the likelihood of computational errors.
The ecorest package can be used to quickly apply approximately 350 of the HSI models developed by the U.S. Fish and Wildlife Service in the 1980s, often referred to as the “blue book” models, but is flexible enough to carry out any user-defined HSI model, including models that appear in the literature or custom models developed for specific projects.
In this tutorial, we will demonstrate how the ecorest package can be used for two different HSI modeling scenarios, one that relies on a blue book model, and another that uses a custom HSI model.
7.3 Run an HSI analysis with ecorest using a blue book model
For this scenario, we will use ecorest to carry out a modeling analysis that was originally implemented by the Sacramento District as part of the Delta Islands and Levees Feasibility Study. In this study, the Marsh Wren HSI model was used to examine the ecosystem effects of intertidal marsh restoration. The overall Feasibility Study can be found here, and the appendix with details on the habitat evaluation procedure can be found here.
To conduct the HSI modeling in this scenario, we first need to add the necessary packages to our R session.
The ecorest package contains a list object called HSImodels. Each of the 349 elements of the list is a dataframe representing one of the U.S. Fish and Wildlife Service’s habitat suitability index (HSI) models. Each HSI model consists of multiple independent variables (e.g., percent canopy cover) and their associated habitat suitability values ranging from 0 to 1. Categorical input variables are coded as letters. All models can be called in R with the corresponding model name using HSImodels$modelname (e.g., HSImodels$barredowl).
To look at the HSI models that are available, run the following code:
names(HSImodels)## [1] "alewifeJuv" "alewifeJuvAndSAEL"
## [3] "alewifeSAEL" "americanalligatorNontidal"
## [5] "americanalligatorTidal" "americanblackduckWinteringEVegWetland"
## [7] "americanblackduckWinteringNCapeCod" "americanblackduckWinteringSCapeCod"
## [9] "americancoot" "americaneiderBreeding"
## [11] "americanoysterGulfofMexModifier" "americanoysterGulfofMexTypical"
## [13] "americanshadEstu" "americanshadRiv"
## [15] "americanwoodcockWinteringForestedDry" "americanwoodcockWinteringForestedMoist"
## [17] "americanwoodcockWinteringForestedWet" "americanwoodcockWinteringShrubDry"
## [19] "americanwoodcockWinteringShrubMoist" "americanwoodcockWinteringShrubWet"
## [21] "arcticgrayling" "arcticgraylingAdultJuv"
## [23] "arcticgraylingSEF" "atlanticcroakerLATideLt0.5m"
## [25] "atlanticcroakerLATideMt0.5m" "atlanticcroakerOtherTideLt0.5m"
## [27] "atlanticcroakerOtherTideMt0.5m" "atlanticcroakerWetlandLATideLt0.5m"
## [29] "atlanticcroakerWetlandLATideMt0.5m" "atlanticcroakerWetlandOtherTideLt0.5m"
## [31] "atlanticcroakerWetlandOtherTideMt0.5m" "bairdssparrow"
## [33] "baldeagleBreeding" "barredowl"
## [35] "beaverLacAreaLt8ha" "beaverLacAreaMtoe8ha"
## [37] "beaverPalu" "beaverRiv"
## [39] "beltedkingfishLenticConstWave" "beltedkingfishLenticNoConstWave"
## [41] "beltedkingfishLotic" "bigmouthbuffaloLacNoSal"
## [43] "bigmouthbuffaloLacSal" "bigmouthbuffaloRivNoSal"
## [45] "bigmouthbuffaloRivSal" "blackbear"
## [47] "blackbelliedwhistlingduck" "blackbrant"
## [49] "blackbullheadLac" "blackbullheadRiv"
## [51] "blackcappedchickadeeFoodCanH" "blackcappedchickadeeFoodCanVol"
## [53] "blackcrappieLacNoSal" "blackcrappieLacSal"
## [55] "blackcrappieRivNoSal" "blackcrappieRivSal"
## [57] "blacknosedaceLac" "blacknosedaceRiv"
## [59] "blackshoulderedkite" "blacktailedprairiedog"
## [61] "bluegillLac" "bluegillRiv"
## [63] "bluegrouse" "blueherringJuv"
## [65] "blueherringJuvAndSAEL" "blueherringSAEL"
## [67] "bluewingedtealBreeding" "bobcatLt4ha"
## [69] "bobcatMtoe4ha" "brewerssparrow"
## [71] "brooktroutLacAllLtoe15C" "brooktroutLacAllMt15C"
## [73] "brooktroutRivAllLtoe15CLtoe5mEC" "brooktroutRivAllLtoe15CMt5mEC"
## [75] "brooktroutRivAllMt15CLtoe5mEC" "brooktroutRivAllMt15CMt5mEC"
## [77] "brownshrimpNGulfofMex" "brownthrasher"
## [79] "browntroutCompLtoe10C" "browntroutCompMt10C"
## [81] "browntroutLimitLtoe10C" "browntroutLimitMt10C"
## [83] "bullfrog" "cactuswren"
## [85] "canvasbackBreeding" "channelcatfishLac"
## [87] "channelcatfishRiv" "chinooksalmonComp5to10CSand"
## [89] "chinooksalmonComp5to10CSilt" "chinooksalmonCompLtoe5CSand"
## [91] "chinooksalmonCompLtoe5CSilt" "chinooksalmonCompMt10CSand"
## [93] "chinooksalmonCompMt10CSilt" "chinooksalmonLimit5to10CSand"
## [95] "chinooksalmonLimit5to10CSilt" "chinooksalmonLimitLtoe5CSand"
## [97] "chinooksalmonLimitLtoe5CSilt" "chinooksalmonLimitMt10CSand"
## [99] "chinooksalmonLimitMt10CSilt" "chumsalmonAlevin"
## [ reached getOption("max.print") -- omitted 251 entries ]Another important object that is contained within the package is the HSImetadata dataframe. This contains information on each HSI model, including how the different suitability variables within a model should be aggregated to calculate an overall HSI value for a project alternative. The following code generates a small subset of the variables and rows of this dataframe.
HSImetadata %>%
dplyr::select(model, submodel, Eqtn) %>% # Select which columns to display
slice(1,4:6) # Select which rows to display## model submodel Eqtn
## 1 alewifeJuv Juvenile life stage model min(CF, CWQ)
## 2 americanalligatorNontidal Nontidal wetland (CCB*CCN)^(1/2)
## 3 americanalligatorTidal Tidal wetlands (CCB*CCN)^(1/2)
## 4 americanblackduckWinteringEVegWetland Estuarine vegetated wetlands CF7.3.1 Explore the Marsh Wren HSI model
To use the Marsh Wren HSI model, we extract it from the HSImodels object.
wren_hsi <- HSImodels$marshwrenPrint the Marsh Wren model to the console to look at its structure:
print(wren_hsi)## emerg.hydrophytes.class emerg.hydrophytes.SIV emerg.herb.can.cov.pct emerg.herb.can.cov.SIV
## 1 a 1.0 0 0.0
## 2 b 0.5 50 0.1
## 3 c 0.1 80 1.0
## 4 d 0.0 100 1.0
## avg.wtr.d.cm avg.wtr.d.cm.SIV woody.can.cov.pct woody.can.cov.SIV
## 1 0 0 0 1
## 2 15 1 100 0
## 3 40 1 NA NA
## 4 NA NA NA NAAs you can see, HSImodels provides a series of suitability curves ordered as parameter breakpoints and associated suitability indices for each parameter. Each column containing information about an environmental variable (e.g., emerg.hydrophytes.class) is followed by a column (with a name ending in .SIV) of suitability index values for the associated environmental predictor.
The marsh wren HSI model can also be viewed graphically. The HSIplotter function creates a JPEG image file in the images folder of your working directory. The user selects the name of the file, but should keep the “.jpg” file extension.
HSIplotter(wren_hsi, "images/marshwren_hsi.jpg")If you open the plot file in your working directory, it should look like this:
7.3.2 Calculate suitability values from environmental variables
To calculate the acres of habitat that a project would produce, we need estimates of the environmental variables that make up the chosen HSI model. For the Marsh Wren, these variables include:
- Growth form of emergent hydrophytes
- Canopy cover of emergent herbaceous vegetation
- Mean water depth
- Canopy cover of woody vegetation
Values for these variables can be imported from a .csv file, or can be entered manually as we do below. Here we take values from the row of the table that corresponds with target year 1 (TY1). Environmental variables must be in the same order as they appear in the HSI model. Categorical variables, such as emergent hydrophyte growth form, are coded as letters.
env_vars <- tibble(v1 = "a",
v2 = 4,
v3 = 99.1,
v4 = 0)In order to calculate habitat units, we also need to provide the project area in acres. Here we provide the area listed for TY1 above.
area = 34Next, we use the SIcalc function to calculate the suitability value for each variable in the model. The SIcalc function computes suitability indices using two inputs: the suitability model (wren_hsi) and the project-specific environmental variables that we’ve provided (env_vars).
## [1] 1.000 0.008 1.000 1.0007.3.3 Aggregate HSI value and habitat units
A crucial part of HSI modeling is determining how the SI values of our different variables will be combined to create the overarching HSI score. SI values can be combined in numerous ways (e.g., geometric mean, minimum value, etc.). Most blue book models, however, identify a species/submodel-specific way to combine SI values, according to the specific attributes of the species. Therefore, when a blue book model is being used, the HSIeqtn function should be used to calculate the cumulative HSI score using the appropriate equation. This function requires the name of the HSI model, the variables for each of the individual metrics that we’ve calculated, and the HSImetadata file, where the equation is stored.
## [1] 0.2Once we have calculated our cumulative HSI score from our separate SI values, we can lastly calculate the habitat units associated with the given project.
Total_hab <- total_HSI*area
print(Total_hab)## [1] 6.87.3.4 Habitat units for multiple years or alternative projects
If multiple project alternatives exist, or there are multiple years of environmental data, we can perform the same operations as above using the below code. Here, a for-loop allows us to iterate through each of the scenarios/years to calculate SI values and habitat units
# Create dataframe of environmental variables for 4 scenarios/years
new_env_vars <- tibble(V1 = rep("a", 4),
V2 = c(5, 8.5, 12, 16),
V3 = c(94.5, 89.9, 85.3, 80.8),
V4 = c(0, 0, 0, 0))
nyears <- nrow(new_env_vars) # Number of years/alternatives
nvars <- ncol(new_env_vars) # Number of environmental variables
new_area = c(68, 102, 136, 170)
# Create an empty matrix and vector to hold the results of the SI and HSI calculations, respectively
new_si_vars <- matrix(0, nrow=nyears, ncol=nvars)
new_total_HSI <- rep(NA, 4)
# Use a for-loop to calculate SI values and HSI scores for each year/scenario
for (i in 1:nyears) {
new_si_vars[i, ] <- SIcalc(wren_hsi, unlist(new_env_vars[i, ]))
new_total_HSI[i] <- HSIeqtn("marshwren", new_si_vars[i, ], HSImetadata)
}
# Create a dataframe of the new SI values
new_si_df <- as_tibble(new_si_vars) %>%
set_names(c("SI_V1","SI_V2","SI_V3","SI_V4"))
# Create an output table with all of the input and output data
output_tab_wren <- new_env_vars %>%
bind_cols(new_si_df) %>%
mutate(output_HSI = round(new_total_HSI, 2),
total_acres = new_area,
HUs = output_HSI*total_acres)
# Print the table using the knitr package
knitr::kable(output_tab_wren,
digits = 2,
padding = 2,
caption = "Table 2. Habitat Suitability Index (HSI) Variables and
Resulting Outputs for the Marsh Wren Model With-Project.")| V1 | V2 | V3 | V4 | SI_V1 | SI_V2 | SI_V3 | SI_V4 | output_HSI | total_acres | HUs |
|---|---|---|---|---|---|---|---|---|---|---|
| a | 5.0 | 94.5 | 0 | 1 | 0.01 | 1 | 1 | 0.22 | 68 | 14.96 |
| a | 8.5 | 89.9 | 0 | 1 | 0.02 | 1 | 1 | 0.26 | 102 | 26.52 |
| a | 12.0 | 85.3 | 0 | 1 | 0.02 | 1 | 1 | 0.29 | 136 | 39.44 |
| a | 16.0 | 80.8 | 0 | 1 | 0.03 | 1 | 1 | 0.32 | 170 | 54.40 |
As you can see, the results we’ve generated using the ecorest package are identical to the ones in the report:
The beauty of carrying this out in R using a script is that we can change an input value and reproduce all the results in a matter of seconds by re-running the code, instead of having to go through all the steps manually. Plus, the ecorest package provides instant access to hundreds of quality-controlled blue book models.
7.4 HSI modeling with custom model
Sometimes, a project will use a custom index model instead of a blue book model to calculate habitat units. For example, for its South San Francisco Bay Shoreline Study, USACE-SPN collaborated with ERDC’s Coastal Ecology and Integrated Ecological Modeling teams to develop a quantitative model for tidal marsh restoration in subsided former baylands. The model is sensitive to multiple restoration measures, including the import of beneficially used material for bathymetry lifts, ecotones and other refuges.
They used the ecorest package for their habitat model calculations.
7.4.1 Explore model and data
Import a .csv of the marsh model, print the names of the variables, and plot the model. Note: for the plotting function to work, the model must be read in using the base R function read.csv instead of the tidyverse function read_csv
# Import marsh HSI model
marsh_model<-read.csv("data/SF-TM-HEI_InputData.csv")
print(names(marsh_model))## [1] "percent.tidal.range"
## [2] "percent.tidal.range.SIV"
## [3] "Average.Marsh.Age"
## [4] "Average.Marsh.Age.SIV"
## [5] "Percent.Marsh.within.500.or.200"
## [6] "Percent.Marsh.within.500.or.200.SIV"
## [7] "Percent.Marsh.Shoreward.Edge.with.Ecotone.Transition"
## [8] "Percent.Marsh.Shoreward.Edge.with.Ecotone.Transition.SIV"
HSIplotter(marsh_model,"images/marsh_model.jpg")
Next, read in the environmental data associated with the alternative project designs.
# Dataset of environmental conditions associated with project alternatives
marsh_alts <- read.csv("data/SF-TM-HEI_FieldData.csv")
print(marsh_alts)## alt V1 V2 V3 V4
## 1 FWOP 0 0 0 0
## 2 Alt1a 20 100 0 0
## 3 Alt1b 50 100 0 0
## 4 Alt1c 100 100 0 0
## 5 Alt2a 20 100 20 100
## 6 Alt2b 50 100 20 100
## 7 Alt2c 100 100 20 100
## 8 Alt3a 20 100 100 100
## 9 Alt3b 50 100 100 100
## 10 Alt3c 100 100 100 100
## 11 Alt3d 100 20 100 100
## 12 Alt3e 100 10 100 100
## 13 Alt3f 100 1 100 100
## 14 Alt3g 100 0 100 100
## 15 Alt3h 100 0 0 07.4.2 Calculate individual suitability indices
Again, we will create an empty matrix to store the results and calculate the SI values for each parameter using SIcalc.
# Create a dataframe of just environmental variables associated with each alternative
env_marsh <- dplyr::select(marsh_alts, -alt)
nalternatives<-length(marsh_alts$alt) # Number of alternatives
# Create empty matrix to store SI values for each alternative
si_marsh <- matrix(NA, nrow=nalternatives, ncol=4)
# Calculate the SI values for each alternative using a for-loop
for(i in 1:nalternatives) {
si_marsh[i, ] <- SIcalc(marsh_model, env_marsh[i, ])
}
# Convert matrix of SI values to dataframe and name the columns
si_marsh_df <- as_tibble(si_marsh) %>%
set_names(paste0("SI_", names(env_marsh)))7.4.3 Calculate the total HSI score
Here, the analysis team developed a custom formula to combine the different HSI scores:
Cube root((Tidal connectivity * Marsh Age) * (Ecotone + Refuges)/2)
In R syntax, that code looks like this
((((SIV2*SIV1)(SIV4+SIV3)/2))^(1/3))
To calculate the cumulative score, we’ll use the mutate function to create a new column in our dataframe of suitability index values, and calculate the total HSI values from the other columns.
7.4.4 Cumulative HSI and habitat units
Next, we’ll calculate the habitat units associated with the project, assuming a total project area of 300 acres.
# Calculate habitat units by multiplying total HSI by project area (300 acres)
area_marsh <- 300
si_total_hu<-si_total %>%
mutate(HUs = area_marsh*total_hsi)Finally, we can create an output table with all of the data.
output_tab <- bind_cols(marsh_alts, si_total_hu) %>%
set_names(c("Alternative", "PTR", "MA", "HTR", "E", "PTR-SI", "MA-SI",
"HTR-SI", "E-SI", "HSI-Total", "Marsh Units"))
knitr::kable(output_tab, caption="Table 3. SF Salt Marsh HSI Values and Marsh Units.",
align="c",
digits = 2)| Alternative | PTR | MA | HTR | E | PTR-SI | MA-SI | HTR-SI | E-SI | HSI-Total | Marsh Units |
|---|---|---|---|---|---|---|---|---|---|---|
| FWOP | 0 | 0 | 0 | 0 | 0.00 | 0.03 | 0.5 | 0.5 | 0.00 | 0.00 |
| Alt1a | 20 | 100 | 0 | 0 | 0.03 | 1.00 | 0.5 | 0.5 | 0.23 | 69.62 |
| Alt1b | 50 | 100 | 0 | 0 | 0.20 | 1.00 | 0.5 | 0.5 | 0.46 | 139.25 |
| Alt1c | 100 | 100 | 0 | 0 | 1.00 | 1.00 | 0.5 | 0.5 | 0.79 | 238.11 |
| Alt2a | 20 | 100 | 20 | 100 | 0.03 | 1.00 | 0.6 | 1.0 | 0.27 | 81.43 |
| Alt2b | 50 | 100 | 20 | 100 | 0.20 | 1.00 | 0.6 | 1.0 | 0.54 | 162.87 |
| Alt2c | 100 | 100 | 20 | 100 | 1.00 | 1.00 | 0.6 | 1.0 | 0.93 | 278.50 |
| Alt3a | 20 | 100 | 100 | 100 | 0.03 | 1.00 | 1.0 | 1.0 | 0.29 | 87.72 |
| Alt3b | 50 | 100 | 100 | 100 | 0.20 | 1.00 | 1.0 | 1.0 | 0.58 | 175.44 |
| Alt3c | 100 | 100 | 100 | 100 | 1.00 | 1.00 | 1.0 | 1.0 | 1.00 | 300.00 |
| Alt3d | 100 | 20 | 100 | 100 | 1.00 | 0.95 | 1.0 | 1.0 | 0.98 | 294.91 |
| Alt3e | 100 | 10 | 100 | 100 | 1.00 | 0.80 | 1.0 | 1.0 | 0.93 | 278.50 |
| Alt3f | 100 | 1 | 100 | 100 | 1.00 | 0.13 | 1.0 | 1.0 | 0.51 | 153.90 |
| Alt3g | 100 | 0 | 100 | 100 | 1.00 | 0.03 | 1.0 | 1.0 | 0.29 | 87.72 |
| Alt3h | 100 | 0 | 0 | 0 | 1.00 | 0.03 | 0.5 | 0.5 | 0.23 | 69.62 |
7.4.5 Linking steps using pipes
The last several operations that we conducted above can be done in a single “tidy” way by using the pipe operator %>%.
The final line of code checks demonstrates that the two output dataframes are identical.
area_marsh = 300
output_tab_2 <- as_tibble(si_marsh) %>%
set_names(paste0("SI_", names(env_marsh))) %>%
mutate(total_hsi = ((((SI_V4+SI_V3)/2)*(SI_V2*SI_V1))^(1/3)),
HUs = area_marsh*total_hsi) %>%
bind_cols(marsh_alts,.) %>% # This code is slightly different so that marsh_alts data will come first; the "." symbol indicates everything that has come before this function
set_names(c("Alternative", "PTR", "MA", "HTR", "E", "PTR-SI", "MA-SI",
"HTR-SI", "E-SI", "HSI-Total", "Marsh Units"))
all.equal(output_tab, output_tab_2)## [1] TRUE7.5 Summary
-
ecorestrepresents a flexible, reproducible tool for carrying out habitat modeling in the R environment, with countless advantages over ad-hoc spreadsheet approaches - The
ecorestpackage can instantly utilize ~ 350 USFWS blue book models -
ecorestalso enables users to quickly implement custom HSI models - Using
ecorestfor HSI modeling can greatly increase efficiency and decrease the likelihood of errors when implementing habitat models