Coastal vulnerability: Calculating habitat rank

Software Support
1

I’m trying to modify the habitat rank results from the coastal vulnerability model (my goal is to recalculate the habitat rank for each shoreline segment without each of the individual habitat types, to see where each habitat type has the strongest effect on the coastal exposure index), but am having trouble recreating the current habitat rank from the model outputs.

I created a csv file with one record for each shoreline segment, containing the final habitat rank (from 1_d_natural_habitats.csv) and the presence/absence of each habitat type (using the extract values to points tool in ArcGIS to extract the values from the 1_d_HAB_XX_[habitat]_influence_on_shore.tif files).

Next, I wrote an R script that calculates the habitat rank for each shoreline segment from the individual habitat ranks, using the equation from the model user’s guide. When I calculate the habitat rank in the same way that (I think) the model does, I don’t get the same results. For example, one shoreline segment is ‘within range’ of three habitat types, with ranks of 1, 4, and 2. I calculate that the final habitat rank should be 2.226092, but the habitat rank output from the model is 1.40884.

I’ve double-checked the individual habitat ranks that I’m using (they match what was used in the model run). I also tried calculating the habitat rank with ‘missing’ habitats (those out of range of the shoreline segment) set to 5 rather than just excluding them from the calculation, but still didn’t get the habitat ranks calculated within the model. In addition, neither my calculated habitat rank nor the model output rank matches the ranks in Appendix B of the User’s Guide (my calculated rank is close, but not exact), even though my protection distance and rank for habitats matches those suggested in the user’s guide.

I’m using version 3.6 of the coastal vulnerability model. Any help understanding why my calculated habitat ranks differ from the model outputs would be great! I can share relevant outputs and the R script if that’s helpful. Thanks!

2

This sounds like a good plan. The habitat rank equation in the User’s Guide is admittedly hard to parse. Is it possible the discrepancy you see is an error in translating the equation into R code? Here’s the python implementation of the equation taken directly from the current version of the model, where row is the series of individual habitat ranks at a shore point.

def _calc_Rhab(row):
    """Equation 4 from User's Guide."""
    sum_sq_rank = 0.0
    min_rank = 5  # 5 is least protection
    for r in row:
        if r < min_rank:
            min_rank = r
        sum_sq_rank += (5 - r)**2

    if sum_sq_rank > 0:
        r_hab_val = max(
            1.0, 4.8 - 0.5 * (
                (1.5 * (5 - min_rank))**2 + sum_sq_rank -
                (5 - min_rank)**2)**0.5)
    else:
        r_hab_val = 5.0
    return r_hab_val

print(_calc_Rhab((1, 4, 2)))

This example yields the 1.40884 value you mentioned.

You probably have your reasons for using version 3.6, but the format of the outputs of version 3.7 are likely easier to work with for this specific task. See intermediate/habitats/habitat_protection.csv which sounds like it is already the table you are creating yourself in the “extract values from raster” step you describe.

If you still haven’t figured out the source of the discrepancy and you would like to post your R script, I’ll take a look!

3

Thanks so much, it’s incredibly helpful to see the original code. I found two discrepancies in my translation of equation 4 into R:

  • I interpreted the maximum rank (to be weighted 1.5 times higher than the other habitat ranks) as the rank with the highest value (closest to 5) instead of the rank with the maximum protective value (lowest numerical value)
  • I had an order of operations issue in the calculation applying the 1.5x weight to that rank - I was applying the exponent before multiplying by 1.5 instead of after.

I’m now able to duplicate the original habitat rank results. As for the model versions, you’re right, the 3.7 output formats would simplify this process. I’ve spent a lot of time tweaking the 3.6 parameters to get the sheltered/exposed shoreline segments right, and when I did a test run in 3.7 it looked like that part of the model has been changed (no depth filter?). Since I don’t have time to figure that out at the moment, I’m sticking with 3.6 for now.

1 Like