<jupyter_start><jupyter_text>A Beginner’s Guide to Carry out Extreme Value Analysis (4) - DDF Intensity Duration Frequency (IDF) or Depth Duration Frequency (DDF) curves are among the most common tools used in water resources management. They are derived from historical rainfall records under the assumption of stationarity. Rainfall IDF/DDF curve information describes the frequency of extreme rainfall events for a variety of durations and intensities.The following presents a primary procedure to create DDF curves.  1. Load precipitation dataThis is a hourly data precipitation data. The data have been convert to annual maxima time series for the durations of 1h, 3h, 6h, 12h, and 24h.<jupyter_code>import os
import math
from pandas import *
import pandas as pd
from datetime import datetime

%matplotlib inline 
import matplotlib.pyplot as plt
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 9
data = pd.read_csv('PluviometriaPaperopoli_4.txt')
data=data.set_index('anno')
data.head()<jupyter_output><empty_output><jupyter_text>In the data, the 1h, 3h, 6h, 12h and 24h are the durations, which is also the accumulatin time for precipitation.  Have a simple check for each duration.<jupyter_code>ax = data.plot(title="Annual Maximum Precipitation (mm)", style="o-")
ax.set_xlabel("Year")<jupyter_output><empty_output><jupyter_text>2. Carry out extreme value analysis for each durationThis is a primary step that has been shown before (part 2). For the purpose of illustration, the GUM parameters are used for further analysis.<jupyter_code>import lmoments
from scipy.stats import genextreme,gumbel_r
from numpy import linspace
from statsmodels.distributions.empirical_distribution import ECDF
from collections import OrderedDict<jupyter_output><empty_output><jupyter_text>Fit Gumble distribution for each duration<jupyter_code>gumfit = OrderedDict()
for name in data.columns.values.tolist():     
    dx = (data[name]).dropna()    
    gumfit[name] = lmoments.pelgum(lmoments.samlmu(dx) )<jupyter_output>lmoments\lmoments.py:235: RuntimeWarning: overflow encountered in double_scalars
  obj+obj, obj-obj, obj*obj, obj**obj, obj/obj<jupyter_text>Sample each time series<jupyter_code>samples = np.linspace(data.min()[0],data.max()[4],100)<jupyter_output><empty_output><jupyter_text>Convert samples into probabilities<jupyter_code>samprobs = pd.DataFrame([lmoments.cdfgum(samples, gumfit["1h"]),
                         lmoments.cdfgum(samples, gumfit["3h"]),
                         lmoments.cdfgum(samples, gumfit["6h"]),
                         lmoments.cdfgum(samples, gumfit["12h"]),
                         lmoments.cdfgum(samples, gumfit["24h"])]).T
samprobs.index = samples
samprobs.columns=["1h","3h","6h","12h","24"]<jupyter_output>lmoments\_otherfunct.py:9: RuntimeWarning: overflow encountered in power
  obj+obj, obj-obj, obj*obj, obj**obj, obj/obj<jupyter_text>Compare with observations<jupyter_code>rcParams['figure.figsize'] = 9, 6

plt.rc('xtick', labelsize=20) #Quest  comandi globali
plt.rc('ytick', labelsize=20) #Questi sono comandi globali
ax = samprobs.plot(color=["red","blue","cyan","green","black"])

ecdf1h = ECDF(data["1h"])
ax.plot(data["1h"],ecdf1h(data["1h"]),'o',c="red",label="1h")

data3h = data["3h"].dropna()
ecdf3h = ECDF(data3h)
ax.plot(data3h,ecdf3h(data3h),"o",c="blue",label="3h")

data6h = data["6h"].dropna()
ecdf6h = ECDF(data6h)
ax.plot(data6h,ecdf6h(data6h),"o",c="cyan")

data12h=data["12h"].dropna()
ecdf12h = ECDF(data12h)
ax.plot(data12h,ecdf12h(data12h),"o",c="green")

data24h=data["24h"].dropna()
ecdf24h = ECDF(data24h)

ax.plot(data24h,ecdf24h(data24h),"o",c="black")
ax.set_title('Lssp 1h')
ax.set_xlabel('Rainfall (mm)')
ax.set_ylabel('P[H<h]')<jupyter_output><empty_output><jupyter_text>3. Fit extreme precipitation into theorectical DDF curve for each return year  Setup some return years (e.g., 2, 10, 20, 50, 100, 200)<jupyter_code>q2   = 1-1/2.
q10  = 1-1/10.
q20  = 1-1/20.
q50  = 1-1/50.
q100 = 1-1/100.
q200 = 1-1/200.

q = [q2,q10,q20,q50, q100, q200]<jupyter_output><empty_output><jupyter_text>Calculate extreme values for return years and durations<jupyter_code>samprobs = OrderedDict()

for name in gumfit.keys():
    samprobs[name] = lmoments.quagum(q, gumfit[name])    

pts = pd.DataFrame.from_dict(samprobs, orient='index')
pts.columns=[2,10,20,50, 100, 200]
pts<jupyter_output><empty_output><jupyter_text>convert index into numbers as axis of plotting<jupyter_code>pts.index = [1,3,6,12,24]
pts
pts.plot(style="o-")
plt.xlim([0,30])
plt.ylim([30,200])
ag = pts.plot(style="o-")
ag.set_yscale('log')
ag.set_xscale('log')
plt.xlim([0,30])
plt.ylim([30,200])
plt.xlim([0.7,30])
plt.ylim([30,200])<jupyter_output><empty_output><jupyter_text>The simplest DDF curve formula is like$$ h(t_p, Tr) = a(T_r) t_p^n$$$$\log h(t_p, Tr) = \log a(T_r) + n \log t_p $$  Fit DDF curves based the above formulanumpy.polyfit is used to fit the above curve for each return year in 10, 20 and 100 years.<jupyter_code>ddfparams = OrderedDict()
for name in pts.columns.values.tolist(): 
    ddfparams[name] = np.polyfit(np.log(pts[name].index),np.log(pts[name]),1)
    ddfparams[name][1] = np.exp( ddfparams[name][1])

fnl = pd.DataFrame.from_dict(ddfparams, orient='index')
fnl.columns=["n","a"]
fnl<jupyter_output><empty_output><jupyter_text>Setup some durations for interpolation and plotting (<=30 hours)<jupyter_code>samdurations = np.linspace(0.1,30,100)
samdurations[-5:]<jupyter_output><empty_output><jupyter_text>Use the ddf regression parameters to calculate ddf values<jupyter_code>def regddf(duration, intercepta, expn):
    return intercepta*duration**expn
fnlt = fnl.T

ddfvalues = OrderedDict()
for name in pts.columns.values.tolist():     
    ddfvalues[name] = regddf(samdurations,fnlt[name]["a"],fnlt[name]["n"])
    
inh = pd.DataFrame.from_dict(ddfvalues, orient='index').T.set_index(samdurations)<jupyter_output><empty_output><jupyter_text>set a understandable column names<jupyter_code>inh.columns = ["Tr = 2","Tr = 10","Tr = 20","Tr = 50","Tr = 100", "Tr = 200"]<jupyter_output><empty_output><jupyter_text>Have a visualization check<jupyter_code>ag = inh.plot()

ag = plt.plot(pts[2],  "o", color="blue")
ag = plt.plot(pts[10], "o", color="orange")
ag = plt.plot(pts[20], "o", color="green")
ag = plt.plot(pts[100],"o", color="purple")
ag = plt.plot(pts[200],"o", color="brown")

plt.xlim([0.7,30])
plt.ylim([50,200])

plt.xlabel("Duration (Hour)")
plt.ylabel("Depth of Precipitation (mm)")
ag = inh.plot()
ag.set_yscale('log')
ag.set_xscale('log')

ag = plt.plot(pts[2],  "o", color="blue")
ag = plt.plot(pts[10], "o", color="orange")
ag = plt.plot(pts[20], "o", color="green")
ag = plt.plot(pts[100],"o", color="purple")
ag = plt.plot(pts[200],"o", color="brown")

plt.xlim([0.7,30])
plt.ylim([50,200])

plt.xlabel("Duration (Hour) in Log Transformation")
plt.ylabel("Depth of Precipitation (mm)")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15404586758978797, 'fraction_numerical': 0.04096350786095485, 'mean_word_length': 6.1410916580844495, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 59, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '248995', 'license': 'MIT', 'n_tokens_mistral': 2840, 'n_tokens_neox': 2613, 'n_words': 558, 'path': 'extremeDDF.ipynb', 'repo_name': 'royalosyin/A-Beginner-s-Guide-to-Carry-out-Extreme-Value-Analysis-with-Codes-in-Python-', 'size': '6933'}