# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# +
import AnalysisFunctions as af

import pandas as pd  
#defaultdict to use nested dictionaries
from collections import defaultdict

import matplotlib.pyplot as plt
import statsmodels.api as sm
import numpy as np

import warnings
warnings.filterwarnings('ignore')

import dill

from IPython.display import IFrame
# -

# # 1. Pre-processing:

# - Choose the initialization time of the simulation, 
# - Organize the runoff and precipitation data in dataframes altogether with their respective observation intervals
# - Calculate the quantiles on the distribution of ensemble forecasts 
# - Calculate the runoff meteorological medians out of the set of realizations

# Dictionary initialization:
df = af.dictionary()

# Decide simulation initialization time:
sim_start = '2018-10-26 12:00:00'

# Some parameters for the basin:
Verzasca_area = 186*1000.0**2 #m2
conv_factor = Verzasca_area/(1000.0*3600.0) 

# +
# Runoff observation: open observation dataframe
obs_pattern = '/home/ciccuz/hydro/prevah/runoff/2605.dat'
obs_columns = ['year','month','day','hour','runoff']
obs_df = pd.DataFrame(pd.read_csv(obs_pattern, names=obs_columns, delim_whitespace=True, header=None))
obs_df['date'] = pd.to_datetime(obs_df[['year', 'month', 'day', 'hour']])


# Precipitation observation: open precipitation observation dataframe obtained by cosmo1+pluviometer 
#data concatenated series before the initialization of the model
prec_obs_df = af.dictionary(pattern="/home/ciccuz/hydro/PrecObs/cosmo1_{simul_time}/{otherstuff}",
                        folders_pattern = '/home/ciccuz/hydro/PrecObs/cosmo1_*')
prec_obs_series= dill.load( open( "prec_obs/prec_obs_series.txt", "rb" ) )
# -

# Extract from the dictionary the dataframe containing all the different realizations of the same event: 
#every ensemble member and parameter set combination for the runoff, every ensemble member for the precipitation.
ens_df_prec = af.ensemble_df(df, sim_start, Verzasca_area, 'P-kor')
ens_df_runoff = af.ensemble_df(df, sim_start, Verzasca_area,'RGES')

# Calculate the quantiles for the variable chosen considering all the different realizations for the 120h ahead.
quant_prec = af.quantiles(ens_df_prec)
quant_runoff = af.quantiles(ens_df_runoff)

# +
# Define the subset of runoff and precipitation observation based on quantiles dataframe date boundaries
obs_indexes_runoff = obs_df.loc[obs_df.index[obs_df['date'] == str(quant_runoff.date[0])] |
        obs_df.index[obs_df['date'] == str(quant_runoff.date[119])]]
obs_indexes_prec = prec_obs_series.loc[prec_obs_series.index[prec_obs_series['date'] == str(quant_runoff.date[0])] |
        prec_obs_series.index[prec_obs_series['date'] == str(quant_runoff.date[119])]]

obs_subset = obs_df.loc[obs_indexes_runoff.index[0]:obs_indexes_runoff.index[1]]
prec_obs_subset = prec_obs_series.loc[obs_indexes_prec.index[0]:obs_indexes_prec.index[1]]

# +
# Meteorological medians:
# Select groups of realizations based on the same ensemble members:
# dictionaries sorted by ensemble members
rm_groups_runoff = af.ens_param_groups(ens_df_runoff)[0]

# Quantiles dictionaries from above rm groups dictionary
quant_rm_dict = lambda: defaultdict(quant_rm_dict)
quant_rm_groups_runoff = quant_rm_dict()

for rm in range(21):
    quant_rm_groups_runoff[rm] = af.quantiles(rm_groups_runoff[rm])

# Construct a dataframe having all the medians obtained for every group of realizations 
# associated to an ens member
rm_medians = pd.DataFrame(index=range(120))

for rm in range(21):
    rm_medians[rm] = quant_rm_groups_runoff[rm]['0.5']
rm_medians['date'] = quant_rm_groups_runoff[rm]['date']

rm_medians.columns = ['rm00','rm01','rm02','rm03','rm04','rm05','rm06','rm07','rm08','rm09','rm10','rm11','rm12',
                     'rm13','rm14','rm15','rm16','rm17','rm18','rm19','rm20','date']

# Quantiles on rm medians:
quant_rm_medians = af.quantiles(rm_medians)
# -

# ## 1.1. Spaghetti plot for the entire set of realizations:

# One forecast initialization comprises 120 hourly values of forecast, the number of runoff forecasts for every initialization is given by the product of the ensemble members of the meteorological model (21 ens) and the set of hydrological parameters that have been made randomly change (25 pin), so 525. The precipitation forecasts are given by the 21 ensemble members composing the meteorological model.

#Spaghetti plot: 
af.spaghetti_plot(ens_df_runoff, ens_df_prec, obs_subset, prec_obs_subset, sim_start)

IFrame('latex/thesis/figs/uncertainty/spaghetti_all.png', width=900, height=600)

# # 2. Separate different sources of uncertainty:

# ## 2.1. Meteorological uncertainty:

# To have a look at how the meteorological uncertainty impact on the forecast take the meteorological medians (i.e. the 21 medians around the 25 different sets of hydrological parameters) and look at the resulting spread of the simulation, compared to the total spread obtained by the 525 ensemble members.

# Quantify the meteorological uncertainty by plotting the range of spread among all the 21 rm medians obtained:
#af.hydrograph(quant_rm_medians, quant_prec, obs_subset, prec_obs_subset, sim_start, medians=True)
af.comparison_meteo_hydrograph(quant_rm_medians, quant_runoff, quant_prec, obs_subset, prec_obs_subset, sim_start)[1]

IFrame('latex/thesis/figs/uncertainty/hydrograph_meteo.png', width=900, height=600)

# ## 2.2. Hydrological parameters uncertainty:

# #ADD DESCRIPTION

IFrame('latex/thesis/figs/uncertainty/hydrograph_hydro.png', width=900, height=600)

# Look for different ensemble realizations how the hydrological spread behaves: detect three realizations having different behaviours and plot the corresponding spread around their medians.

# #ADD DESCRIPTION

# +
rm_high = 4
rm_medium = 9
rm_low = 7

af.hydrograph_rms(rm_high, rm_medium, rm_low, ens_df_prec, quant_rm_groups_runoff, quant_runoff, 
                  obs_subset, prec_obs_subset, sim_start)
# -

IFrame('latex/thesis/figs/uncertainty/hydro_unc_spread.png', width=900, height=350)

# Quantify the hydrological uncertainty considering the spread between the quantiles around every meteorological median: for every meteo median, in every hourly point, report the total spread range and the IQR and normalize it with the median value of discharge in that point.

# # 3. Peak-box approach:
# Application of the algorithm developed to construct multiple peak-boxes related to different peak events, as an extension of the algorithm developed in Zappa et al. (2013). For more details take a look at the [dedicated  repository](https://github.com/agiord/peakbox)

import individual_scripts.peakbox_v5 as pbk

pbk.peak_box_multipeaks(rm_medians1, obs_subset1, sim_start)

IFrame('latex/thesis/figs/Peakbox/new/pb1.png', width=750, height=600)

# ## 3.1 Peak-box evaluation

# ### - Forecast sharpness

IFrame('latex/thesis/figs/Peakbox/sharp_verif/sharpness_v4_1.png', width=950, height=500)

# ### - Peak median verification

IFrame('latex/thesis/figs/Peakbox/sharp_verif/verification_v4_2.png', width=950, height=500)

# ### - Events detection

IFrame('latex/thesis/figs/Peakbox/sharp_verif/hitmiss_v2.png', width=950, height=380)

# # 4. Cluster analysis

import individual_scripts.cluster_funct as cl

# Perform a hierarchical agglomerative complete-linkage cluster analysis on the meteorological precipitation ensemble members and extract a restricted number of representative members. Goal: check whether the representative members extracted from the few clusters (3-5-7) we obtain are able to give similar/worst/better forecasts than by using the entire set of ensemble forecasts (21), and see which part of the spread are able to cover.

#Plot the dendrogram: visualization of the clustering algorithm applied
cl.clustered_dendrogram(ens_df_prec.drop('date', axis=1), sim_start)
plt.show()

IFrame('latex/thesis/figs/cluster/runoff/dendrogram.png', width=950, height=550)

# +
#Choose the number of clusters (3 or 5):
Nclusters = 3

#extract the representative members
RM = cl.clustered_RM(ens_df_prec.drop('date', axis=1), sim_start, Nclusters)

#extract the sub-dataframe for prec and runoff forecasts containing only the members related to the new extracted representative members:
clust_ens_df_prec = pd.DataFrame()
clust_ens_df_runoff = pd.DataFrame()

for rm_index in range(Nclusters):
    clust_ens_df_prec = pd.concat([clust_ens_df_prec, ens_df_prec.loc[:, ens_df_prec.columns == f'rm{RM[rm_index]:02d}_pin01']], axis=1, sort=False)
    for pin in range(1,26):
        clust_ens_df_runoff = pd.concat([clust_ens_df_runoff, ens_df_runoff.loc[:, ens_df_runoff.columns == f'rm{RM[rm_index]:02d}_pin{pin:02d}']], axis=1, sort=False)

clust_ens_df_prec = pd.concat([clust_ens_df_prec, ens_df_prec.date], axis=1)
clust_ens_df_runoff = pd.concat([clust_ens_df_runoff, ens_df_runoff.date], axis=1)

# Cluster quantiles:
clust_quant_prec = af.quantiles(clust_ens_df_prec)
clust_quant_runoff = af.quantiles(clust_ens_df_runoff)
# -

# Hydrograph plot of clustered forecasts
cl.cluster_hydrograph(clust_quant_runoff, clust_quant_prec, quant_runoff, quant_prec, obs_subset, prec_obs_subset, sim_start, Nclusters=Nclusters)[2]

IFrame('latex/thesis/figs/cluster/runoff/cluster_hydrograph_3_NEW.png', width=950, height=550)

IFrame('latex/thesis/figs/cluster/runoff/cluster_hydrograph_5_NEW.png', width=950, height=550)

IFrame('latex/thesis/figs/cluster/runoff/cluster_hydrograph_7_NEW.png', width=950, height=550)

# ## 4.1 Clustering performance

IFrame('latex/thesis/figs/cluster/cluster_coverage_wIntervals.png', width=700, height=450)

IFrame('latex/thesis/figs/cluster/runoff/ROCa_runoff_TOTvsCLUSTERS.png', width=850, height=650)

StarcoderdataJupyterScriptsDedupFiltered