Initialize repo with Belize HRA/CV workflows in R by jagoldstein · Pull Request #1 · natcap/HRA_CV_plugin

jagoldstein · 2026-07-15T23:25:30Z

Summary

This PR initializes the repository by uploading the preliminary R workflows developed for the Belize context. These scripts establish the baseline data handling and logical workflows for both the Habitat Risk Assessment (HRA) and Coastal Vulnerability (CV) InVEST® models.

Key Changes

Uploaded foundational Belize-specific R scripts for HRA and CV postprocessing, summarizing, and integration.
Established the core baseline logic and spatial data structures for the project.

Next Steps

These workflows serve as the logical blueprint. Once this initial code review is approved and merged, the next phases will involve:

Translating the R script logic into structured Python modules.
Generalizing the codebase to fit a reusable InVEST plugin framework.

…ures and update initial HRA script

… and postprocess CV results

jagoldstein · 2026-07-16T20:30:54Z

@dcdenu4 from you here, I'm hoping for feedback, flags, or points I should be thinking about or looking out for as I prepare to translate this R code into Python with the eventual goal of building an InVEST plugin (first a Belize-specific one, but ultimately a generalized version). Are there best practices or recommended approaches I should be using? How is the depth and breadth of my annotated comments? What am I missing or should be thinking about? e.g.: the workflow is currently spread across 4 sequential scripts, should it all be compiled into 1? I would also welcome general review for code that could be streamlined or simplified.

@jade-md from you I'm looking more for a review of this revised R code and sanity-check (second set of eyes) to make sure it's faithfully executing the steps from your 4 original scripts and that I didn't break or miss something.

Thank you both!

dcdenu4

@jagoldstein, ultimately I think it would be helpful for me to understand how the HRA and CV models are run across scenarios before these scripts. When thinking of this as a plugin, I feel like there are a few options. The first is to turn these scripts into a post processing plugin as they are. The second is to create a workflow plugin that actually runs both the HRA and CV model and does this processing as part of a single "run".

My suggestion would be to translate these four scripts to their Python version and check results are comparable. For Python, this will likely entail moving the raster processing into Pygeoprocessing calls. For the vector manipulation, there are a lot of examples throughout InVEST of these types of operations.

In practice, these scripts probably don't need to all be separate. They seem to share some common information and variables that could nice to centralize and reuse. Though they could be broken into an HRA and CV module to help organize.

A plugin will have one entry point called "execute" that is run when someone hits the "run" button in the Workbench. You can do whatever you want within that "execute" function. So initially, you could imagine each of these scripts being run in sequence within the "execute" function.

I imagine that AI would be very good at helping translate from R to Python. If you do use AI in some way as a helper, I might suggest asking it to use the InVEST API as a style reference and to prefer to use Pygeoprocessing for handling raster and vector operations.

dcdenu4 · 2026-07-17T12:25:31Z

+library(reshape)
+```
+
+Note: The working directory needs to reflect your local machine. Set the single root path below to point at your local copy of the "Science" folder; every input and output path in this script is built from it, so there is only one line to change.


I suspect that we'll want to be more flexible in the plugin version. InVEST models now output an output registry, which is a JSON file mapping output paths. That could be an input into another part of a model process.

dcdenu4 · 2026-07-17T12:30:18Z

+    risk_layer <- reclass_risk == risk_value
+
+    # Visual check for this habitat/risk-class combination
+    plot(risk_layer, main = paste(habitat, risk_label))


Will this kind of run time plotting still be needed? Or could we include plots as part of a Report for users to check?

good question!... for @jade-md

given this script is to inform the MRV system for Belize, the vision is to produce graphics that will live an an onlione dashboard. IRD and the government of BZ is on the hook to produce the dashboard. In the interim, i think it would be good for us to output some graphics that decision makers can copy and paste in deck or other outlets to communicate the status of the ecosystems and people-related KPIs. That being said if you have suggestions on how to improve, i am open

That makes a lot of sense. So, not necessarily needing to display the plots during the model runs, but saving plots as outputs to be used later. That seems totally reasonable.

dcdenu4 · 2026-07-17T12:48:56Z

+
+# C/ Recompute the habitat rank for each CV point
+
+Based on equation 17 in the InVEST UG, recompute the habitat rank (`R_hab`) that incorporates the HRA-derived protective values, labeled here as `r_hra`. [ R_hab = 4.8 − 0.5·√( (1.5·(5−R_min))² + Σ(5−R_k)² − (5−R_min)² ) ]. Import the protective-value shapefile from section B, compute the equation 17 terms, and derive `r_hra`. Also add a no-habitat rank (`r_noh`, set to 5) to assess the effect of habitat.


I wonder if we could just import the HRA function that does this in the HRA model and use that instead of having to maintain a separate version of the implemented equation that could become desynced.

i love that idea Doug!

dcdenu4 · 2026-07-17T12:50:39Z

+         driver = "ESRI Shapefile", delete_layer = TRUE)
+```
+
+# D/ Recompute the exposure index for each CV point with the corrected Rhab


A lot of these sections and maybe even chunks of code within will likely be their own functions in Python or Pygeoprocessing calls.

jade-md · 2026-07-21T22:20:23Z

+
+# B/ Habitat 2016 - Stressors 2016: Calculate the habitat area under each risk class per coastal planning unit (CPU)
+
+Note that to calculate the habitat area we need to multiply the raster values by 0.25 because: 1) we sum the grid cells from HRA at 500m x 500m = 250,000 m2 and 2) we want to convert those values to km2. This factor is stored above as `cell_area_km2`.


This is where we may want the user to be able to define it if the resolution of HRA changes

jade-md · 2026-07-21T22:24:28Z

+habitats_dir <- file.path(wd_root, "01_hra_model/5_hra_outputs/2016_habitats")
+
+habitats <- c("corals", "seagrass", "mangroves")
+years    <- c("2016", "2020", "2025")


For the NP3 project we will only have 2016 and 2026 as the years of analysis... and the idea is to set up 2026 as the new baseline for the MRV system

jade-md · 2026-07-21T22:27:06Z

+wd_root <- "G:/Shared drives/GreenFin Belize/Science"
+habitats_dir <- file.path(wd_root, "01_hra_model/5_hra_outputs/2016_habitats")
+
+func_weights <- c("No Risk" = 1, "Low" = 1, "Medium" = 0.5, "High" = 0)


These assumptions are the ones we have been using for Belize over the years... but these could change for other places or over time as more research progresses... maybe something to consider to be users defined in the plugin

jade-md

Thank you Jesse for doing this and improving the flow of the codes already!! it is so cool to see this evolving!!
i made some comments flaggin places that we want to allow user to change parameters and not code them in hard.
I also made comments on places that are very belize specific
and places where we will only have 2 scenarios (2016 + 2026). For NP3 the goal is to set up 2026 as a baseline only. But i think by coding the comparison to 2016 we can provide the country w the tool to track change over time... although i would like to check in w Jess and Katie to see if we want to track change over time knowing the habitat maps change or if we want to compare to with and no habitat and see if the metric itself changes to avoid the habitat mapping challenge... to be continued

jade-md · 2026-07-22T20:05:29Z

+reclass_df <- c(0, 1,
+                1, 1,
+                2, 2,
+                3, 4)


This is where we could also have user defined option to choose the habitat rpotective value based on the risk

jade-md · 2026-07-22T20:06:45Z

+
+```{r Reclassify risk levels into habitat protective values and export as rasters}
+hab_risk <-raster("G:/Shared drives/GreenFin Belize/Science/01_hra_model/5_hra_outputs/2016_habitats/2016/outputs/RECLASS_RISK_seagrass.tif")
+reclass_df <- c(0, 4,


another user defined section for converting seagrass risk to seagrass protective role

jade-md · 2026-07-22T20:07:11Z

+
+```{r Reclassify risk levels into habitat protective values and export as rasters}
+hab_risk <-raster("G:/Shared drives/GreenFin Belize/Science/01_hra_model/5_hra_outputs/2016_habitats/2016/outputs/RECLASS_RISK_mangroves.tif")
+reclass_df <- c(0, 1,


Same here for mangroves risk to habitat protective role

jade-md · 2026-07-22T20:07:47Z

+linkages_dir   <- file.path(wd_root, "02_cv_model/hra_cv_linkages")
+out_dir        <- file.path(wd_root, "02_cv_model/cv_outputs/2016_habitats/2016")
+
+# Each habitat's protective distance (m), raster file, and output column name


The buffer size as well are distances we use for belize but that could be redefined by users

jade-md · 2026-07-22T20:09:06Z

+
+For coral, divide the coral habitat by the fringing and barrier reefs for CV. Import a coral map which distinguishs between fringing and barrier reefs, where nearshore reef = 1 and barrier reef = 2. In one map, we set the coral habitat that is barrier reef as NA and in the other we set the coral habitat that is nearshore as NA. Then we print the maps.
+
+```{r Export rasters with the habitat protective values & correct for the barrier reef}


this step splits thew reef into nearshore and barrier for belize - super belize specific

jade-md · 2026-07-22T20:14:20Z

+
+# Each year's population raster (2025 uses 2020)
+ppp_files <- c("2016" = "bz_ppp_2016_2km_mean.tif",
+               "2020" = "bz_ppp_2020_2km_mean.tif",


we will drop this scenario later

jade-md · 2026-07-22T20:14:45Z

+delta_by_cpu <- list()
+
+for (y in c("2016", "2020", "2025")) {
+


remove 2020 and label 2016 and 2026 for NP3

jade-md · 2026-07-22T20:15:15Z

+```{r Compute km and people no longer at high risk due to habitat, per CPU}
+high_delta_by_cpu <- list()
+
+for (y in c("2016", "2020", "2025")) {


drop 2020 here too

jade-md · 2026-07-22T20:15:41Z

+  # Keep points that were high-risk (eic_hig == 1) AND improved due to habitat
+  pts <- pts[pts$eic_hig == 1 & pts$risk_dlt %in% improved, ]
+
+  # Extract people per point; scale by 25 (0.25 km point, 2 km mean raster)


this could be another user defined factor as it will vary based on the CV scale and population layer

jade-md · 2026-07-22T20:16:38Z

+  pts <- pts[pts$eic_hig == 1 & pts$risk_dlt %in% improved, ]
+
+  # Extract people per point; scale by 25 (0.25 km point, 2 km mean raster)
+  ppl_raster <- raster(file.path(ppp_dir, ppp_files[[y]]))


for this layer Jess asked me to convert the pop count to density by using focal stats over 2km cicle radius. Maybe this is a step we want to include in here to be able to input the worlpop layer directly?

jagoldstein added 3 commits July 15, 2026 10:23

add script to postprocess and summarize HRA results

25f4663

add script to postprocess risk summaries into combined tables and fig…

6eb6ad7

…ures and update initial HRA script

Add 2 remaining scripts from Belize: convert HRA outputs to CV inputs…

4d44825

… and postprocess CV results

jagoldstein requested review from dcdenu4 and jade-md July 15, 2026 23:25

dcdenu4 reviewed Jul 17, 2026

View reviewed changes

dcdenu4 approved these changes Jul 17, 2026

View reviewed changes

jade-md reviewed Jul 21, 2026

View reviewed changes

jade-md reviewed Jul 22, 2026

View reviewed changes


		# C/ Recompute the habitat rank for each CV point

		Based on equation 17 in the InVEST UG, recompute the habitat rank (`R_hab`) that incorporates the HRA-derived protective values, labeled here as `r_hra`. [ R_hab = 4.8 − 0.5·√( (1.5·(5−R_min))² + Σ(5−R_k)² − (5−R_min)² ) ]. Import the protective-value shapefile from section B, compute the equation 17 terms, and derive `r_hra`. Also add a no-habitat rank (`r_noh`, set to 5) to assess the effect of habitat.


		# B/ Habitat 2016 - Stressors 2016: Calculate the habitat area under each risk class per coastal planning unit (CPU)

		Note that to calculate the habitat area we need to multiply the raster values by 0.25 because: 1) we sum the grid cells from HRA at 500m x 500m = 250,000 m2 and 2) we want to convert those values to km2. This factor is stored above as `cell_area_km2`.


		For coral, divide the coral habitat by the fringing and barrier reefs for CV. Import a coral map which distinguishs between fringing and barrier reefs, where nearshore reef = 1 and barrier reef = 2. In one map, we set the coral habitat that is barrier reef as NA and in the other we set the coral habitat that is nearshore as NA. Then we print the maps.

		```{r Export rasters with the habitat protective values & correct for the barrier reef}

		delta_by_cpu <- list()

		for (y in c("2016", "2020", "2025")) {

Uh oh!

Conversation

jagoldstein commented Jul 15, 2026

Summary

Key Changes

Next Steps

Uh oh!

jagoldstein commented Jul 16, 2026

Uh oh!

dcdenu4 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jade-md left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants