Lab 3: Depth-Damage Models

DataFrames and Distributions

Lab
Author

CEVE 421/521

Published

Fri., Jan. 26

Overview

Today, we’re going to be working with depth-damage functions. This will give us practice:

  • working with and manipulating tabular data
  • writing functions

In addition, the depth-damage function you choose / build will be a building block for your final project.

Setup

As before:

  1. Clone the repository for this lab to your computer and open it in VS Code.
  2. In the Julia REPL, activate and then instantiate the project environment.
  3. Check that you can preview the project by running quarto preview template.qmd in the terminal (not Julia REPL). If that doesn’t work, open the Julia REPL, enter package mode with ], and run build IJulia.
  4. If that doesn’t work, ask for help! The way VS Code looks for Python on your computer can be weird and counterintuitive.

Load packages

As usual, we load all required packages at the top of the notebook, in one place.

using CSV
using DataFrames
using DataFramesMeta
using Distributions
using Interpolations
using Plots
using StatsPlots
using Unitful

Plots.default(; margin=6Plots.mm)
1
This updates the default margin in our plots so that axis labels don’t get cut off.

Depth-damage data

Today we’re going to work with deterministic depth-damage functions from the HAZUS model developed by the US Army Corps of Engineers. Please see the data source for more information and links. We’re going to work with the haz_fl_dept file today. We can read it in as before:

haz_fl_dept = CSV.read("data/haz_fl_dept.csv", DataFrame)
first(haz_fl_dept, 3)
1
This lets us view the first three rows
3×52 DataFrame
Row Column1 Occupancy DmgFnId Source Description ft04m ft03m ft02m ft01m ft00 ft01 ft02 ft03 ft04 ft05 ft06 ft07 ft08 ft09 ft10 ft11 ft12 ft13 ft14 ft15 ft16 ft17 ft18 ft19 ft20 ft21 ft22 ft23 ft24 Comment ft00_5 ft01_5 ft02_5 ft03_5 ft04_5 ft05_5 ft06_5 ft07_5 ft08_5 ft09_5 ft10_5 ft11_5 ft12_5 ft13_5 Source_Table Occupy_Class Cover_Class
Int64 String7 Int64 String31 String? String3 String3 String3 String3 Int64 Int64 Int64 Int64 Int64 Int64 Int64 Int64 Int64 Int64 Int64 String3 String3 String3 String3 String3 String3 String3 String3 String3 String3 String3 String3 String3 String3 String? String3 String3 String3 String3 String3 String3 String3 String3 String3 String3 String3 String3 String3 String3 String31 String15 String15
1 1 RES1 183 USACE - Wilmington two story, Pile foundation, structure 3 3 3 4 7 9 15 23 25 29 31 32 34 35 37 40 44 48 51 55 57 60 62 65 68 78 89 100 100 missing NA NA NA NA NA NA NA NA NA NA NA NA NA NA flBldgStructDmgFn RES Bldg
2 2 RES1 184 USACE - Wilmington two story, Structure 0 0 1 1 7 9 15 18 22 25 27 29 32 34 36 40 43 47 50 54 57 60 62 65 68 71 74 77 77 missing NA NA NA NA NA NA NA NA NA NA NA NA NA NA flBldgStructDmgFn RES Bldg
3 3 RES1 105 FIA one floor, no basement, Structure, A-Zone 0 0 0 0 18 22 25 28 30 31 40 43 43 45 46 47 47 49 50 50 50 51 51 52 52 53 53 54 54 missing NA NA NA NA NA NA NA NA NA NA NA NA NA NA flBldgStructDmgFn RES Bldg

We can see that there are a lot of columns. Each depth-damage curve is a row, and each depth (or metadata) is a column.

There’s a lot to explore with these dataframes. One thing we can do in VS Code is to click on the Julia tab in VS Code (three circles on the far left of your window). This will show you a list of all the variables in your workspace. Click on haz_fl_dept and then click on the “preview” icon that pops up (looks like a newspaper emoji). This lets you interactively explore the DataFrame.

Parsing

We’d like to be able to use the depth-damage functions in this file. However, the depths are stored in a somewhat annoying format (e.g., “ft04m” means -4 feet). To make life simple, I’ve created some functionality in the depthdamage.jl file that you can use. We can load it as follows:

include("depthdamage.jl")
DepthDamageData

The main thing that we’ll use is called DepthDamageData. This is a data structure or type that stores the depth-damage data, as well as any relevant metadata. If you’ve created a class in a language like C++ or Python, it’s the same idea. I’ve also defined a constructor that takes in the row of a DataFrame and creates a DepthDamageData object, to make life easy.

I’ll show you how to do this for an illustrative depth-damage function from the New Orleans USACE.

demo_row = @rsubset(
    haz_fl_dept, :Description == "one story, Contents, fresh water, short duration"
)[
    1, :,
]
dd = DepthDamageData(demo_row)
DepthDamageData(Quantity{Float64, 𝐋, Unitful.FreeUnits{(ft,), 𝐋, nothing}}[-4.0 ft, -3.0 ft, -2.0 ft, -1.0 ft, 0.0 ft, 1.0 ft, 2.0 ft, 3.0 ft, 4.0 ft, 5.0 ft  …  15.0 ft, 16.0 ft, 17.0 ft, 18.0 ft, 19.0 ft, 20.0 ft, 21.0 ft, 22.0 ft, 23.0 ft, 24.0 ft], [0.0, 0.0, 0.0, 0.0, 0.0, 42.0, 63.0, 82.0, 85.0, 91.0  …  91.0, 91.0, 91.0, 91.0, 91.0, 91.0, 91.0, 91.0, 91.0, 91.0], String7("RES1"), "57", String31("USACE - New Orleans"), "one story, Contents, fresh water, short duration", "missing")

This prints out a bunch of data. We can see that it has the following fields, which should broadly match with our DataFrame:

fieldnames(typeof(dd))
(:depths, :damages, :occupancy, :dmg_fn_id, :source, :description, :comment)

Plotting

Now that we’ve created a DepthDamageData object, we can plot it. When we plot things with units, the Unitful package (as long as we are using it) knows how to handle them.

scatter(
    dd.depths,
    dd.damages;
    xlabel="Flood Depth at House",
    ylabel="Damage (%)",
    label="$(dd.description) ($(dd.source))",
    legend=:bottomright,
    size=(700, 500),
)

Interpolating

This is great. However, what if we want to estimate damage between the points? We need a way to interpolate. We can do this using the Interpolations package!

itp = let
    depth_ft = ustrip.(u"ft", dd.depths)
    damage_frac = dd.damages
    Interpolations.LinearInterpolation(
        depth_ft,
        damage_frac;
        extrapolation_bc=Interpolations.Flat(),
    )
end
1
I really like these let...end blocks and use them quite a bit. The main thing to know is that all the variables defined inside the let block are only available inside the let block. Once we get to the end of the block, they vanish! This keeps us from defining tons of variables that get in each others’ way.
2
The Interpolations package doesn’t take units on its input, so we convert the input (which can be of any length unit) to feet before passing it in. If our depths are in meters or millimeters, it won’t be a problem – the ustrip function will convert to feet and then turn them into scalars.
3
Interpolations requires us to specify how to extrapolate. We choose Flat(), meaning that anything below the lowest value in the table will be assumed to have the same damage as the lowest value in the table and anything above the highest value in the table will be assumed to have the same damage as the highest value in the table.

Now we can use this interpolation function to estimate damage at any depth.

let
    dmg_fn(x) = itp(ustrip.(u"ft", x))
    dmg_fn.([3.1u"ft", 2.2u"m", 91.4u"inch"])
end
1
Convert the input to feet
2
Estimate damage at 3.1 feet, 2.2 meters, and 91.4 inches 
3-element Vector{Float64}:
 82.30000000000001
 91.0
 91.0

Packaging

To make life simple, we can define a function that takes in some depths and some damages and returns a function that can be used to estimate damage at any depth.

function get_depth_damage_function(
    depth_train::Vector{<:T}, dmg_train::Vector{<:AbstractFloat}
) where {T<:Unitful.Length}

    # interpolate
    depth_ft = ustrip.(u"ft", depth_train)
    interp_fn = Interpolations.LinearInterpolation(
        depth_ft,
        dmg_train;
        extrapolation_bc=Interpolations.Flat(),
    )

    damage_fn = function (depth::T2) where {T2<:Unitful.Length}
        return interp_fn(ustrip.(u"ft", depth))
    end
    return damage_fn
end
1
The Interpolations package doesn’t take units on its input, so we convert the input (which can be of any length) to feet before passing it in. If our depths are in meters or millimeters, it won’t be a problem – the ustrip function will convert to feet and then turn them into scalars.
2
Interpolations requires us to specify how to extrapolate. We choose Flat(), meaning that anything below the lowest value in the table will be assumed to have the same damage as the lowest value in the table and anything above the highest value in the table will be assumed to have the same damage as the highest value in the table.
3
This is a bit confusing. We are defining a function, inside of a function.
4
We return the function that we just defined. So when we call this function, we get a function – we in turn need to call that function on something else. 
get_depth_damage_function (generic function with 1 method)
damage_fn = get_depth_damage_function(dd.depths, dd.damages)
#16 (generic function with 1 method)

Now damage_fn is a function. It takes in a depth, with some type of length unit defined using Unitful, and returns the damage in percent. We can use this to plot a depth-damage curve:

p = let
    depths = uconvert.(u"ft", (-7.0u"ft"):(1.0u"inch"):(30.0u"ft"))
    damages = damage_fn.(depths)
    scatter(
        depths,
        damages;
        xlabel="Flood Depth",
        ylabel="Damage (%)",
        label="$(dd.description) ($(dd.source))",
        legend=:bottomright,
        size=(800, 400),
        linewidth=2,
    )
end
p
1
We create a vector of depths from -7 feet to 30 feet, in 1 inch increments. We use uconvert to convert the units to feet (by default, Unitful converts to meters when we add together length units).
2
Our damage_fn is defined to take in a single scalar. To make predictions about a Vector of depths, we use . to broadcast the function over the vector.

Of course, if we use plot instead of scatter, then we get a line plot which is automatically smooth.

Expected damages

Now that we have a depth-damage function, we can combine it with a probability distribution of hazard to assess the annual expected damages. First, we need to come up with a distribution of hazard! We’re not going to go into extreme value statistics today. Instead, we’re going to consider a simple distribution that is often used to model extreme events: the generalized extreme value distribution. We’ll consider hypotheticalparameter values as an illustration. If you want to adjust them, go for it! This is the distribution of the maximum water level in a given year, in feet, at our gauge, for a single year (i.e., we’re not considering sea-level rise).

gauge_dist = GeneralizedExtremeValue(5, 1.5, 0.1)
GeneralizedExtremeValue{Float64}(μ=5.0, σ=1.5, ξ=0.1)

We can see some quantiles of the distribution:

quantile.(gauge_dist, [0.5, 0.9, 0.99, 0.999])
4-element Vector{Float64}:
  5.559968481853559
  8.785530774057731
 13.761464356944844
 19.927437691344867

and we can plot it using StatsPlots:

p1 = plot(
    gauge_dist;
    label="Gauge Distribution",
    xlabel="Water Level (ft)",
    ylabel="Probability Density",
    legend=:topright,
    linewidth=2,
)

Our building might be above the gauge. We can correct for this by adding an offset to the location parameter of the gauge. Again, we are taking hypothetical values here!

offset = 2.5 # house is 2.5 feet above gauge
house_dist = GeneralizedExtremeValue(gauge_dist.μ - offset, gauge_dist.σ, gauge_dist.ξ)
GeneralizedExtremeValue{Float64}(μ=2.5, σ=1.5, ξ=0.1)

We can plot this

plot!(p1, house_dist; label="House Distribution", linewidth=2)
1
This adds to the existing plot

We can see that the distribution of hazard at the house is shifted left. That is, the house is less likely to experience extreme water levels than the gauge. This makes sense – it’s higher up!

Important: this is a very simple approach and only makes sense if the house is near the gauge.

Now that we have a distribution of hazard, we can combine it with our depth-damage function to estimate the expected annual damages. A very simple Monte Carlo algorithm is

  1. Sample N values from the distribution of hazard
  2. For each value, estimate the damage using the depth-damage function
  3. Average the damages

Instructions

Edit the template.qmd file to complete the following tasks.

  1. Pick a site. For your final project, we will develop a decision support tool around the question of whether a building subject to coastal flooding should be elevated. For future analyses to run smoothly, you should pick a building that is near a water gauge with a long record – I suggest Sewells Point, VA or Galveston Pier 21, TX, but feel free to find another gauge. You can change sites later, but you’ll need to re-do this analysis. Once you have chosen your site, make sure there’s a long record of extreme water levels (not many gaps) by clicking Tides/Water Levels and then Extreme Water Levels. You should get something that looks like this. Indicate which site you have chosen.
  2. Pick a building. Find the gauge on Google Maps, and then find a building that is near the gauge. Define what you mean by “near”, but the building should not be too far inland. Indicate which building you have chosen.
  3. Find the building’s elevation. Find the elevation of the building. You can do this with USGS data following these instructions. Record your estimate.
  4. Find the building’s depth-damage data. Find a depth-damage function for your building. Use one of the depth-damage functions in haz_fl_dept.csv. Read the documemtation and make sure you understand what it means. Explain why this is an appropriate depth-damage function to use.
  5. Build the depth-damage curve. Use the tools we have built in this workshop to build a depth-damage curve for your building. Plot the curve in 1 inch increments from -10 to 30 feet, as above. Explain what it means.
  6. Implement the Monte Carlo algorithm described above to estimate the expected annual damages. Use 1,000,000 samples. Explain what this means.

Finally, add any plots or discussion that you think are relevant! For example, consider looking at plausible alternative depth-damage functions and the sensitivity of your results to the choice of depth-damage function.