<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Assignment CIE 5703 - week 6
Import Libraries
Step1:  Charlotte rain gauge dataset 15 min data from 2003 - 2014
Show the data source locations. Red dots are available gauges, blue cross denotes the selected station. No specific criteria was chosen to select that specific station.
Plot the positions of gauges. Optional code! Requires additional sources. Not easily copy-paste-able!
Step2:  Read in data
Step3:  Format data to year, month, day, hour, min and rainfall & select only ONE rain gauge
Step4:  Plot rain data as read
Step5:  Replace invalid data with NaNs and plot again
Step6:  Resample the 10-min dataset to 24h accumulated rainfall data
Step7:  Resample 15 min data to 1h accumulated dataset
Step8:  Select only summer months (April - Sept)
Step9:  Select only winter months (Oct - Mar)
Step10:  Resample 15 min dataset to monthly accumulated dataset
Step11:  Answering the assignments
1. General statistics for 24-hour and 15-min datasets
Step12:  Histogram of the data
Step13:  Histogram of the data without zeros
Step14:  24h accumulated dataset
Mean, standard deviation and skewness of 24h accumulated dataset
Step15:  Histogram of the dataset
Step16:  Histogram without zeros
Step17:  2. a. Analysis of seasonal cycles
Step18:  Or on a yearly scale
Step19:  2. b. Analysis of diurnal cycles
Step20:  Neglecting events < 1mm/h
Step21:  Neglecting events < 3mm/h
Step22:  2. c. Variation of diurnal cycles with seasons
Step23:  Neglecting events <1mm/hour
Step24:  Neglecting events <3mm/hour
Step25:  Merge hourly winter data
Step26:  Neglecting events <1mm/h
Step27:  Neglecting events <3mm/h
Step28:  2. d. Diurnal cycles of intense storm events
Step29:  Amount of hourly events
Step30:  10 mm/h events in summer periods
Step31:  Amount of hourly events
Step32:  3. Fit GEV-distribution for POT values in the time series
3. a. Create plots
Step33:  3. c. Compute rainfall amounts associated with return periods of 1 year, 10 years and 100 years
Step34:  Update 10.10.2017
Block maxima & GEV
Step35:  GEV and block maxima of monthly maxima of 1h data
Step36:  POT & GPD
Step37:  GPD and POT of data>10mm/h
Step38:  Boxplot of POT values
Step39:  Number of occurences per hour
<ASSISTANT_TASK:>
Python Code:
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
%matplotlib inline
plt.style.use('ggplot')
Explanation: Assignment CIE 5703 - week 6
Import Libraries
End of explanation
from mpl_toolkits.basemap import Basemap
def get_basemap(_resolution):
    return Basemap(projection='merc', llcrnrlat=25, urcrnrlat=38, llcrnrlon=275, urcrnrlon=285, lat_ts=35.,
                   resolution=_resolution)
positions = pd.read_csv('./Raw_RG_Data/RG_lat_lon.csv', header=None)
positions.columns=['lat', 'lon']
plt.figure(figsize=(24,12))
m = get_basemap('h')
m.drawcoastlines()
m.drawcountries()
m.fillcontinents(color = 'gray')
m.drawmapboundary()
for index, row in positions.iterrows():
    x,y = m(row['lon']+360, row['lat'])
    m.plot(x, y, 'ro', markersize=6)
x,y = m(positions['lon'][0]+360, positions['lat'][0])
m.plot(x, y, 'bx', markersize=6)
m.drawstates()
m.drawrivers()
plt.show()
Explanation: Charlotte rain gauge dataset 15 min data from 2003 - 2014
Show the data source locations. Red dots are available gauges, blue cross denotes the selected station. No specific criteria was chosen to select that specific station.
Plot the positions of gauges. Optional code! Requires additional sources. Not easily copy-paste-able!
End of explanation
charlotte_rainfall = pd.read_csv('./charlotte_rg_2003-2014.csv', header = None)
#charlotte_rainfall = pd.read_csv('./Raw_RG_Data/Charlotte_CRN_gage_2003.csv', header = None)
#for i in range(2004,2014):
#    cur_rainfall = pd.read_csv('./Raw_RG_Data/Charlotte_CRN_gage_%d.csv' % i, header = None)
#    charlotte_rainfall = charlotte_rainfall.append(cur_rainfall, ignore_index=True)
Explanation: Read in data
End of explanation
#charlotte_rainfall = charlotte_rainfall.iloc[:,:6]
charlotte_rainfall.columns = ["year","month","day", "hour", "min", "Rainfall"]
charlotte_rainfall.loc[:,'dt'] = pd.to_datetime(dict(year=charlotte_rainfall['year'], month=charlotte_rainfall['month'], day=charlotte_rainfall['day'], hour=charlotte_rainfall['hour'], minute=charlotte_rainfall['min']))
charlotte_rainfall.index=charlotte_rainfall['dt']
charlotte_rainfall.head()
Explanation: Format data to year, month, day, hour, min and rainfall & select only ONE rain gauge
End of explanation
plt.plot(charlotte_rainfall['dt'], charlotte_rainfall["Rainfall"])
plt.ylabel('mm/15min')
plt.gcf().autofmt_xdate()
Explanation: Plot rain data as read
End of explanation
charlotte_rainfall["Rainfall"] = charlotte_rainfall["Rainfall"].replace(-99, np.nan)
plt.plot(charlotte_rainfall['dt'], charlotte_rainfall["Rainfall"])
plt.ylabel('mm/15min')
plt.gcf().autofmt_xdate()
charlotte_rainfall.head()
Explanation: Replace invalid data with NaNs and plot again
End of explanation
charlotte_24h_rainfall = pd.DataFrame()
charlotte_24h_rainfall['mean_rain'] = charlotte_rainfall.Rainfall.resample('D').mean()
charlotte_24h_rainfall['accum_rain'] = charlotte_rainfall.Rainfall.resample('D').sum()
charlotte_24h_rainfall.head()
plt.plot(charlotte_24h_rainfall["accum_rain"])
plt.ylabel('mm/24h')
plt.gcf().autofmt_xdate()
plt.plot(charlotte_24h_rainfall["mean_rain"])
plt.ylabel(r'mm/15min ($\varnothing$ of 24h)')
plt.gcf().autofmt_xdate()
Explanation: Resample the 10-min dataset to 24h accumulated rainfall data
End of explanation
charlotte_1h_rainfall = pd.DataFrame()
charlotte_1h_rainfall['mean_rain'] = charlotte_rainfall.Rainfall.resample('H').mean()
charlotte_1h_rainfall['accum_rain'] = charlotte_rainfall.Rainfall.resample('H').sum()
charlotte_1h_rainfall.head()
plt.plot(charlotte_1h_rainfall["accum_rain"])
plt.ylabel('mm/h')
plt.gcf().autofmt_xdate()
plt.plot(charlotte_1h_rainfall["mean_rain"])
plt.ylabel(r'mm/15min ($\varnothing$ of 1h)')
plt.gcf().autofmt_xdate()
Explanation: Resample 15 min data to 1h accumulated dataset
End of explanation
charlotte_summer_1h_rainfall = charlotte_1h_rainfall.loc[(charlotte_1h_rainfall.index.month>=4) & (charlotte_1h_rainfall.index.month<=9)]
plt.plot(charlotte_summer_1h_rainfall["accum_rain"])
plt.ylabel('mm/h')
plt.gcf().autofmt_xdate()
Explanation: Select only summer months (April - Sept)
End of explanation
mask_start = (charlotte_1h_rainfall.index.month >= 1) & (charlotte_1h_rainfall.index.month <= 3)
mask_end = (charlotte_1h_rainfall.index.month >= 10) & (charlotte_1h_rainfall.index.month <= 12)
mask = mask_start | mask_end
charlotte_winter_1h_rainfall = charlotte_1h_rainfall.loc[mask]
plt.plot(charlotte_winter_1h_rainfall["accum_rain"])
plt.ylabel('mm/h')
plt.gcf().autofmt_xdate()
charlotte_winter_1h_rainfall.head()
Explanation: Select only winter months (Oct - Mar)
End of explanation
charlotte_monthly_rainfall = pd.DataFrame()
charlotte_monthly_rainfall['mean_rain'] = charlotte_rainfall.Rainfall.resample('M').mean()
charlotte_monthly_rainfall['accum_rain'] = charlotte_rainfall.Rainfall.resample('M').sum()
plt.plot(charlotte_monthly_rainfall["accum_rain"])
plt.ylabel('mm/month')
plt.gcf().autofmt_xdate()
plt.plot(charlotte_monthly_rainfall["mean_rain"])
plt.ylabel(r'mm/15min ($\varnothing$ per month)')
plt.gcf().autofmt_xdate()
Explanation: Resample 15 min dataset to monthly accumulated dataset
End of explanation
print('Mean: %s' % str(charlotte_rainfall.Rainfall.mean()))
print('Std: %s' % str(charlotte_rainfall.Rainfall.std()))
print('Skew: %s' % str(charlotte_rainfall.Rainfall.skew()))
Explanation: Answering the assignments
1. General statistics for 24-hour and 15-min datasets: compute mean, standard deviation, skewness; plot histograms
15 min dataset
Mean, standard deviation and skewness of the 15 min dataset
End of explanation
charlotte_rainfall.Rainfall.hist(bins = 100)
plt.xlabel('mm/15min')
plt.gca().set_yscale("log")
Explanation: Histogram of the data
End of explanation
cur_data = charlotte_rainfall.Rainfall.loc[charlotte_rainfall.Rainfall>0]
hist_d = plt.hist(cur_data, bins=100)
plt.xlabel('mm/15min')
plt.gca().set_yscale("log")
Explanation: Histogram of the data without zeros
End of explanation
print('Mean: %s' % str(charlotte_24h_rainfall.accum_rain.mean()))
print('Std: %s' % str(charlotte_24h_rainfall.accum_rain.std()))
print('Skew: %s' % str(charlotte_24h_rainfall.accum_rain.skew()))
Explanation: 24h accumulated dataset
Mean, standard deviation and skewness of 24h accumulated dataset
End of explanation
charlotte_24h_rainfall.accum_rain.hist(bins = 100)
plt.xlabel('mm/24h')
plt.gca().set_yscale("log")
charlotte_24h_rainfall.mean_rain.hist(bins = 100)
plt.xlabel(r'mm/15min ($\varnothing$ per 24h)')
plt.gca().set_yscale("log")
Explanation: Histogram of the dataset
End of explanation
cur_data = charlotte_24h_rainfall.accum_rain.loc[charlotte_24h_rainfall.accum_rain>0]
hist_d = plt.hist(cur_data, bins=100)
plt.xlabel('mm/24h')
plt.gca().set_yscale("log")
Explanation: Histogram without zeros
End of explanation
charlotte_monthly_rainfall['mon'] = charlotte_monthly_rainfall.index.month
charlotte_monthly_rainfall['year'] = charlotte_monthly_rainfall.index.year
charlotte_monthly_rainfall.boxplot(column=['accum_rain'], by='mon', sym='+')
plt.ylabel('mm/month')
Explanation: 2. a. Analysis of seasonal cycles: create boxplots for monthly totals across all years
Boxplot of monthly totals
End of explanation
charlotte_monthly_rainfall.dropna().boxplot(column=['accum_rain'], by='year', sym='+')
plt.ylabel('mm/month')
plt.gcf().autofmt_xdate()
Explanation: Or on a yearly scale:
End of explanation
charlotte_1h_rainfall['hour'] = charlotte_1h_rainfall.index.hour
charlotte_1h_rainfall.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')
Explanation: 2. b. Analysis of diurnal cycles: create boxplots for hourly totals for entire dataseries
End of explanation
cur_df = charlotte_1h_rainfall.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')
Explanation: Neglecting events < 1mm/h
End of explanation
cur_df = charlotte_1h_rainfall.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')
Explanation: Neglecting events < 3mm/h
End of explanation
pd.options.mode.chained_assignment = None  # default='warn'
charlotte_summer_1h_rainfall['hour'] = charlotte_summer_1h_rainfall.index.hour
charlotte_summer_1h_rainfall.boxplot(column=['accum_rain'], by='hour', sym='+')
Explanation: 2. c. Variation of diurnal cycles with seasons: create boxplots for hourly totals for summer season (April – September) and for winter season (October-March)
Merge summer hourly data
End of explanation
cur_df = charlotte_summer_1h_rainfall.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')
Explanation: Neglecting events <1mm/hour
End of explanation
cur_df = charlotte_summer_1h_rainfall.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')
Explanation: Neglecting events <3mm/hour
End of explanation
charlotte_winter_1h_rainfall['hour'] = charlotte_winter_1h_rainfall.index.hour
charlotte_winter_1h_rainfall.boxplot(column=['accum_rain'], by='hour', sym='+')
plt.ylabel('mm/h')
Explanation: Merge hourly winter data
End of explanation
cur_df = charlotte_winter_1h_rainfall.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')
Explanation: Neglecting events <1mm/h
End of explanation
cur_df = charlotte_winter_1h_rainfall.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')
Explanation: Neglecting events <3mm/h
End of explanation
charlotte_1h_exceeds = charlotte_1h_rainfall.accum_rain[charlotte_1h_rainfall.accum_rain>10]
Explanation: 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
End of explanation
print(len(charlotte_1h_exceeds))
y = np.array(charlotte_1h_exceeds)
N = len(y)
x = range(N)
width = 1
plt.bar(x, y, width)
plt.ylabel('mm/h')
Explanation: Amount of hourly events
End of explanation
charlotte_1h_exceeds_summer = charlotte_summer_1h_rainfall.accum_rain[charlotte_summer_1h_rainfall.accum_rain>10]
y = np.array(charlotte_1h_exceeds_summer)
N = len(y)
x = range(N)
width = 1
plt.bar(x, y, width)
plt.ylabel('mm/h')
Explanation: 10 mm/h events in summer periods
End of explanation
print(len(charlotte_1h_exceeds_summer))
Explanation: Amount of hourly events
End of explanation
plt.plot(charlotte_1h_exceeds)
plt.gcf().autofmt_xdate()
charlotte_1h_exceeds.hist(bins=100)
from scipy.stats import genextreme
x = np.linspace(0, 80, 1000)
y = np.array(charlotte_1h_exceeds[:])
np.seterr(divide='ignore', invalid='ignore')
genextreme.fit(y)
pdf = plt.plot(x, genextreme.pdf(x, *genextreme.fit(y)))
pdf_hist = plt.hist(y, bins=50, normed=True, histtype='stepfilled', alpha=0.8)
Explanation: 3. Fit GEV-distribution for POT values in the time series
3. a. Create plots: histogram and GEV fit and interpret
End of explanation
genextreme.ppf((1-1/1), *genextreme.fit(y))
genextreme.ppf((1-1/10), *genextreme.fit(y))
genextreme.ppf((1-1/100), *genextreme.fit(y))
Explanation: 3. c. Compute rainfall amounts associated with return periods of 1 year, 10 years and 100 years
End of explanation
from scipy.stats import genpareto
temp_monthly = charlotte_1h_rainfall.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, 100, 1000)
pdf = plt.plot(x, genextreme.pdf(x, *genextreme.fit(block_max_y)))
pdf_hist = plt.hist(block_max_y, bins=50, normed=True, histtype='stepfilled', alpha=0.8)
Explanation: Update 10.10.2017
Block maxima & GEV
End of explanation
genextreme.fit(block_max_y)
genextreme.ppf((1-1/10), *genextreme.fit(block_max_y))
Explanation: GEV and block maxima of monthly maxima of 1h data
End of explanation
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)
Explanation: POT & GPD
End of explanation
genpareto.fit(y)
genpareto.ppf((1-1/10), *genpareto.fit(y))
Explanation: GPD and POT of data>10mm/h
End of explanation
event_occurences = pd.DataFrame(charlotte_1h_exceeds)
event_occurences['hour'] = event_occurences.index.hour
event_occurences.boxplot(column=['accum_rain'], by='hour', sym='+')
Explanation: Boxplot of POT values
End of explanation
event_occurences.hour.value_counts(sort=False)
# plt.plot(asd.hour.value_counts(sort=False))
cur_hist = plt.hist(event_occurences.hour, bins=24, histtype='stepfilled')
plt.xticks(range(24))
plt.xlabel('hour')
Explanation: Number of occurences per hour
End of explanation
<END_TASK>