Skip to content

Commit

Permalink
add option to give already resamples ams to the ams_analysis function
Browse files Browse the repository at this point in the history
  • Loading branch information
MAfarrag committed Jun 21, 2023
1 parent 5fe16c0 commit 8eac5d4
Showing 1 changed file with 27 additions and 20 deletions.
47 changes: 27 additions & 20 deletions statista/eva.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,10 +12,11 @@

def ams_analysis(
time_series_df: DataFrame,
ams: bool = False,
ams_start: str = "A-OCT",
save_plots: bool = False,
save_to: str = None,
filter_out: Union[bool, float, int] = False,
filter_out: Union[float, int] = None,
distribution: str = "GEV",
method: str = "lmoments",
estimate_parameters: bool = False,
Expand All @@ -31,6 +32,8 @@ def ams_analysis(
----------
time_series_df : [DataFrame]
DataFrame containing multiple time series to do the statistical analysis on.
ams: [bool]
True if the the given time series is annual mean series. Default is False.
ams_start: [str]
The beginning of the year which is used to resample the time series to get the annual maximum series.
Default is"A-OCT".
Expand Down Expand Up @@ -129,30 +132,32 @@ def ams_analysis(

for i in gauges:
QTS = time_series_df.loc[:, i]
# The time series is resampled to the annual maxima, and turned into a
# numpy array.
# The time series is resampled to the annual maxima, and turned into a numpy array.
# The hydrological year is 1-Nov/31-Oct (from Petrow and Merz, 2009, JoH).
ams = QTS.resample(ams_start).max().values
if not ams:
ams_df = QTS.resample(ams_start).max().values
else:
ams_df = QTS.values

if not isinstance(filter_out, bool):
ams = ams[ams != filter_out]
if filter_out is not None:
ams_df = ams_df[ams_df != filter_out]

if estimate_parameters:
# TODO: still to be tested and prepared for GEV
# estimate the parameters through an optimization
# alpha = (np.sqrt(6) / np.pi) * ams.std()
# beta = ams.mean() - 0.5772 * alpha
# alpha = (np.sqrt(6) / np.pi) * ams_df.std()
# beta = ams_df.mean() - 0.5772 * alpha
# param_dist = [beta, alpha]
threshold = np.quantile(ams, quartile)
threshold = np.quantile(ams_df, quartile)
if distribution == "GEV":
dist = GEV(ams)
dist = GEV(ams_df)
param_dist = dist.estimateParameter(
method="optimization",
ObjFunc=Gumbel.ObjectiveFn,
threshold=threshold,
)
else:
dist = Gumbel(ams)
dist = Gumbel(ams_df)
param_dist = dist.estimateParameter(
method="optimization",
ObjFunc=Gumbel.ObjectiveFn,
Expand All @@ -162,12 +167,12 @@ def ams_analysis(
# estimate the parameters through maximum liklehood method
try:
if distribution == "GEV":
dist = GEV(ams)
dist = GEV(ams_df)
# defult parameter estimation method is maximum liklihood method
param_dist = dist.estimateParameter(method=method)
else:
# A gumbel distribution is fitted to the annual maxima
dist = Gumbel(ams)
dist = Gumbel(ams_df)
# defult parameter estimation method is maximum liklihood method
param_dist = dist.estimateParameter(method=method)
except Exception as e:
Expand Down Expand Up @@ -198,10 +203,10 @@ def ams_analysis(
Qrp = dist.theporeticalEstimate(param_dist[0], param_dist[1], F)

# to get the Non Exceedance probability for a specific Value
# sort the ams
ams.sort()
# sort the ams_df
ams_df.sort()
# calculate the F (Exceedence probability based on weibul)
cdf_Weibul = PlottingPosition.weibul(ams)
cdf_Weibul = PlottingPosition.weibul(ams_df)
# Gumbel.probapilityPlot method calculates the theoretical values
# based on the Gumbel distribution
# parameters, theoretical cdf (or weibul), and calculate the confidence interval
Expand Down Expand Up @@ -240,10 +245,12 @@ def ams_analysis(
statistical_properties.loc[i, "max"] = QTS.max()
statistical_properties.loc[i, "t_beg"] = QTS.index.min()
statistical_properties.loc[i, "t_end"] = QTS.index.max()
statistical_properties.loc[i, "nyr"] = (
statistical_properties.loc[i, "t_end"]
- statistical_properties.loc[i, "t_beg"]
).days / 365.25
if not ams:
statistical_properties.loc[i, "nyr"] = (
statistical_properties.loc[i, "t_end"]
- statistical_properties.loc[i, "t_beg"]
).days / 365.25

for irp, irp_name in zip(Qrp, rp_name):
statistical_properties.loc[i, irp_name] = irp

Expand Down

0 comments on commit 8eac5d4

Please sign in to comment.