<jupyter_start><jupyter_text>Assignment CIE 5703 - week 6  Import libraries<jupyter_code>import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
%matplotlib inline
plt.style.use('ggplot')<jupyter_output><empty_output><jupyter_text>Rotterdam rain gauge dataset 10 min data from 2003 - 2013 Read in data<jupyter_code>data = pd.read_csv('rotterdam_rg_2003-2014.csv', skipinitialspace=True)<jupyter_output><empty_output><jupyter_text>Convert the dates to a readable format...<jupyter_code>dates = data['Datum']
time = data['Tijd']
dates = dates.map(str)
date_split = dates.str.extract('(.{4})(.{2})(.{2})', expand=True)

time=time.apply(lambda x: '%04d' % x)
time_split = time.str.extract('(.{2})(.{2})', expand=True)

date_split.loc[:,3] = time_split.loc[:,0]
date_split.loc[:,4] = time_split.loc[:,1]

data.loc[:,'dt'] = pd.to_datetime(dict(year=date_split[0], month=date_split[1], day=date_split[2], hour=date_split[3], minute=date_split[4]))
data.index=data['dt']
data.head()<jupyter_output><empty_output><jupyter_text>Plot all data<jupyter_code>plt.plot(data['dt'], data["R10m"])
plt.ylabel('mm/10min')
plt.gcf().autofmt_xdate()<jupyter_output><empty_output><jupyter_text>Resample the 10 min dataset to hourly accumulated data<jupyter_code>data_1h = pd.DataFrame()
data_1h['mean_rain'] = data.R10m.resample('H').mean()
data_1h['accum_rain'] = data.R10m.resample('H').sum()
data_1h.tail()
plt.plot(data_1h["accum_rain"])
plt.ylabel('mm/1h')
plt.gcf().autofmt_xdate()
plt.plot(data_1h["mean_rain"])
plt.ylabel(r'mm/10min ($\varnothing$ of 1h)')
plt.gcf().autofmt_xdate()<jupyter_output><empty_output><jupyter_text>Resample 10 min dataset to 24 accumulated data<jupyter_code>data_24h = pd.DataFrame()
data_24h['mean_rain'] = data.R10m.resample('D').mean()
data_24h['accum_rain'] = data.R10m.resample('D').sum()
plt.plot(data_24h["accum_rain"])
plt.ylabel('mm/24h')
plt.gcf().autofmt_xdate()
plt.plot(data_24h["mean_rain"])
plt.ylabel(r'mm/10min ($\varnothing$ of 24h)')
plt.gcf().autofmt_xdate()<jupyter_output><empty_output><jupyter_text>Select summer and winter months as separate datasets<jupyter_code>data_summer_1h = data_1h.loc[(data_1h.index.month>=4) & (data_1h.index.month<=9)]

mask_start = (data_1h.index.month >= 1) & (data_1h.index.month <= 3)
mask_end = (data_1h.index.month >= 10) & (data_1h.index.month <= 12)
mask = mask_start | mask_end

data_winter_1h = data_1h.loc[mask]
plt.plot(data_summer_1h["accum_rain"])
plt.ylabel('mm/h')
plt.gcf().autofmt_xdate()
plt.plot(data_winter_1h["accum_rain"])
plt.ylabel('mm/h')
plt.gcf().autofmt_xdate()<jupyter_output><empty_output><jupyter_text>Resample 10 min rain data to monthly accumulated data<jupyter_code>data_monthly = pd.DataFrame()
data_monthly['mean_rain'] = data.R10m.resample('M').mean()
data_monthly['accum_rain'] = data.R10m.resample('M').sum()
plt.plot(data_monthly["accum_rain"])
plt.ylabel('mm/month')
plt.gcf().autofmt_xdate()
plt.plot(data_monthly["mean_rain"])
plt.ylabel(r'mm/10min ($\varnothing$ per month)')
plt.gcf().autofmt_xdate()<jupyter_output><empty_output><jupyter_text>Answering the assignments  1. General statistics for 24-hour and 10-min datasets: compute mean, standard deviation, skewness; plot histograms  10 min dataset<jupyter_code>print('Mean: %s' % str(data.R10m.mean()))
print('Std: %s' % str(data.R10m.std()))
print('Skew: %s' % str(data.R10m.skew()))
data.R10m.hist(bins = 100)
plt.xlabel('mm/10min')
plt.gca().set_yscale("log")
cur_data = data.R10m.loc[data.R10m>0]
hist_d = plt.hist(cur_data, bins=100)
plt.xlabel('mm/10min')
plt.gca().set_yscale("log")<jupyter_output><empty_output><jupyter_text>24 h dataset<jupyter_code>print('Mean: %s' % str(data_24h.accum_rain.mean()))
print('Std: %s' % str(data_24h.accum_rain.std()))
print('Skew: %s' % str(data_24h.accum_rain.skew()))
data_24h.accum_rain.hist(bins = 100)
plt.gca().set_yscale("log")
plt.xlabel('mm/24h')<jupyter_output><empty_output><jupyter_text>Ignoring zeros:<jupyter_code>cur_data = data_24h.accum_rain.loc[data_24h.accum_rain>0]
hist_d = plt.hist(cur_data, bins=100)
plt.gca().set_yscale("log")
plt.xlabel('mm/24h')
data_24h.mean_rain.hist(bins = 100)
plt.xlabel(r'mm/10min ($\varnothing$ per 24h)')
plt.gca().set_yscale("log")
selected_monthly_data = data_monthly
#selected_monthly_data = data_monthly[(data_monthly.index >= '2004-01-01')]
selected_monthly_data.head()<jupyter_output><empty_output><jupyter_text>2. a. Analysis of seasonal cycles: create boxplots for monthly totals across all year<jupyter_code>pd.options.mode.chained_assignment = None  # default='warn'
selected_monthly_data['mon'] = selected_monthly_data.index.month
selected_monthly_data['year'] = selected_monthly_data.index.year
selected_monthly_data.boxplot(column=['accum_rain'], by='mon', sym='+')
plt.ylabel('mm/month')<jupyter_output><empty_output><jupyter_text>Or per year:<jupyter_code>selected_monthly_data.boxplot(column=['accum_rain'], by='year', sym='+')
plt.ylabel('mm/month')
plt.gcf().autofmt_xdate()<jupyter_output><empty_output><jupyter_text>2. b. Analysis of diurnal cycles: create boxplots for hourly totals for entire dataseries<jupyter_code>data_1h['hour'] = data_1h.index.hour
data_1h.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>Neglecting events <1mm/h<jupyter_code>cur_df = data_1h.copy()
cur_df.loc[cur_df.accum_rain<1, 'accum_rain'] = np.nan
cur_df.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>Neglecting events < 3mm/h<jupyter_code>cur_df = data_1h.copy()
cur_df.loc[cur_df.accum_rain<3, 'accum_rain'] = np.nan
cur_df.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>2. c. Variation of diurnal cycles with seasons: create boxplots for hourly totals for summer season (April – September) and for winter season (October-March)<jupyter_code>data_summer_1h['hour'] = data_summer_1h.index.hour
data_summer_1h.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>Neglecting events <1mm/h<jupyter_code>cur_df = data_summer_1h.copy()
cur_df.loc[cur_df.accum_rain<1, 'accum_rain'] = np.nan
cur_df.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>Neglecting events <3mm/h<jupyter_code>cur_df = data_summer_1h.copy()
cur_df.loc[cur_df.accum_rain<3, 'accum_rain'] = np.nan
cur_df.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>Repeat with winter<jupyter_code>data_winter_1h['hour'] = data_winter_1h.index.hour
data_winter_1h.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>Neglecting events <1mm/h<jupyter_code>cur_df = data_winter_1h.copy()
cur_df.loc[cur_df.accum_rain<1, 'accum_rain'] = np.nan
cur_df.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>Neglecting events <3mm/h<jupyter_code>cur_df = data_winter_1h.copy()
cur_df.loc[cur_df.accum_rain<3, 'accum_rain'] = np.nan
cur_df.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>2. d. Diurnal cycles of intense storm events: Count nr of exceedances above 10 mm/h threshold for each hour of the day, for entire data series and for summer months only Show rainfall events > 10mm /h over entire 1h accumulated dataset<jupyter_code>rotterdam_1h_exceeds = data_1h.accum_rain[data_1h.accum_rain>5]
y = np.array(rotterdam_1h_exceeds)
N = len(y)
x = range(N)
plt.bar(x, y)
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>Amount of events<jupyter_code>print(len(rotterdam_1h_exceeds))<jupyter_output>183<jupyter_text>Events in summer<jupyter_code>rotterdam_1h_exceeds_summer = data_summer_1h.accum_rain[data_summer_1h.accum_rain>10]

y = np.array(rotterdam_1h_exceeds_summer)
N = len(y)
x = range(N)
plt.bar(x, y)
plt.ylabel('mm/h')<jupyter_output><empty_output><jupyter_text>Amount of events<jupyter_code>print(len(rotterdam_1h_exceeds_summer))<jupyter_output>32<jupyter_text>3. Fit GEV-distribution for POT values in the time series  3. a. Create plots: histogram and GEV fit and interpret<jupyter_code>exceed_hist = plt.hist(rotterdam_1h_exceeds, bins=100)
from scipy.stats import genextreme
x = np.linspace(5, 30, 1000)
y = np.array(rotterdam_1h_exceeds[:])
np.seterr(divide='ignore', invalid='ignore')
genextreme.fit(y)
np.seterr(divide='ignore', invalid='ignore')
pdf = plt.plot(x, genextreme.pdf(x, *genextreme.fit(y)))
pdf_hist = plt.hist(y, bins=50, normed=True, histtype='stepfilled', alpha=0.8)<jupyter_output><empty_output><jupyter_text>3. c. Compute rainfall amounts associated with return periods of 1 year, 10 years and 100 years<jupyter_code>genextreme.ppf((1-1/1), *genextreme.fit(y))
genextreme.ppf((1-1/10), *genextreme.fit(y))
genextreme.ppf((1-1/100), *genextreme.fit(y))<jupyter_output><empty_output><jupyter_text>Lower since we're missing out the heavy rainfall events in 2015 and 2016 (only 10 year dataset)  Update 10.10.2017  Block maxima & GEV<jupyter_code>from scipy.stats import genpareto
temp_monthly = data_1h.groupby(pd.TimeGrouper(freq='M'))
block_max_y = np.array(temp_monthly.accum_rain.max())
print(block_max_y)
print(len(block_max_y))
x = np.linspace(0, 30, 1000)
pdf_bm = plt.plot(x, genextreme.pdf(x, *genextreme.fit(block_max_y)))
pdf_hist_bm = plt.hist(block_max_y, bins=100, normed=True, histtype='stepfilled', alpha=0.8)<jupyter_output><empty_output><jupyter_text>GEV and block maxima of monthly maxima of 1h data<jupyter_code>genextreme.fit(block_max_y)
genextreme.ppf((1-1/10), *genextreme.fit(block_max_y))<jupyter_output><empty_output><jupyter_text>GPD & POT<jupyter_code>pdf_bm = plt.plot(x, genpareto.pdf(x, *genpareto.fit(y)))
pdf_hist_bm = plt.hist(y, bins=100, normed=True, histtype='stepfilled', alpha=0.8)<jupyter_output><empty_output><jupyter_text>GPD and POT of data>10mm/h<jupyter_code>genpareto.fit(y)
genpareto.ppf((1-1/10), *genpareto.fit(y))<jupyter_output><empty_output><jupyter_text>Boxplot of POT values<jupyter_code>event_occurences = pd.DataFrame(rotterdam_1h_exceeds)
event_occurences['hour'] = event_occurences.index.hour
event_occurences.boxplot(column=['accum_rain'], by='hour', sym='+')
event_occurences.hour.value_counts(sort=False)
cur_hist = plt.hist(event_occurences.hour, bins=24, histtype='stepfilled')
plt.xticks(range(24))
plt.xlabel('hour')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16570400530252816, 'fraction_numerical': 0.030868288987785247, 'mean_word_length': 9.039923954372624, 'pattern_counts': {'":': 0, '<': 146, '<?xml version=': 0, '>': 141, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '33', 'dir': 'jupyter-structured-clean-dedup', 'id': '304645', 'license': 'MIT', 'n_tokens_mistral': 4691, 'n_tokens_neox': 4332, 'n_words': 604, 'path': 'week_6_Rotterdam.ipynb', 'repo_name': 'kriete/cie5703_notebooks', 'size': '10561'}