Merge pull request #91 from MET-OM/dev_cca

Dev cca

Author: KonstantinChri

metocean_stats/plots/extreme.py

@@ -1043,7 +1043,7 @@ def plot_storm_surge_for_given_hs(data: pd.DataFrame, var_surge: str, var_hs: st
     return fig
 
 
-def plot_cca_profiles(data,var='current_speed_',month=None,percentile=None,return_period=None,distribution='GUM',method='default',threshold='default',unit_var='m/s',unit_lev='m',output_file='plot_cca_profiles.png'):
+def plot_cca_profiles(data,var='current_speed_',month=None,percentile=None,return_period=None,distribution='GUM',method='default',threshold='default',max='worst_case',unit_var='m/s',unit_lev='m',output_file='plot_cca_profiles.png'):
     """
     This function plots the CCA profiles for a specific percentile or return period
     data: dataframe
@@ -1052,18 +1052,20 @@ def plot_cca_profiles(data,var='current_speed_',month=None,percentile=None,retur
     percentile: is the percentile associated with the worst case scenario
     return_period: a return-period e.g., 10 for a 10-yr return period
     distrbution, method, and threshold: only used if retrun_period is specified
+    max: by default, the highest curve will always be the worst-case scenario 
     unit_var: is a string with the unit of the variable var, default is m/s
     unit_lev: is a string with the units of the vertical levels, default is m
     output_file: name of the figure file
     with the dimensions (vertical levels of the profile, vertical level of the worst case scenario)
+    Function written by clio-met
     """
     if ((percentile is None) and (return_period is None)):
         print('Please specify either a percentile or a return period in years')
         sys.exit()
     if not(percentile is None):
         lev,woca,cca=cca_profiles(data,var=var,month=month,percentile=percentile)
     if not(return_period is None):
-        lev,woca,cca=cca_profiles(data,var=var,month=month,return_period=return_period,distribution=distribution,method=method,threshold=threshold)
+        lev,woca,cca=cca_profiles(data,var=var,month=month,return_period=return_period,distribution=distribution,method=method,threshold=threshold,max=max)
     import matplotlib as mpl
     n_lines = len(lev)
     cmap = mpl.colormaps['jet_r']

metocean_stats/stats/extreme.py

@@ -1326,7 +1326,7 @@ def get_joint_2D_contour(data=pd.DataFrame,var1='hs', var2='tp', periods=[50,100
         })
     return contours, data_2D
 
-def cca_profiles(data,var='current_speed_',month=None,percentile=None,return_period=None,distribution='GUM',method='default',threshold='default'):
+def cca_profiles(data,var='current_speed_',month=None,percentile=None,return_period=None,distribution='GUM',method='default',threshold='default',max='worst_case'):
     import sys
     """
     This function calculates the CCA profiles for a specific percentile or return period
@@ -1335,7 +1335,8 @@ def cca_profiles(data,var='current_speed_',month=None,percentile=None,return_per
     month: gives the month of interest (e.g. January or Jan), default (None) takes all months
     percentile: is the percentile associated with the worst case scenario
     return_period: a return-period e.g., 10 for a 10-yr return period
-    distrbution, method, and threshold: only used if retrun_period is specified
+    distrbution, method, and threshold: only used if return_period is specified
+    max: by default, if some points of the profile exceed the worst-case curve, they will be brought back down to the worst case value (if not wanted, write 'None')
     output: list of vertical levels used, 1-d array for the worst case scenario (the percentiles or return values), the array with the profiles
     with the dimensions (vertical levels of the profile, vertical level of the worst case scenario)
     Function written by clio-met based on 'Turkstra models of current profiles by Winterstein, Haver and Nygaard (2009)'
@@ -1389,15 +1390,18 @@ def cca_profiles(data,var='current_speed_',month=None,percentile=None,return_per
         nlevels=np.where(np.isnan(wcs))[0][0]
     # Calculate the current at the other depths when curr_ref = percentile or return-period value  
     cca_prof=np.zeros((nlevels,nlevels))
-    for d in range(nlevels):
+    for d in range(nlevels): # Loop over the worst case level
         prof_others=df_sel.drop(columns=[list_col[d]]).to_numpy() # extraxt currents as a matrix for all levels except the worst-case one
         prof_wcd=df_sel[list_col[d]].to_numpy() # extract currents at the worst-case level
         list_ind=np.arange(0,nlevels,1).tolist()
         del list_ind[d] # remove the index of the worst-case level (index d)
         cca_prof[d,d]=wcs[d]
         i=0
-        for dd in list_ind:
+        for dd in list_ind: # Loop over the other levels
             cca_prof[dd,d]=np.mean(prof_others[:,i],axis=0)+np.corrcoef(prof_wcd,prof_others[:,i])[0][1]*np.std(prof_others[:,i],axis=0)*((wcs[d]-np.mean(prof_wcd))/np.std(prof_wcd))
+            if max=='worst_case':
+                if cca_prof[dd,d]>wcs[dd]:
+                    cca_prof[dd,d]=wcs[dd]
             i=i+1
     return levels[0:nlevels],wcs[0:nlevels],cca_prof
 

metocean_stats/tables/extreme.py

@@ -733,7 +733,7 @@ def table_storm_surge_for_rv_hs(data: pd.DataFrame, var_hs='HS',var_tp='TP',var_
     return df
 
 
-def table_cca_profiles(data,var='current_speed_',month=None,percentile=None,return_period=None,distribution='GUM',method='default',threshold='default',output_file='table_cca_profiles.csv'):
+def table_cca_profiles(data,var='current_speed_',month=None,percentile=None,return_period=None,distribution='GUM',method='default',threshold='default',max='worst_case',output_file='table_cca_profiles.csv'):
     """
     This function returns a table containing the CCA profiles for a specific percentile or return period
     data: dataframe
@@ -742,16 +742,18 @@ def table_cca_profiles(data,var='current_speed_',month=None,percentile=None,retu
     percentile: is the percentile associated with the worst case scenario
     return_period: a return-period e.g., 10 for a 10-yr return period
     distrbution, method, and threshold: only used if retrun_period is specified
+    max: by default, the highest curve will always be the worst-case scenario
     output_file: name of the csv file
     with the dimensions (vertical levels of the profile, vertical level of the worst case scenario)
+    Function written by clio-met
     """
     if ((percentile is None) and (return_period is None)):
         print('Please specify either a percentile or a return period in years')
         sys.exit()
     if not(percentile is None):
         lev,woca,cca=cca_profiles(data,var=var,month=month,percentile=percentile)
     if not(return_period is None):
-        lev,woca,cca=cca_profiles(data,var=var,month=month,return_period=return_period,distribution=distribution,method=method,threshold=threshold)
+        lev,woca,cca=cca_profiles(data,var=var,month=month,return_period=return_period,distribution=distribution,method=method,threshold=threshold,max=max)
     list_ind=[]
     columns=['Depth [m]']+[d for d in lev]+['Worst case']
     table = np.zeros((len(lev),len(lev)+2))