GFDRR / CCDR-tools

Geoanalytics for climate and disaster risk screening
https://gfdrr.github.io/CCDR-tools/
12 stars 8 forks source link

Plot Exceedance Frequency Curve #2

Open matamadio opened 2 years ago

matamadio commented 2 years ago

Use the results to plot a chart of Annual Exceedence Probability and related EAI at the end, above or below map:

immagine

As in: immagine

matamadio commented 1 year ago

Update on previous with better EAI calculation

immagine

The calculation performed in the table depends on the number (var) of RPs used. In the analysis:

Exp_EAI: Sum(RPi_exp_imp*(1/RPi - 1/RPi+1)

e.g. with RP selection [5, 10, 25, 50]:
RPi = 5
RPi+1 = 10

See example for Exp=Pop

  RP5_pop_imp RP10_pop_imp RP20_pop_imp RP50_pop_imp RP75_pop_imp RP100_pop_imp RP200_pop_imp RP250_pop_imp RP500_pop_imp RP1000_pop_imp POP_EAI
Impact 0.87 65.575 67.475 68.881 69.497 70.07 70.727 71.058 72.037 72.869  
EAI 0.09 3.28 2.02 0.46 0.23 0.35 0.07 0.14 0.07 0.07 6.79
Calc: B2*(1/5-1/10) C2*(1/10-1/20) D2*(1/20-1/50) E2*(1/50-1/75) F2*(1/75-1/100) G2*(1/100-1/200) H2*(1/200-1/250) I2*(1/250-1/500) J2*(1/500-1/1000) K2*(1/1000)  
                       

For each ADM, we are interested to keep in the final table only the columns:

Note that the EAI value changes with the number of RP selected, where the most reliable is 10; a selection of only 3 RP leads to results with a ratio of 56% compared to 10RP EAI in this example over DOM. immagine

artessen commented 1 year ago

Here, a first version of a Python notbook to create the plot function for EAI using matplotlib and seaborn Sample file to run the code here: NGA_RSKv2.csv

# Importing the required libraries
import datetime, os, shutil, gc, socket
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter, StrMethodFormatter
import seaborn as sns

# Reading the data for plotting
dfData = pd.read_csv('NGA_RSKv2.csv') 

# Defining the title and sub-title
title    = 'Flood x Population EAI'
subtitle = 'Exceedance frequency curve'

# Defining the x- and y-axis data and text
x  = dfData[['Freq']].values.squeeze()
y  = dfData[['Pop_impact']].values.squeeze()
xlbl = 'Frequency'
ylbl = 'Impact (population)'

# Defining if plotting total EAI
txtTotal = True
xpos = 0.70
totEAI = np.sum(dfData[['Pop_EAI']].values.squeeze())

# Defining the plot style
sns.set_style('whitegrid')
blue, = sns.color_palette('muted', 1)
xFreqAsRP = True
yOnRight = False

# Getting the labels for RPs
rp_lbl = dfData[['RP']].values.squeeze().tolist()
rp_lbl = ['RP '+str(item) for item in rp_lbl]
x_offset = -0.003
y_offset = 20000
if yOnRight: x_offset = -1*x_offset

# Creating the plot
fig, ax = plt.subplots()
ax.plot(x, y, color=blue, lw=1, marker="o")
if xFreqAsRP:
    for i in range(len(rp_lbl)):
        ax.text(x[i]+x_offset, y[i]+y_offset, rp_lbl[i], color='#4D4D4D',
                bbox=dict(facecolor='#FAFAFA', edgecolor='#4D4D4D', boxstyle='round,pad=0.1', alpha=0.4))
ax.fill_between(x, 0, y, alpha=.3)
ax.set(xlim=(0, max(x)), ylim=(0, max(y)*1.05), xticks=np.linspace(0, max(x), 6))
ax.ticklabel_format(style='plain', axis='y')
ax.yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}'))
ax.set_xlabel(xlbl)
ax.set_ylabel(ylbl)
ax.set_title(str(title+'\n'+subtitle))
if yOnRight:
    ax.invert_xaxis()
    ax.yaxis.tick_right()
    ax.yaxis.set_label_position('right')
if txtTotal:
    vtext = 'EAI = '+'{:,}'.format(totEAI)
    if yOnRight: xpos = 0.04
    ax.text(xpos, 0.9, vtext, fontsize=14, color='black', ha='left', va='bottom', transform=plt.gca().transAxes, 
            bbox=dict(facecolor='#FAFAFA', edgecolor='#4D4D4D', boxstyle='square,pad=0.4', alpha=0.6));

# Show the graph
plt.show()
matamadio commented 1 year ago

Great!

immagine

Added it here: https://github.com/GFDRR/CCDR-tools/blob/main/Top-down/plots/Exceedance%20Frequency%20Curve.py

Would be nice to integrate directly in the script or results dasbhoard.

matamadio commented 1 year ago

@artessen please check latest EAI calculation in the code; in the KHM output, it seems to be simply sum(RPi_imp/RPi) as in the first post; while it should use exceedence frequency instead as explained in updated approach. I.e. RP5_EAI = RP5_IMP * (1/5-1/10)

matamadio commented 1 year ago

To add script in the repo, and possibly integrate in parallel code output @artessen