R/FloodVulnerabilityStudy.R

In essatech/MNAI.CPBT: Municipal Natural Assets Initiative - Coastal Protection Benefit Toolbox

#' @title Flood Vulnerability Study
#'
#' @description Calculates the number of flooded structures, median flood
#' depth and total structural damage across a range of water level elevations.
#'
#' @param Bldgs Spatial polygon of building footprints of class sf and dataframe.
#' See data(Bldgs) for an example input format.
#' @param group_col (Optional) Column name of group column name in building
#' data frame. If used results will be summarized by region.
#' @param TopoBathy TopoBathy digital elevation model of class RasterLayer.
#' @param start_elevation Starting lower water level elevation for sequence.
#' @param end_elevation Ending upper level elevation for sequence.
#' @param step_size Water level elevation step size for sequence.
#' @param simulation_name Character. Simulation name for export.
#' @param simulation_name Character. Simulation name for export.
#' @param dir_output Local output directory for CPBT simulation results.
#'
#' @details The flood vulnerability is a computationally intensive yet simple
#' tool to allow communities to better understand the vulnerability of their
#' communities to coastal flooding. The function works by gradually raising
#' the water level by a fixed amount and then recalculating the structural
#' damage cost of a flood at that water level. Users define a range of
#' elevation e.g., 0 – 5 meters above chart datum and a step size e.g.,
#' 0.15 m and then the function calculates to total damage cost at each
#' 0.15 interval. A html report is exported with three plots showing the
#' depth-damage curves used in the assessment, the number of affected
#' structures at each water level and the total damage cost at each
#' water level. The intent of this tool is to look at overall damage cost
#' to different storms. If the number of affected structures and total
#' damage costs quickly climb upwards at low water levels, then users can
#' assume that the area is highly vulnerable to coastal flooding.
#' Conversely, if damage costs do not show any significantly values
#' until the water level reaches extreme levels, then users can
#' conclude that the community is naturally resilient to flooding.
#' Users can also evaluate tipping points at which the damage cost
#' climbs significantly for a given water level. For example, if a
#' community is expect to have no structural damage below flood water
#' levels of 2.1 m CD, but then experience millions of dollars in
#' damages between water levels of 2.1 – 2.5m then users can conclude
#' that any forces pushing the water levels above 2 m represent a
#' significant concern for their communities.
#'
#' @return Exports results to a specified folder file path.
#'
#' @export
FloodVulnerabilityStudy <- function(
  Bldgs = NA,
  group_col = NA,
  TopoBathy = NA,
  start_elevation = 2,
  end_elevation = 10,
  step_size = 0.2,
  simulation_name = "My Simulation",
  dir_output = getwd()
) {


  print(nrow(Bldgs))

  HAZUS <- HAZUS


  #----------------------------
  # Process by region
  #----------------------------

  if(is.na(group_col)){
    Bldgs$GROUP <- 1
  } else {
    gval <- Bldgs[,c(group_col)]
    sf::st_geometry(gval) <- NULL
    gval <- as.character(gval[,1])
    Bldgs$GROUP <- gval
  }



  # Extract mean elev for each bldg
  bdsp <- methods::as(Bldgs, "Spatial")

  print("Extract elev ...")
  val <- raster::extract(TopoBathy, bdsp, fun=mean)

  bdsp$elev <- val

  bdsp$Detail_ID <- 1:nrow(bdsp)






  #----------------------------
  # Loop through depth sequence
  #----------------------------

  ugroups <- unique(bdsp$GROUP)
  all_g <- list()

  for(g in 1:length(ugroups)) {

  this_group <- ugroups[g]

  bdsp2 <- bdsp[which(bdsp$GROUP == this_group), ]

  dseq <- seq(start_elevation, end_elevation, by=step_size)
  output <- list()

  for(i in 1:length(dseq)){

    this_depth <- dseq[i]

    depth_datum <- this_depth

    this_set <- bdsp2[which(bdsp2$elev < depth_datum), ]


    # nothing flooded
    if(nrow(this_set) == 0){

      add_row <- data.frame(
        depth = this_depth,
        n_blds = 0,
        total_cost = 0,
        median_depth = 0
      )

      output[[i]] <- add_row
      next

    }

    this_set$DamageCost <- NA
    this_set$DamagePct <- NA
    this_set$DamageDepth <- NA


    # Perform rolling join to match closest value from lookup table

    for(b in 1:nrow(this_set)) {

      this_b <- this_set[b, ]

      flood_depth <- depth_datum - this_b$elev

      # Curve for bld
      curv <- HAZUS[which(HAZUS$DDID == this_b$DDID),]

      # Reference DD curve
      dmg_est <- stats::approx(x = curv$depth_m,
                        y = curv$damage,
                        xout = flood_depth,
                        rule = 2)

      dmg_est <- round(dmg_est$y, 3)




      if(is.na(dmg_est)){

        this_set$DamageCost[which(this_set$Detail_ID == this_b$Detail_ID)] <- NA
        this_set$DamagePct[which(this_set$Detail_ID == this_b$Detail_ID)] <- NA
        this_set$DamageDepth[which(this_set$Detail_ID == this_b$Detail_ID)] <- NA
        next

      }



      # Fix extreme values
      dmg_est <- ifelse(flood_depth > 7, max(curv$damage, na.rm=TRUE), dmg_est)
      dmg_est <- ifelse(flood_depth < 0.001, min(curv$damage, na.rm=TRUE), dmg_est)

      if(is.na(dmg_est) | dmg_est < 0 | dmg_est > 1){
        stop("bad damage estimate")
      }

      # Calculate damage cost
      damage_cost <- this_b$VAL*dmg_est

      this_set$DamageCost[which(this_set$Detail_ID == this_b$Detail_ID)] <- damage_cost
      this_set$DamagePct[which(this_set$Detail_ID == this_b$Detail_ID)] <- dmg_est
      this_set$DamageDepth[which(this_set$Detail_ID == this_b$Detail_ID)] <- flood_depth

    }


    blds_affected <- nrow(this_set)
    max_depth <- max(this_set$DamageDepth, na.rm=TRUE)
    median_depth <- stats::median(this_set$DamageDepth, na.rm=TRUE)
    median_cost <- stats::median(this_set$DamageCost, na.rm=TRUE)
    median_damage_pct <- stats::median(this_set$DamagePct, na.rm=TRUE)
    total_cost <- sum(this_set$DamageCost, na.rm=TRUE)

    add_row <- data.frame(
      depth = this_depth,
      n_blds = blds_affected,
      total_cost = total_cost,
      median_depth = median_depth
    )

    output[[i]] <- add_row


  }
  #======================================================


  dd_curves <- do.call("rbind", output)
  dd_curves$cost_millions <- dd_curves$total_cost/1000000
  dd_curves$depth_datum <- dd_curves$depth

  dd_curves$GROUP <- this_group
  all_g[[g]] <- dd_curves

  }

  dat <- do.call("rbind", all_g)






  #=====================================================
  # Prep OUTPUT FOLDER
  #=====================================================

  # Create Output Directory
  substrRight <- function(x, n){
    substr(x, nchar(x)-n+1, nchar(x))
  }
  rR  <- substrRight(dir_output, 1)
  iR <- ifelse(rR == "/", "", "/")

  path_output <- paste0(dir_output, iR, simulation_name, "/")

  dir.create(path_output)



  #=====================================================
  # Export Plots
  #=====================================================
  ddids <- unique(Bldgs$DDID)
  dir.create(paste0(path_output, "depth-damage-curves"))


  for(i in 1:length(ddids)){
    tid <- ddids[i]
    td <- HAZUS[which(HAZUS$DDID == tid),]

    grDevices::png(filename = paste0(path_output, "depth-damage-curves/",tid,".png"),
                   width=5,
                   height=5,
                   units="in",
                   res=300)
    graphics::plot(td$depth_m, td$damage*100, type="b",
    xlab="Flood Water Depth (m)",
    ylab="Structure Damage Percent (%)",
    main = paste0("Curve ID: ", tid))
    grDevices::dev.off()

  }





  #=====================================================
  # Export Plots
  #=====================================================
  ugroup <- unique(dat$GROUP)
  dir.create(paste0(path_output, "flooding"))

  for(i in 1:length(ugroup)){

    tid <- ugroup[i]
    td <- dat[which(dat$GROUP == tid),]

    grDevices::png(filename = paste0(path_output, "flooding/",tid,".png"),
                   width=5,
                   height=8,
                   units="in",
                   res=300)
    graphics::par(mfrow=c(3,1))

    graphics::plot(td$depth, td$n_blds, type="l",
         xlab="Flood Water Depth (m CD)",
         ylab="Number of Structures",
         main = paste0("", tid))

    graphics::plot(td$depth, td$median_depth, type="l",
         xlab="Flood Water Depth (m CD)",
         ylab="Median Flood Depth",
         main = paste0("", tid))

    graphics::plot(td$depth, td$cost_millions, type="l",
         xlab="Flood Water Depth (m CD)",
         ylab="Cost ($ Millions)",
         main = paste0("", tid))


    grDevices::dev.off()

  }



  utils::write.csv(dat, row.names = FALSE,
                   paste0(path_output, "flood-study-data.csv")
                   )








}