Repository: MET-OM/metocean-stats (language: python)

Issue:
merge main into dev



Discussion / hints:

Fix patch:
diff --git a/docs/2D_contours.png b/docs/2D_contours.png
new file mode 100644
index 0000000..a856f71
Binary files /dev/null and b/docs/2D_contours.png differ
diff --git a/docs/index.rst b/docs/index.rst
index b84d6ed..ce724c7 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -13,7 +13,20 @@ The package contains functions that:
 
 Installing **metocean-stats**
 =============================================
-1. Install anaconda3 or miniconda3
+Alternative 1: Using Mambaforge (alternative to Miniconda)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+1. Install `mambaforge <https://mamba.readthedocs.io/en/latest/installation.html>`_ (`download <https://github.com/conda-forge/miniforge#mambaforge>`_)
+2. Set up a *Python 3* environment for metocean-stats and install metocean-stats
+
+.. code-block:: bash
+
+   $ mamba create -n metocean-stats python=3 metocean-stats
+   $ conda activate metocean-stats
+
+Alternative 2: Using Mambaforge (alternative to Miniconda) and Git
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+1. Install `mambaforge <https://mamba.readthedocs.io/en/latest/installation.html>`_ (`download <https://github.com/conda-forge/miniforge#mambaforge>`_)
 2. Clone metocean-stats:
 
 .. code-block:: bash
@@ -25,16 +38,17 @@ Installing **metocean-stats**
 
 .. code-block:: bash
 
-  $ conda config --add channels conda-forge
-  $ conda env create -f environment.yml
+  $ mamba env create -f environment.yml
   $ conda activate metocean-stats
-  $ pip install -e .
-  
-To update the enviroment using a new environment.yml, run:
+  $ pip install --no-deps -e .
+
+This installs the metocean-stats as an editable package. Therefore, you can directly make changes to the repository or fetch the newest changes with :code:`git pull`. 
+
+To update the local conda environment in case of new dependencies added to environment.yml:
 
 .. code-block:: bash
 
-   $ conda env update --file environment.yml --prune
+  $ mamba env update -f environment.yml
 
 Download metocean data
 =====================================
@@ -104,7 +118,7 @@ Create table with monthly percentiles:
    des , 0.4 , 1.3 , 1.4 , 2.8 , 3.5 
    Annual , 0.4 , 1.3 , 1.4 , 2.8 , 3.5
 
-Create table with monthly min, mean, and max values:
+Create table with monthly min, mean, and max based on monthly values:
 
 .. code-block:: python
 
@@ -224,8 +238,20 @@ Create time series of Hs with return levels using POT and Annual Maximum(GEV):
            output_file='return_levels_GEV.png')
 
 .. image:: return_levels_GEV.png
+  :width: 700   
+
+Plot joint Hs-Tp contours for different return periods (to get the contour values use: extreme_stats.get_joint_2D_contour):
+
+.. code-block:: python
+   
+   extreme_stats.plot_joint_2D_contour(data=ds.data,var1='hs',var2='tp', periods=[50,100], 
+                                       output_file='2D_contours.png')
+
+.. image:: 2D_contours.png
   :width: 700
 
+
+
 Profile Statistics
 =====================================
 To generate profile stastistics, import profile_stats module: 
@@ -249,6 +275,9 @@ Estimate and plot mean wind profile:
 .. image:: wind_profile.png
   :height: 500
 
+
+
+
 .. toctree::
    :maxdepth: 2
    :caption: Contents:
diff --git a/environment.yml b/environment.yml
index 681c495..17d0e73 100644
--- a/environment.yml
+++ b/environment.yml
@@ -12,3 +12,5 @@ dependencies:
   - windrose
   - seaborn>=0.12.2
   - metocean-api
+  - pyextremes
+  - virocon
diff --git a/metocean_stats/stats/extreme_stats.py b/metocean_stats/stats/extreme_stats.py
index f075de9..22c7253 100644
--- a/metocean_stats/stats/extreme_stats.py
+++ b/metocean_stats/stats/extreme_stats.py
@@ -1,24 +1,24 @@
+from typing import Dict, Sequence
 import pandas as pd
 import numpy as np
 import matplotlib.pyplot as plt
+import xarray as xr
 
-def return_levels_pot(data,var,threshold=None,periods=[50,100], output_file='return_levels_pot.pdf'):
+def return_levels_pot(data,var,threshold=None,periods=[50,100], output_file=False):
     """
     Subroutine written by clio-met (July 2023)
     Purpose: calculates the return levels (rl) from a time series (ts) and given threshold for different return periods (periods)
     - Inputs:
     1) ts is a dataframe with time enabled containing the hourly time series of the variable of interest (ex: wind speed)
     2) threshold is a scalar 
-    2) periods is a numpy array containing the return periods (in years) of interest (ex: np.array([20,30,50,100]))
+    3) periods is a numpy array containing the return periods (in years) of interest (ex: np.array([20,30,50,100]))
     - Follows the POT method 
     """
     if threshold == None:
         threshold = get_threshold_os(data=data, var=var)
     else:
         pass
-    # Required modules
     periods = np.array(periods, dtype=float)
-    # Take yearly maximum
     ts = data[var]
     # select the values above the threshold in the timeseries ts
     it=np.where(ts.values>=threshold)[0]
@@ -43,19 +43,19 @@ def return_levels_pot(data,var,threshold=None,periods=[50,100], output_file='ret
     lambd=np.mean(ns_yr.values,axis=0)
     # Fit the Generalized Pareto distribution to the selected values
     from scipy.stats import genpareto
-    shape,loc,scale = genpareto.fit(sel_val.values,loc=threshold)
-    rl = genpareto.isf(1/(lambd*periods), shape, loc, scale)
+    shape,loc,scale = genpareto.fit(sel_val.values-threshold)
+    rl = genpareto.isf(1/(lambd*periods), shape, loc, scale) + threshold
     # If the return level is unrealistic (usually very high), set rl to a special number
     rl=np.where(rl>=3*threshold,np.nan,rl)
     if output_file == False:
         pass
     else:
         plot_return_levels(data,var,rl,periods,output_file,it_selected_max=ts.iloc[it_selected_max].index)
-
+        #probplot(data=sel_val.values, sparams=(shape, loc, scale))
     return rl
 
 
-def return_levels_annual_max(data,var='hs',periods=[50,100,1000],method='GEV',output_file='return_levels_annual_max.pdf'): 
+def return_levels_annual_max(data,var='hs',periods=[50,100,1000],method='GEV',output_file=False): 
     """
     The function is written by dung-manh-nguyen and KonstantinChri.
     It calulates return value estimates for different periods using different methods 
@@ -93,8 +93,9 @@ def return_levels_annual_max(data,var='hs',periods=[50,100,1000],method='GEV',ou
         pass
     else:
         plot_return_levels(data,var,rl,periods,output_file,it_selected_max)
+        #probplot(data=data[var].loc[it_selected_max], sparams=(shape, loc, scale))
 
-
+    
     return rl
 
 def get_threshold_os(data,var):
@@ -109,6 +110,64 @@ def get_threshold_os(data,var):
     min_ym=yearly_max.min()
     return min_ym
 
+
+def return_levels_weibull_2p(data,var,periods=[50,100,1000],threshold=None,r="48H"):
+    import scipy.stats as stats
+    from pyextremes import get_extremes
+    # how to use, to test this 
+    #return_period = np.array([1,10,100,10000],dtype=float)
+    #return_values = RVE_Weibull_2P(df.hs,return_period,threshold=6.2)
+    #print (return_values)
+    if threshold == None:
+        threshold = get_threshold_os(data=data, var=var)
+    else:
+        pass
+    extremes = get_extremes(data[var], method="POT", threshold=threshold, r=r)
+    
+    # Fit a 2-parameter Weibull distribution to the data
+    shape, loc, scale = stats.weibull_min.fit(extremes, floc=0)
+       
+    years = data.index.year
+    yr_num = years[-1]-years[0]+1
+    length_data = extremes.shape[0]
+    time_step = yr_num*365.2422*24/length_data # in hours 
+    return_period = np.array(periods)*24*365.2422/time_step # years is converted to K-th
+    
+    # Calculate the 2-parameter Weibull return value
+    return_value = stats.weibull_min.isf(1/return_period, shape, loc, scale)
+    
+    for i in range(len(return_value)) : 
+        return_value[i] = round(return_value[i],1)
+    
+    return return_value
+
+def return_levels_exp(data,var='hs',periods=[50,100,1000],threshold=None, r="48H"):
+    from scipy.stats import expon
+    from pyextremes import get_extremes
+    # how to use this function 
+    #return_periods = np.array([1.6, 10, 100, 10000]) # in years
+    #return_values = RVE_EXP(df.hs,return_periods,4)
+    #print (return_values)
+    if threshold == None:
+        threshold = get_threshold_os(data=data, var=var)
+    else:
+        pass
+    extremes = get_extremes(data[var], method="POT", threshold=threshold, r=r)
+    loc, scale = expon.fit(extremes)
+    #print (loc,scale)
+    
+    years = data.index.year
+    yr_num = years[-1]-years[0]+1
+    length_data = extremes.shape[0]
+    interval = yr_num*365.2422*24/length_data # in hours 
+    
+    return_periods = np.array(periods)*24*365.2422/interval # years is converted to K-th
+        
+    return_levels = expon.isf(1/return_periods, loc, scale)
+    
+    return return_levels
+
+
 def plot_return_levels(data,var,rl,periods, output_file,it_selected_max=[]):
     color = plt.cm.rainbow(np.linspace(0, 1, len(rl)))
     fig, ax = plt.subplots(figsize=(12,6))
@@ -128,4 +187,146 @@ def plot_return_levels(data,var,rl,periods, output_file,it_selected_max=[]):
     plt.savefig(output_file)
     plt.close()
     
-    return
\ No newline at end of file
+    return
+
+
+def probplot(data, sparams):    
+    import scipy.stats as stats
+    stats.probplot(data,sparams=sparams, dist=stats.genpareto,fit=True, plot=plt)
+    plt.grid()
+    plt.axline((data[0], data[0]), slope=1,label='y=x')
+    plt.legend()
+    plt.tight_layout()
+    plt.show()
+    return
+
+def return_value_plot(obs_return_periods, obs_return_levels,model_return_periods, model_return_levels, model_method=''):    
+    fig, ax = plt.subplots()   
+    ax.scatter(obs_return_periods, obs_return_levels, marker="o",s=20,lw=1,facecolor="k",edgecolor="w",zorder=20)
+    ax.plot(model_return_periods, model_return_levels,color='red',label=model_method)
+    ax.semilogx()
+    ax.grid(True, which="both")
+    ax.set_xlabel("Return period")
+    ax.set_ylabel("Return levels")
+    plt.legend()
+    plt.tight_layout()
+    plt.show()
+    #alpha = 0.95
+    #cil, ciu = np.quantile(a=rl_sample, q=[(1 - alpha) / 2, (1 + alpha) / 2])
+    return
+
+
+
+
+def get_empirical_return_levels(data,var,method="POT",block_size="365.2425D"):
+    """
+    Return an estimation of empirical/observed [return levels, return periods] 
+    """
+    from pyextremes import get_extremes, get_return_periods
+    if method == 'BM':
+        extremes = get_extremes(ts=data[var],method=method,block_size=block_size)
+        rp = get_return_periods(ts=data[var],extremes=extremes,extremes_method=method,extremes_type="high",return_period_size=block_size)
+    elif method == 'POT':
+        extremes = get_extremes(ts=data[var],method=method,threshold=get_threshold_os(data=data, var=var))
+        rp = get_return_periods(ts=data[var],extremes=extremes,extremes_method=method,extremes_type="high")
+    return rp['return period'], rp[var]
+
+
+def diagnostic_return_level_plot(data,var,model_method,periods=np.arange(0.1,1000,0.1),threshold=None):
+    if model_method == 'POT':
+        return_value_plot(obs_return_periods=get_empirical_return_levels(data,var)[0], obs_return_levels=get_empirical_return_levels(data,var)[1],
+                          model_return_periods=periods, model_return_levels=return_levels_pot(data,var,threshold=threshold,periods=periods), model_method=model_method)
+    elif model_method == 'GEV':
+        return_value_plot(obs_return_periods=get_empirical_return_levels(data,var)[0], obs_return_levels=get_empirical_return_levels(data,var)[1],
+                          model_return_periods=periods, model_return_levels=return_levels_annual_max(data,var,method='GEV',periods=periods), model_method=model_method)
+    elif model_method == 'GUM':
+        return_value_plot(obs_return_periods=get_empirical_return_levels(data,var)[0], obs_return_levels=get_empirical_return_levels(data,var)[1],
+                          model_return_periods=periods, model_return_levels=return_levels_annual_max(data,var,method='GUM',periods=periods), model_method=model_method)
+    elif model_method == 'Weibull_2P':
+        return_value_plot(obs_return_periods=get_empirical_return_levels(data,var)[0], obs_return_levels=get_empirical_return_levels(data,var)[1],
+                          model_return_periods=periods, model_return_levels=return_levels_weibull_2p(data,var,threshold=threshold,periods=periods), model_method=model_method)
+    elif model_method == 'EXP':
+        return_value_plot(obs_return_periods=get_empirical_return_levels(data,var)[0], obs_return_levels=get_empirical_return_levels(data,var)[1],
+                          model_return_periods=periods, model_return_levels=return_levels_exp(data,var,threshold=threshold,periods=periods), model_method=model_method)
+        
+
+
+def get_joint_2D_contour(data=pd.DataFrame,var1='hs', var2='tp', periods=[50,100]) -> Sequence[Dict]:
+    """Compute a joint contour for the given return periods. 
+    Input:
+        data: pd.DataFrame
+        var1: e.g., 'hs'
+        var2: e.g., 'tp'
+        return_periods: A list of return periods in years, default [50,100]
+    Output:
+         contours : list of joint contours, 
+         i.e.,  number of contours based on given return_periods 
+    """
+    from virocon import get_OMAE2020_Hs_Tz,get_DNVGL_Hs_U,get_OMAE2020_Hs_Tz,get_OMAE2020_V_Hs, GlobalHierarchicalModel,IFORMContour, ISORMContour     
+    # Define 2D joint distribution model
+    dist_descriptions, fit_descriptions, _ = get_OMAE2020_Hs_Tz()
+    model = GlobalHierarchicalModel(dist_descriptions)
+    dt = (data.index[1]-data.index[0]).total_seconds()/3600  # duration in hours
+    data_2D =  np.transpose(np.array([data[var1],data[var2]]))
+    model.fit(data_2D, fit_descriptions=fit_descriptions)
+    contours = []
+    for rp in periods:
+        alpha = 1 / (rp * 365.25 * 24 / dt)
+        contour = IFORMContour(model, alpha)
+        coords = contour.coordinates
+        x = coords[:, 1].tolist()
+        y = coords[:, 0].tolist()
+        contours.append({
+            'return_period': rp,
+            'x': x,
+            'y': y
+        })
+    return contours, data_2D
+
+def plot_joint_2D_contour(data=pd.DataFrame,var1='hs', var2='tp', periods=[50,100], output_file='2D_contours.png') -> Sequence[Dict]:
+    contours, data_2D = get_joint_2D_contour(data=data,var1=var1,var2=var2, periods=periods)
+    """Plot joint contours for the given return periods. 
+    Input:
+        data: pd.DataFrame
+        var1: e.g., 'hs'
+        var2: e.g., 'tp'
+        return_periods: A list of return periods in years, default [50,100]
+        output_file: Path to save the plot 
+    Output:
+         figure with 2Djoint contours
+    """
+    # Plot the contours and save the plot
+
+    fig, ax = plt.subplots()
+    if len(data_2D)>0:
+        import seaborn as sns
+        sns.set_theme(style="ticks")
+        sns.scatterplot(x=data_2D[:,1], y=data_2D[:,0], ax=ax, color='black')
+    else:
+        pass
+
+    color = plt.cm.rainbow(np.linspace(0, 1, len(periods)))
+    # Compute an IFORM contour with a return period of N years
+    #loop over contours, get index
+
+    for i,contour in enumerate(contours):
+        # Plot the contour
+        x = []
+        x.extend(contour['x'])
+        x.append(x[0])
+
+        y = []
+        y.extend(contour['y'])
+        y.append(y[0])
+
+        rp = contour['return_period']
+        ax.plot(x, y, color=color[i],label=str(rp)+'y')
+
+        ax.set_xlabel(var2,fontsize='16')
+        ax.set_ylabel(var1,fontsize='16')
+
+    plt.grid()
+    plt.title(output_file.split('.')[0],fontsize=18)
+    plt.legend()
+    plt.savefig(output_file,dpi=300)
+    return fig, ax
diff --git a/metocean_stats/stats/general_stats.py b/metocean_stats/stats/general_stats.py
index 0a2730c..8faf8bd 100644
--- a/metocean_stats/stats/general_stats.py
+++ b/metocean_stats/stats/general_stats.py
@@ -227,7 +227,7 @@ def table_monthly_percentile(data,var,output_file='var_monthly_percentile.txt'):
 def table_monthly_min_mean_max(data, var,output_file='montly_min_mean_max.txt') :  
     """
     The function is written by dung-manh-nguyen and KonstantinChri.
-    It calculates min, mean, max in each month  and the omni values 
+    It calculates monthly min, mean, max based on monthly maxima. 
     data : panda series
     var  : variable 
     output_file: extension .txt for latex table or .csv for csv table
@@ -256,7 +256,7 @@ def table_monthly_min_mean_max(data, var,output_file='montly_min_mean_max.txt')
         min_year = annual_max.min()
         mean_year = annual_max.mean()
         max_year = annual_max.max()
-        f.write('Annual & ' + str(min_year) + ' & ' + str(round(mean_year,1)) + ' & ' + str(max_year) + ' \\\\' + '\n')
+        f.write('Annual Max. & ' + str(min_year) + ' & ' + str(round(mean_year,1)) + ' & ' + str(max_year) + ' \\\\' + '\n')
             
         f.write('\hline' + '\n')
         f.write('\end{tabular}' + '\n')
diff --git a/version.py b/version.py
index e74fdb9..2147b70 100644
--- a/version.py
+++ b/version.py
@@ -1,4 +1,4 @@
-__version__ = "0.1.4"
+__version__ = "0.1.5"
 
 
 def git_describe():


Test patch:
diff --git a/tests/test_stats.py b/tests/test_stats.py
index 6e1c4ba..7962651 100644
--- a/tests/test_stats.py
+++ b/tests/test_stats.py
@@ -38,4 +38,9 @@ def test_return_levels_annual_max(ds=ds):
     extreme_stats.return_levels_annual_max(data=ds.data, var='hs', periods=[20,100], output_file=False) 
 
 def test_return_levels_annual_max(ds=ds):
-    profile_stats.mean_profile(data = ds.data, vars = ['wind_speed_10m','wind_speed_20m','wind_speed_50m','wind_speed_100m','wind_speed_250m','wind_speed_500m','wind_speed_750m'],height_levels=[10,20,50,100,250,500,750], perc = [25,75], output_file=False)
\ No newline at end of file
+    profile_stats.mean_profile(data = ds.data, vars = ['wind_speed_10m','wind_speed_20m','wind_speed_50m','wind_speed_100m','wind_speed_250m','wind_speed_500m','wind_speed_750m'],height_levels=[10,20,50,100,250,500,750], perc = [25,75], output_file=False)
+
+def test_joint_2D_contour(ds=ds):
+    extreme_stats.plot_joint_2D_contour(data=ds.data,var1='hs',var2='tp', periods=[20,100], output_file='test.png')
+    os.remove('test.png')
+