Return period estimation RS flood 2024

main.py

# Fitting a Gumbel distribution using scipy's built-in functions for simplicity and stability
from scipy.stats import gumbel_r

# Provided water heights including the latest years up to 2024
water_heights = [
    2.60, 1.48, 0.98, 1.99, 1.45, 1.51, 2.50, 1.53, 2.00, 1.69, 1.52, 1.34, 2.05, 2.13, 1.19,
    2.60, 1.91, 1.78, 0.98, 1.49, 2.21, 1.60, 1.58, 1.68, 1.75, 1.61, 1.31, 2.60, 1.56, 3.20,
    2.05, 2.35, 1.70, 1.84, 1.34, 1.70, 1.64, 3.24, 2.51, 1.43, 1.60, 2.24, 4.75, 2.33, 1.60,
    1.90, 1.26, 1.55, 1.67, 1.68, 1.71, 1.91, 2.10, 2.06, 2.52, 2.91, 1.80, 2.32, 2.08, 2.00,
    1.99, 1.77, 2.16, 1.25, 2.67, 1.73, 2.72, 2.61, 3.13, 1.18, 1.36, 1.71, 1.72, 2.21, 1.93,
    1.48, 1.64, 1.84, 2.13, 1.19, 1.66, 1.58, 1.54, 1.97, 2.32, 2.56, 1.96, 1.73, 2.36, 1.98,
    2.00, 2.22, 1.45, 1.94, 2.07, 1.86, 1.96, 1.62, 1.96, 1.97, 1.46, 1.86, 2.40, 2.46, 1.74,
    1.56, 2.10, 1.38, 2.44, 1.82, 2.23, 1.62, 2.04, 1.66, 2.24, 2.11, 2.94, 5.33, 3.46, 2.65
]


# Fitting Gumbel distribution to the data
gumbel_params = gumbel_r.fit(water_heights)  # Fit the Gumbel distribution
loc_gumbel, scale_gumbel = gumbel_params

# Function to calculate the return period
def return_period_gumbel(value, loc, scale):
    # Calculate the exceedance probability
    exceedance_prob = 1 - gumbel_r.cdf(value, loc=loc, scale=scale)
    # Calculate the return period
    return 1 / exceedance_prob

# Calculate the return period for the extreme event of 5.33 meters
return_period_gumbel_533 = return_period_gumbel(5.33, loc_gumbel, scale_gumbel)

# Bootstrap to estimate confidence intervals
bootstrap_estimates_gumbel = []
bootstrap_samples = 1000
for _ in range(bootstrap_samples):
    sample = np.random.choice(water_heights, size=len(water_heights), replace=True)
    loc_sample, scale_sample = gumbel_r.fit(sample)
    rp_sample = return_period_gumbel(5.33, loc_sample, scale_sample)
    bootstrap_estimates_gumbel.append(rp_sample)

# Compute the 95% confidence interval from the bootstrap results
ci_lower_gumbel = np.percentile(bootstrap_estimates_gumbel, 2.5)
ci_upper_gumbel = np.percentile(bootstrap_estimates_gumbel, 97.5)

(loc_gumbel, scale_gumbel, return_period_gumbel_533, (ci_lower_gumbel, ci_upper_gumbel))