Weather windows with multiple (max 3) variables

.github/workflows/tests.yml

@@ -11,7 +11,7 @@ on:
 
 jobs:
   tests:
-    name: pytests (${{ matrix.os }}, Miniforge3)
+    name: pytests (${{ matrix.os }} - ${{ matrix.python-version }}, Miniforge3)
     runs-on: ${{ matrix.os }}-latest
     defaults:
       run:
@@ -40,7 +40,7 @@ jobs:
           python setup.py bdist_wheel
 
       - name: Upload wheels
-        if: matrix.python-version == '3.10' && matrix.os == 'ubuntu'
+        if: matrix.python-version == '3.12' && matrix.os == 'ubuntu'
         uses: actions/upload-artifact@v4
         with:
           name: wheels
@@ -57,7 +57,7 @@ jobs:
           name: wheels
       - uses: conda-incubator/setup-miniconda@v3
         with:
-          python-version: '3.10'
+          python-version: '3.12'
           miniforge-version: latest
           miniforge-variant: Miniforge3
       - name: Publish to PyPi

docs/index.rst

@@ -344,20 +344,22 @@ Joint Distribution Hs-Tp Multi Plot with binned Var3
 .. image:: files/Hs.Tp.joint.distribution.multi.binned.var3.png
    :width: 500
 
-Monthly Weather Window Plot
----------------------------
+Monthly Weather Window Multiple Variables Plot (Waiting time in Days)
+----------------------------------------------
 
 .. code-block:: python
 
    plots.plot_monthly_weather_window(
        df, 
-       var='HS', 
-       threshold=4, 
-       window_size=12, 
-       output_file='NORA10_monthly_weather_window4_12_plot.png'
+       var=['W10','HS','TP'], 
+       threshold=[12.35,2,8], 
+       window_size=62,
+       timestep=3,
+       add_table=True,
+       output_file='NORA10_monthly_weather_window_mv_plot.png'
    )
 
-.. image:: files/NORA10_monthly_weather_window4_12_plot.png
+.. image:: files/NORA10_monthly_weather_window_mv_plot.png
    :width: 500
 
 Number Of Hours Per Year Below A Threshold Plot

environment.yml

@@ -2,7 +2,7 @@ name: metocean-stats
 channels:
   - conda-forge
 dependencies:
-  - python=>3.9
+  - python>=3.10
   - xarray
   - matplotlib>=3.1
   - numpy>=2.0

metocean_stats/plots/general.py

@@ -293,8 +293,8 @@ def plot_directional_stats(data: pd.DataFrame, var: str, step_var: float, var_di
     if output_file != "": plt.savefig(output_file)
     return fig
 
-def plot_monthly_weather_window(data: pd.DataFrame, var: str,threshold=5, window_size=12,add_table=True, output_file: str = 'monthly_weather_window_plot.png'):
-    results_df = stats.calculate_monthly_weather_window(data=data, var=var, threshold=threshold, window_size=window_size)
+def plot_monthly_weather_window(data: pd.DataFrame, var: str,threshold=5, window_size=12, timestep=3, add_table=True, output_file: str = 'monthly_weather_window_plot.png'):
+    results_df = tables.table_monthly_weather_window(data=data, var=var, threshold=threshold, window_size=window_size, timestep=timestep)
     # Plot the results
     fig, ax = plt.subplots(figsize=(12, 6))
     results_df.T.plot(marker='o')
@@ -351,6 +351,50 @@ def plot_monthly_weather_window(data: pd.DataFrame, var: str,threshold=5, window
 #     return df2 
 
 
+def plot_monthly_weather_window_MultipleVariables(data: pd.DataFrame, var: str,threshold=[5], window_size=12, timestep=3, add_table=True, output_file: str = 'monthly_weather_window_plot.png'):
+    # var is a list of variables (max 3) as well as thresholds (one for each variable)
+    # adjusted by clio-met
+    results_df = tables.table_monthly_weather_window_MultipleVariables(data=data, var=var, threshold=threshold, window_size=window_size, timestep=timestep)
+    # Plot the results
+    fig, ax = plt.subplots(figsize=(12, 6))
+    results_df.T.plot(marker='o')
+    lines = results_df.T.plot(marker='o')
+    title=[]
+    for i in range(len(var)):
+        title.append(var[i]+' < '+str(threshold[i]))
+    title=', '.join(title)
+    title=title+' for ' + str(window_size)+' hours'
+    plt.title(title)
+    plt.xlabel('Month')
+    plt.ylabel('Duration [days]')
+    plt.legend()
+    plt.grid(True)
+    plt.xticks(rotation=45)
+
+    if add_table:
+        # Get legend colors and labels
+        legend_colors = [line.get_color() for line in lines.get_lines()]
+        # Add the table of results_df under the plot
+        plt.xticks([])
+        plt.xlabel('')
+        plt.legend('',frameon=False)
+        cell_text = []
+        for row in range(len(results_df)):
+            cell_text.append(results_df.iloc[row].values)
+        table = plt.table(cellText=cell_text, colLabels=results_df.columns, rowLabels=results_df.index, loc='bottom')
+        table.auto_set_font_size(False)
+        table.set_fontsize(10)
+        table.scale(1, 1.5)
+        # Make the blue color value
+        cell_dict = table.get_celld()
+        for i in range(1,len(legend_colors)+1):
+            cell_dict[(i, -1)].set_facecolor(legend_colors[i-1])
+    plt.tight_layout()
+    if output_file != "": plt.savefig(output_file)
+
+    return fig, table
+
+
 def plot_profile_stats(data,var=['W10','W50','W80','W100','W150'], z=[10, 50, 80, 100, 150],reverse_yaxis=False, output_file='stats_profile.png'):
     df = tables.table_profile_stats(data=data, var=var, z=z, output_file=None)
     df = df.drop(['Std.dev', 'Max Speed Event'],axis=1)

metocean_stats/stats/general.py

@@ -420,24 +420,6 @@ def scatter(df,var1,var2,location,regression_line,qqplot=True):
     return fig
 
 
-
-
-def calculate_monthly_weather_window(data: pd.DataFrame, var: str,threshold=5, window_size=12,output_file: str = None):
-    results = []
-    months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
-    for month in range(1, 13):
-        avg_duration, p10, p50, p90 = weather_window_length(time_series=data[var],threshold=threshold,op_duration=window_size,timestep=3,month=month)
-        results.append((p10,p50, avg_duration, p90))
-    results_df = pd.DataFrame(results, columns=['P10', 'P50', 'Mean', 'P90'], index=months).T.round(2)
-    if output_file:
-        # Save results to CSV
-        results_df.to_csv('monthly_weather_window_results.csv')
-
-    return results_df
-
-
-
-
 def calculate_weather_window(data: pd.DataFrame, var: str,threshold=5, window_size=12):
     """
     Calculate the mean and percentiles of consecutive periods where a variable
@@ -556,6 +538,78 @@ def weather_window_length(time_series,threshold,op_duration,timestep,month=None)
         p90 = np.percentile(wt1[mon_s0==month],90)
     return mean, p10, p50, p90
 
+
+def weather_window_length_MultipleVariables(df,vars,threshold,op_duration,timestep,month=None):
+    """
+    Usage:
+    df: input df containing timeseries for different variables
+    vars: list of strings: variables to consider (up to 3)
+    threshold: list of floats: thresholds to use for each variable included (same unit as timeseries)
+    op_duration: duration of operation in hours
+    timestep: time resolution of time_series in hours
+    month: default is all year
+    Returns specific weather windows duration in hours
+    Generalization of weather_window_length to multiple conditions
+    Modified by clio-met
+    """
+    month_ts = df.index.month
+    if (len(vars)!=len(threshold)):
+        print('Error: vars must be the same length as threshold')
+        print(sys.exit())
+    if len(vars)==1:
+        ts_mask=np.where(df[vars[0]]<threshold[0],1,0)
+    elif len(vars)==2:
+        ts_mask=np.where(((df[vars[0]]<threshold[0]) & (df[vars[1]]<threshold[1])),1,0)
+    elif len(vars)==3:
+        ts_mask=np.where(((df[vars[0]]<threshold[0]) & (df[vars[1]]<threshold[1]) & (df[vars[2]]<threshold[2])),1,0)
+    else:
+        print('Error: only 3 variables can be given')
+        sys.exit()
+    od=int(op_duration/timestep)
+    ts=np.zeros((len(ts_mask)-od+1))
+    for i in range(len(ts_mask)-od+1):
+        ts[i]=np.sum(ts_mask[i:i+od])
+    s0=np.where(ts==od)[0].tolist()
+    mon_s0=month_ts[0:s0[-1]+1]
+    wt=[]
+    for s in range(len(s0)):
+        if s0[s]==0:
+            wt.append(0)
+        elif s==0:
+            diff=s0[s]
+            a=0
+            wt.append(diff)
+            while (diff!=0):
+                diff=s0[s]-a-1
+                wt.append(diff*timestep)
+                a=a+1
+        else:
+            diff=s
+            a=0
+            while (diff!=0):
+                diff=s0[s]-s0[s-1]-a-1
+                wt.append(diff*timestep)
+                a=a+1
+    wt1=(np.array(wt)+op_duration)/24
+    # Note that we this subroutine, we stop the calculation of the waiting time
+    # at the first timestep of the last operating period found in the timeseries
+    if month is None:
+        mean = np.mean(wt1)
+        p10 = np.percentile(wt1,10)
+        p50 = np.percentile(wt1,50)
+        p90 = np.percentile(wt1,90)
+        p95 = np.percentile(wt1,95)
+        max = np.max(wt1)
+    else:
+        mean = np.mean(wt1[mon_s0==month])
+        p10 = np.percentile(wt1[mon_s0==month],10)
+        p50 = np.percentile(wt1[mon_s0==month],50)
+        p90 = np.percentile(wt1[mon_s0==month],90)
+        p95 = np.percentile(wt1[mon_s0==month],95)
+        max = np.max(wt1[mon_s0==month])
+    return mean, p10, p50, p90, p95, max
+
+
 def pressure_surge(df,var='MSLP'):
     surge_max = (min(df.MSLP)-np.mean(df.MSLP))*(-0.01)
     surge_min = (max(df.MSLP)-np.mean(df.MSLP))*(-0.01)

metocean_stats/tables/general.py

@@ -601,18 +601,34 @@ def monthly_directional_percentiles(
 
     return monthly_tables
 
-
-def table_monthly_weather_window(data: pd.DataFrame, var: str,threshold=5, window_size=12,output_file: str = None):
+# For one variable
+#def table_monthly_weather_window(data: pd.DataFrame, var: str,threshold=5, window_size=12, timestep=3, output_file: str = None):
+#    results = []
+#    months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
+#    for month in range(1, 13):
+#        avg_duration, p10, p50, p90 = stats.weather_window_length(time_series=data[var],month=month ,threshold=threshold,op_duration=window_size,timestep=timestep)
+#        results.append((p10,p50, avg_duration, p90))
+#    results_df = pd.DataFrame(results, columns=['P10', 'P50', 'Mean', 'P90'], index=months).T.round(2)
+#    if output_file:
+#        # Save results to CSV
+#        results_df.to_csv('monthly_weather_window_results.csv')
+#    
+#    return results_df
+
+# For multivariable
+def table_monthly_weather_window(data: pd.DataFrame, var: str, threshold: float, window_size=12, timestep=3, output_file: str = None):
+    # var should be a list of variables, and threshold should be a list of thresholds
+    # more outputs than table_monthly_weather_window
+    # Written by clio-met
     results = []
     months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
     for month in range(1, 13):
-        avg_duration, p10, p50, p90 = stats.weather_window_length(time_series=data[var],month=month ,threshold=threshold,op_duration=window_size,timestep=3)
-        results.append((p10,p50, avg_duration, p90))
-    results_df = pd.DataFrame(results, columns=['P10', 'P50', 'Mean', 'P90'], index=months).T.round(2)
+        avg_duration, p10, p50, p90, p95, max = stats.weather_window_length_MultipleVariables(data,vars=var,threshold=threshold,op_duration=window_size,timestep=timestep,month=month)
+        results.append((avg_duration, p10, p50, p90, p95, max))
+    results_df = pd.DataFrame(results, columns=['Mean', 'P10', 'P50', 'P90', 'P95', 'Max'], index=months).T.round(1)
     if output_file:
         # Save results to CSV
-        results_df.to_csv('monthly_weather_window_results.csv')
-
+        results_df.to_csv('monthly_weather_window_results_mv.csv')
     return results_df
 
 def table_profile_stats(data: pd.DataFrame, var: str, z=[10, 20, 30], var_dir=None, output_file='table_profile_stats.csv'):

tests/test_plots.py

@@ -66,8 +66,8 @@ def test_plot_joint_distribution_Hs_Tp(ds=ds):
         raise ValueError("FigValue is not correct")
 
 def test_plot_monthly_weather_window(ds=ds):
-    output_file = 'test_monthly_weather_window.png'
-    fig, table = plots.plot_monthly_weather_window(ds, var='HS', threshold=4, window_size=12, output_file=output_file)
+    output_file = 'test_monthly_weather_window_mv.png'
+    fig, table = plots.plot_monthly_weather_window(ds, var=['HS','TP'], threshold=[2,8], timestep=3, window_size=12, output_file=output_file)
     if os.path.exists(output_file):
         os.remove(output_file)
     if table._loc == 17:

tests/test_tables.py

@@ -411,14 +411,13 @@ def test_table_nb_hours_below_threshold(ds=ds):
         pass
     else:
         raise ValueError("Shape is not correct")
-
 
-def test_table_weather_window_thresholds(ds=ds):
+def test_table_monthly_weather_window(ds=ds):
     output_file = 'table_ww_threshold.csv'
-    df = tables.table_weather_window_thresholds(ds,var='HS',threshold=[0.5,1,2],op_duration=[6,12,24,48],output_file=output_file)
+    df = tables.table_monthly_weather_window(ds,var=['HS','TP'],threshold=[2,8],window_size=24,timestep=3,output_file=output_file)
     if os.path.exists(output_file):
         os.remove(output_file)
-    if df.shape == (4, 4):
+    if df.shape == (6, 12):
         pass
     else:
         raise ValueError("Shape is not correct")