RP Methods - Tabular Admin 2 Pipeline

FloodScan HDX

Work in Progress

1 Intro

1.1 Key Recommendations

  1. Internal pipeline should include RP calculations made using Log-Pearson Type III (LP3).The value should be kept as the raw calculated RP value. Additionally, we may consider a. including columns to hold RPs calculated from Gumbel, GEV, and Empirical Distributions or b. have functions at our disposal to calculate these as needed.
  2. HDX pipeline should use RP value from LP3 calculation. We will provide raw figures from 1-100 and classify values above the lower an upper bounds as <=1 & >100, respectively.
  3. Team should revisit RP calculation for data sets with small historical records (i.e IRI seasonal forecast)

1.2 Takeaways

  • The resulting differences in RP calculations from different methods on FloodScan data is smaller at lower RP values (i.e < 5 year RP). The differences get larger at higher return periods, and the resulting values are sensitive to the shape of the distributions.
  • In some cases the distribution-based methods are either required or more appropriate. These include:
    • calculating return period values with low sample size. The closer the return period of interest gets to the length of the historical record the more important this becomes.
    • calculating return values for return periods > historical record. This cannot be done empirically.
    • providing error terms/confidence metrics.
    • Converting historical “raw” values to return periods. This requires a method that can appropriately extrapolate return periods beyond the historical record.
  • The classic empirical method can be thought of as a percentile algorithm. There are several percentile calculation algorithms that give slightly different results. When using Empirical methods, team should standardize.

1.3 Background

As part of our AER FloodScan-HDX pipeline we will be providing a tabular (excel/csv) data set with Admin 2 level zonal statistics of simplified AER FloodScan raster data. SFED raster pixel values represent estimated flood fraction (%). Therefore, our principle zonal statistic calculated will most likely be mean Standard Flood Event Depiction (SFED) calculated at the Admin 2 level to provide an average flood fraction for the Admin 2 zone.

Data sharing limitations prevent us from sharing the entire historical record of data. Therefore, we would like to provide some key metrics so that the user can better contextualize the values with respect to the historical record/distribution. In addition to a time series of mean SFED values for each day & admin unit, we have proposed providing:

  1. Return period estimations per row.
  2. Historical baseline estimates per row.

Below is a blank template of what is proposed for the tabular data set.

This document explores methods to calculate return periods and return values and was originally based upon this proposal.

1.3.1 RP Methods Covered

This document explores various options/methods for the calculation to provide a basis for a decision. While flood fraction is a unique case, these decisions also relate to general team discussions around standardizing return period calculations where possible.

There are various methods to calculate return period levels which include both empirical and distribution based methods. Here are the methods explored in this document.

  1. Empirical methods
  2. Probability distribution based methods

2 Method Comparison

Here we compare all the methods to eachother. The {extRemes} package is used to model the GEV and Gumbel distribution fits. The LP3 distribution is first fit following the instructions in this video (adapted to R code). We then fit the LP3 distribution to the data using L-moments using the {lmom} R package to compare the LP3 results.

2.1 LP3, Gumbel, GEV, Empirical

2.2 Can Distributions Exceed 100%?

2.3 LP3 vs Empirical

“According to the U.S. Water Advisory Committee on Water Data (1982), the Log-Pearson Type III Distribution is the recommended technique for flood frequency analysis” OSU

2.4 LP3 Method Comparison

Below we compare the “manual” calculation implemented following the instructions available in the Youtube video to an implementation using L-moments where we calculate the LP3 CDF function and then the RP and RP values with help from the {lmom} R package.

The results are very similar and we believe either approach can be justified.

# A tibble: 18 × 2
   adm1_en           max
   <chr>           <dbl>
 1 Awdal            35.7
 2 Bakool           21.4
 3 Banadir          15.7
 4 Bari             69.6
 5 Bay             164. 
 6 Galgaduud        65.9
 7 Gedo             58.4
 8 Hiraan           31.9
 9 Lower Juba      114. 
10 Lower Shabelle  311. 
11 Middle Juba     149. 
12 Middle Shabelle  43.1
13 Mudug            30.2
14 Nugaal           45.5
15 Sanaag          105. 
16 Sool             88.7
17 Togdheer         61.8
18 Woqooyi Galbeed  21.2
# A tibble: 3 × 3
# Groups:   adm1_en [1]
  adm1_en adm1_pcode  params
  <chr>   <chr>        <dbl>
1 Bari    SO16       0.00395
2 Bari    SO16       0.00339
3 Bari    SO16       1.84

3 Example of Admin 1 Stats

To give an example of what the admin tabular data would look like, below we show a subset of one administrative zones tabular data.

Most The resulting historical data set is strongly skewed towards 1. So much that you cannot really see the distribution above 1 even when log scaled. Therefore, below we show a distribution plot of all historical RP values above 1.5 using the different methods.

4 Discussion

TBD

  • “According to the U.S. Water Advisory Committee on Water Data (1982), the Log-Pearson Type III Distribution is the recommended technique for flood frequency analysis” OSU
  • Can make a final decision on LP3 method (moments vs conventional)
  • Should decide on imputation method for when 0 is present in the data. It’s not that common because we take the annual maxima first, but will occur. In current analysis we are just imputing with second the lowest value.

5 Appendix

5.1 Formulas

5.1.1 Gumbel Formula

CDF

\[ f(x) = \frac{1}{\beta} \exp\left(- \left( \frac{x - \mu}{\beta} \right) - \exp\left(- \frac{x - \mu}{\beta} \right) \right) \] PDF

\[ F(x) = \exp\left( - \exp\left( - \frac{x - \mu}{\beta} \right) \right) \]

5.1.2 GEV Formula

CDF

\[ F(x) = \exp\left( - \left( 1 + \xi \frac{x - \mu}{\beta} \right)^{-1/\xi} \right), \quad 1 + \xi \frac{x - \mu}{\beta} > 0 \]

PDF

\[ f(x) = \frac{1}{\beta} \left( 1 + \xi \frac{x - \mu}{\beta} \right)^{-1 - 1/\xi} \exp\left( - \left( 1 + \xi \frac{x - \mu}{\beta} \right)^{-1/\xi} \right) \]

5.1.3 LP3 RP Formula

\[ X_T = 10^{\left( \mu + \frac{2}{g} \left( \left( \frac{Z_T - \frac{g}{6}}{\frac{g}{6}} + 1 \right)^3 - 1 \right) \cdot \sigma \right)} \]

Where:

  • \(X_T\) is the quantile for the return period \(T\),
  • \(\mu\) is the mean of the log-transformed data,
  • \(\sigma\) is the standard deviation of the log-transformed data,
  • \(g\) is the skewness of the log-transformed data,
  • \(Z_T\) is the standard normal quantile corresponding to the exceedance probability \(P_T = \frac{1}{T}\).

5.2 Filtering out low-duration peaks

TBD

There was a suggestion to filter out low duration peaks (< 2 or 3days) before running RP calculation. I started thinking about this, but not sure how to implement it in a way that would be generalization or how necessary it is. Some thoughts:

It’s the max values per year/admin that we are concerned with as those are the input for the RP calculations, but we don’t currently have an algorithm picked out that would generalize well to this across thousands of admin 2 units. I thought through some possibilities below, and so far think it’s overly complicated and would cause problems with an automated pipeline for a variety of reasons:

  1. We don’t have a threshold that constitutes whether or not a peak is sustained.
  2. It would be a cause for concern if the max peak days per year were false positives. I could imagine false positive noise scattered through the time series, but would think if they were high enough to create a yearly maxima that would indicate some more systemic issue?
  3. Using a n-day average to calculate the max peak per year would dampen the effect of extreme daily value which could be real, especially in flash flood scenarios which FloodScan does seem to pick up well in some contexts? Also it could be confusing if the admin stats are provided at daily time-steps, but the RP calculations is based on n-day average value.

Theoretically we could attempt something like the below (but just writing it out seems problematic):

  1. Take the take an n-day (3 day?) rolling centered mean across the time series smooth the values.
  2. calculate the max daily value per year + admin, max n-day mean value per year + admin
  3. We could then see if the n-day max value aligns with the max 3-day mean value. We could say it does if they are with n (5?) days of each other. If they are not, we could use the next highest daily peak that DOES align with the smoothed peak.
    • Seems problematic as the appropriate n could change based on admin. There will be a large number of admin + year combinations that will be very flat with almost all 0 values.
    • additionally there will be many instances when the same exact max value occurs more than 1x per year. For simple unfiltered RP calculation this is fine, but when trying to consider if either are not a sustained peak it adds a complication.

5.3 RPs & Seasonality - General

It was decided not to build seasonality into the RP calculations and rather handle this component in our baseline calculations which the user can use to better understand seasonal anomalies. Nonetheless, I’ve tried to summarize some of the discussions/thinking below.

We were originally wondering whether it would make sense to incorporate seasonality into RP calculations. Therefore, rather than taking annual maxima for the entire year one could choose the window of interest i.e July-Sep and take the maxima from that range per year. We often do this when designing Anticipatory Action (AA) framework thresholds to assign risk/probabilities for certain windows/seasons. The benefit of adding seasonality to the RP metric in the simplified FloodScan data would be that it would identify out-of-season anomalies (potentially more flash flood type events).

RP typically is used to understand extreme event levels with respect to there probability of occurrence. For example, “First Street’s data suggests that 17.7 million properties nationwide are at risk in a 100-year event.” This means that 17.7 million properties would be at risk if a storm level is as high as one that typically occurs every 100 years and has 1 % chance of occurring any given year.

You will notice that seasonality is not typically built into the example above which is the a pretty standard use-case. So let’s see how the communication/interprtation wouuld change if we play out a hypothetical scenario where we calculate RP from a smaller defined window of 10 days:

Calculating FloodScan RP for today, 17 September: we calculate the RP value based on only historical data from 12 September - 22 September. The “results” is a 1 in 5 year RP level event. The interpretation is: “Sep 17 had a flood fraction level of 1 in 5 years for the dates of 12-22 September.” While I could see this being useful to some niche analyses and technical audiences. It’s not a simple way to communicate and use the concept of RPs to a more general audience.

Therefore, we think it is better to communicate seasonal outliers with season anomalies which can be calculated by the user from baseline_avg value that will be provided.

5.4 Original Questions

  • what range to to calculate the RP value. The first step in calculating RP values is to take the annual maxima of the data. For seasonal analyses we look for this annual maxima for just the season of interest. In this floodscan product we don’t have any inherent seasons of interest, however in the raster product to calculate historical baseline we do apply a smoothing function of +/-n days (likely 10). In theory we could recalculate the RP values based on only that window rather than the full year.
  • Calculations should be based on full year annual maxima
  • Utility of RP column especially since most values are 1. Can we calculte values below 1?
  • Which EVD distribution to use?
  • LP3 as it is a hydrological standard and performs with FloodScan when compared to Gumbel and GEV

The plots below show that different empirical methods do diverge in there estimation of return period levels depending on the characteristics of the data.

One thing that is not immediately evident from the plot is that the “Weibell (manual) based methodology” is nearly equivalent to the quantile() function used for the rest of the methods when the type algorithm is set to 3 which represents the weibell formula which is the default method in minitab software. Below is a table looking at the correlation between Weibell manual and quantile(..., type = 3). I am not totally sure why it is not a perfect correlation, but perhaps small differences in interpolation method and rounding.

5.5 Additional Reading to Check out