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

Issue:
Dev1
-Add plot and table functions for number of hours below threshold
-add documentation
-update Readme


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 020927d..15418a9 100644
--- a/README.md
+++ b/README.md
@@ -1,14 +1,38 @@
-# metocean-stats
+<h1>
+  met<img src="https://github.com/MET-OM/metocean-stats/blob/main/docs/files/logo.png" width="20" height="20" />cean-stats
+ </h1>
 
 [![Tests (python)](https://github.com/MET-OM/metocean-stats/actions/workflows/tests.yml/badge.svg)](https://github.com/MET-OM/metocean-stats/actions/workflows/tests.yml)
 [![Documentation Status](https://readthedocs.org/projects/metocean-stats/badge/?version=latest)](https://metocean-stats.readthedocs.io/en/latest/?badge=latest)	[![Documentation Status](https://readthedocs.org/projects/metocean-stats/badge/?version=latest)](https://metocean-stats.readthedocs.io/en/latest/?badge=latest)
 
-## Purpose
-Metocen stastistics
+## Purpose 🎯
+Metocean-stats is a tool for comprehensive statistics and visualization of metocean data (wind, waves, ocean currents, tide levels, air and water temperature, sea ice, and more). The tool is also compatible with WEkEO Jupyter Lab, allowing seamless integration and use.
 
-## Quick Installation
+![metocean-stats](https://github.com/MET-OM/metocean-stats/blob/main/docs/files/readme_plots.png)
 
+⚠️ Please note that while metocean-stats is designed to provide accurate statistics and visualizations of metocean data, the user is ultimately responsible for correctly using the tool and interpreting the results. The tool is not liable for any errors or actions taken based on the information provided by the tool. Users should cross-verify the results and use them with discretion.
+
+## Quick Installation 💻
+
+To get started with metocean-stats, you can install it with pip or conda:
+
+```shell
 $ pip install metocean-stats 
+```
+
+or
+
+```shell
+$ conda install conda-forge::metocean-stats
+```
+
+## Docs 📚
+Looking for more details? Our documentation (under development) provides information about the package, its features, and how to use it. Check it out [here](https://metocean-stats.readthedocs.io/en/latest/index.html).
+
+Let's dive into the world of metocean data! 🌊🏄‍♀️
+
+<img src="https://github.com/MET-OM/metocean-stats/blob/main/docs/files/animation.gif" width="400" height="400" />
+
+
+
 
-## Docs
-Information about the package can be found at the [documentation](https://metocean-stats.readthedocs.io/en/latest/index.html).
diff --git a/docs/files/animation.gif b/docs/files/animation.gif
new file mode 100644
index 0000000..a007732
Binary files /dev/null and b/docs/files/animation.gif differ
diff --git a/docs/files/logo.png b/docs/files/logo.png
new file mode 100644
index 0000000..2efc4a1
Binary files /dev/null and b/docs/files/logo.png differ
diff --git a/docs/files/mean_air_temperature_map.png b/docs/files/mean_air_temperature_map.png
new file mode 100644
index 0000000..cacfc12
Binary files /dev/null and b/docs/files/mean_air_temperature_map.png differ
diff --git a/docs/files/nb_hour_below_thr.csv b/docs/files/nb_hour_below_thr.csv
new file mode 100644
index 0000000..60b95d2
--- /dev/null
+++ b/docs/files/nb_hour_below_thr.csv
@@ -0,0 +1,11 @@
+HS,Minimum,Mean,Maximum
+<1,771,1227,1812
+<2,3783,4454,5025
+<3,6036,6582,7086
+<4,7260,7717,8055
+<5,7971,8279,8472
+<6,8397,8538,8670
+<7,8583,8676,8745
+<8,8694,8734,8778
+<9,8736,8756,8784
+<10,8742,8763,8784
diff --git a/docs/files/nb_hour_below_thr.png b/docs/files/nb_hour_below_thr.png
new file mode 100644
index 0000000..6d7b6fa
Binary files /dev/null and b/docs/files/nb_hour_below_thr.png differ
diff --git a/docs/files/readme_plots.png b/docs/files/readme_plots.png
new file mode 100644
index 0000000..1760da1
Binary files /dev/null and b/docs/files/readme_plots.png differ
diff --git a/docs/files/weather_window_thresholds.csv b/docs/files/weather_window_thresholds.csv
new file mode 100644
index 0000000..ea36f6e
--- /dev/null
+++ b/docs/files/weather_window_thresholds.csv
@@ -0,0 +1,5 @@
+Operation duration [h],HS<0.5m,HS<1m,HS<2m
+6,130.66,16.63,2.2
+12,185.78,20.07,3.26
+24,433.91,31.09,5.73
+48,702.67,68.74,11.06
diff --git a/docs/index.rst b/docs/index.rst
index 4b9e351..e428af9 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -356,6 +356,54 @@ Monthly Weather Window Plot
 .. image:: files/NORA10_monthly_weather_window4_12_plot.png
    :width: 500
 
+Number Of Hours Per Year Below A Threshold Plot
+---------------------------
+
+.. code-block:: python
+
+   plots.plot_nb_hours_below_threshold(
+       df, 
+       var='HS', 
+       thr_arr=(np.arange(0.05,20.05,0.05)).tolist(),
+       output_file='number_hours_per_year.png'
+   )
+
+.. image:: files/nb_hour_below_thr.png
+   :width: 500
+
+Number Of Hours Per Year Below A Threshold Table
+---------------------------
+
+.. code-block:: python
+
+   tables.table_nb_hours_below_threshold(
+       df,
+       var='HS', 
+       threshold=[1,2,3,4,5,6,7,8,9,10],
+       output_file='number_hours_per_year.csv'
+   )
+
+.. csv-table:: 
+   :header-rows: 1
+   :file: files/nb_hour_below_thr.csv
+
+All-Year Round Weather Window For Hs Under A Threshold Table
+---------------------------
+
+.. code-block:: python
+
+   tables.table_weather_window_thresholds(
+       df,
+       var='HS', 
+       threshold=[0.5,1,2],
+       op_duration=[6,12,24,48],
+       output_file='weather_window_thresholds.csv'
+   )
+
+.. csv-table:: 
+   :header-rows: 1
+   :file: files/weather_window_thresholds.csv
+
 Monthly Return Periods Table
 ----------------------------
 
@@ -1585,5 +1633,34 @@ Plot extreme wind at 100 m height based on NORA3 data:
 .. image:: files/extreme_wind_map100m.png
   :width: 500
 
+Plot mean 2-m air temperature based on NORA3 data:
+
+.. code-block:: python
+
+   plot_mean_air_temperature_map(product='NORA3', 
+                                 title='Mean 2-m air temperature 1991-2020 (NORA3)', 
+                                 set_extent=[-25,-10,60.5,68], 
+                                 unit='degC')
+
+.. image:: files/mean_air_temperature_map.png
+  :width: 500
+
+
+Auxiliary Functions
+===================
+
+.. code-block:: python
+
+   from metocean_stats.stats.aux_funcs import *
+
+Convert lat/lon coordinates from degrees/minutes/seconds to decimals:
+
+.. code-block:: python
+
+   lat = degminsec_to_decimals(60,30,00)
+   
+returns lat = 60.5
+
+
 .. toctree::
-   :maxdepth: 1
+   :maxdepth: 1
\ No newline at end of file
diff --git a/examples/nb_hour_below_thr.csv b/examples/nb_hour_below_thr.csv
new file mode 100644
index 0000000..60b95d2
--- /dev/null
+++ b/examples/nb_hour_below_thr.csv
@@ -0,0 +1,11 @@
+HS,Minimum,Mean,Maximum
+<1,771,1227,1812
+<2,3783,4454,5025
+<3,6036,6582,7086
+<4,7260,7717,8055
+<5,7971,8279,8472
+<6,8397,8538,8670
+<7,8583,8676,8745
+<8,8694,8734,8778
+<9,8736,8756,8784
+<10,8742,8763,8784
diff --git a/examples/nb_hour_below_thr.png b/examples/nb_hour_below_thr.png
new file mode 100644
index 0000000..6d7b6fa
Binary files /dev/null and b/examples/nb_hour_below_thr.png differ
diff --git a/metocean_stats/maps/maps.py b/metocean_stats/maps/maps.py
index 59818e0..7bdda16 100644
--- a/metocean_stats/maps/maps.py
+++ b/metocean_stats/maps/maps.py
@@ -128,6 +128,7 @@ def plot_extreme_wave_map(return_period=50, product='NORA3', title='empty title'
     gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, x_inline=False, y_inline=False, linewidth=0.5, color='black', alpha=1, linestyle='--')
     gl.right_labels = False
     gl.top_labels = False
+    gl.rotate_labels=False
 
     cb = plt.colorbar(hb, ax=ax, orientation='vertical', shrink=0.7, pad=0.05)
     cb.set_label(ds['hs'].attrs.get('standard_name', 'hs') + ' [' + ds['hs'].attrs.get('units', 'm') + ']', fontsize=14)
@@ -186,6 +187,7 @@ def plot_extreme_wind_map(return_period=50, product='NORA3', z=0, title='empty t
     gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, x_inline=False, y_inline=False, linewidth=0.5, color='black', alpha=1, linestyle='--')
     gl.right_labels = False
     gl.top_labels = False
+    gl.rotate_labels=False
 
     cb = plt.colorbar(hb, ax=ax, orientation='vertical', shrink=0.7, pad=0.05)
     cb.set_label(ds['wind_speed'].attrs.get('standard_name', 'wind_speed') + ' [' + ds['wind_speed'].attrs.get('units', 'm/s') + ']', fontsize=14)
@@ -195,3 +197,75 @@ def plot_extreme_wind_map(return_period=50, product='NORA3', z=0, title='empty t
         plt.title(title, fontsize=16)
     plt.savefig(output_file, dpi=300)
     return fig
+
+# Function to plot mean air temperature map
+def plot_mean_air_temperature_map(product='NORA3', title='empty title', set_extent=[0, 30, 52, 73], unit='degC', mask_land=True, output_file='mean_air_temperature_map.png'):
+    """
+    Plots an extreme wind speed map based on the specified return level and dataset.
+
+    Parameters:
+    - product: str, optional, default='NORA3'
+        The dataset to use for the wind speed data. Currently only 'NORA3' is supported.
+    - title: str, optional, default='empty title'
+        The title to display on the map.
+    - set_extent: list, optional, default=[0, 30, 52, 73]
+        The map extent specified as [lon_min, lon_max, lat_min, lat_max].
+    - unit: str, optional, default='degC' (Celsius degrees), option 'K' for Kelvin
+    - mask_land. bool, optional, default=True for masked land
+    - output_file: str, optional, default='extreme_wind_map.png'
+        The filename to save the plot.
+
+    Returns:
+    - fig: matplotlib.figure.Figure
+        The figure object containing the plot.
+    """    
+    if product == 'NORA3':
+        ds = xr.open_dataset(f'https://thredds.met.no/thredds/dodsC/nora3_subset_stats/atm/CF_Overall_Mean_AirTemperature2m_1991_2020.nc')
+        #standard_lon, standard_lat, _ = get_transformed_coordinates(ds, 'x', 'y', projection_type='lambert_conformal')
+        hs_flat = ds['air_temperature_2m'].isel(time=0).values.flatten()
+        lon_flat = ds['longitude'].values.flatten()
+        lat_flat = ds['latitude'].values.flatten()
+    else:
+        print(product, 'is not available')
+        return
+
+    #mask = ~np.isnan(hs_flat)
+    #hs_flat = hs_flat[mask]
+    #lon_flat = standard_lon.flatten()[mask]
+    #lat_flat = standard_lat.flatten()[mask]
+    if unit=='degC':
+        vmin=-14
+        vmax=14
+    elif unit=='K':
+        hs_flat=hs_flat+273.15
+        vmin=260
+        vmax=290
+    else:
+        print('Chosen unit not correct, should be K or degC')
+        return
+
+    cmap = ListedColormap(plt.cm.get_cmap('RdYlBu_r', int((vmax-vmin)/1))(np.linspace(0, 1, int((vmax-vmin)/1))))
+
+    fig = plt.figure(figsize=(9, 10))
+    ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=set_extent[0]+(set_extent[1]-set_extent[0])/2))
+    ax.coastlines(resolution='10m', zorder=3)
+    if mask_land:
+        ax.add_feature(cfeature.BORDERS, linestyle=':', zorder=3)
+
+    hb = ax.hexbin(lon_flat, lat_flat, C=hs_flat, gridsize=180, cmap=cmap, vmin=vmin, vmax=vmax, edgecolors='white', transform=ccrs.PlateCarree(), reduce_C_function=np.mean, zorder=1)
+    ax.add_feature(cfeature.LAND, color='darkkhaki', zorder=2)
+
+    ax.set_extent(set_extent, crs=ccrs.PlateCarree())
+    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, x_inline=False, y_inline=False, linewidth=0.5, color='black', alpha=1, linestyle='--')
+    gl.right_labels = False
+    gl.top_labels = False
+    gl.rotate_labels=False
+
+    cb = plt.colorbar(hb, ax=ax, orientation='vertical', shrink=0.7, pad=0.05)
+    cb.set_label(ds['air_temperature_2m'].attrs.get('standard_name', 'air_temperature') + ' [' + unit + ']', fontsize=14)
+    plt.tight_layout()
+
+    if title:
+        plt.title(title, fontsize=16)
+    plt.savefig(output_file, dpi=300)
+    return fig
diff --git a/metocean_stats/plots/extreme.py b/metocean_stats/plots/extreme.py
index aeb7bad..2a67f6a 100644
--- a/metocean_stats/plots/extreme.py
+++ b/metocean_stats/plots/extreme.py
@@ -189,7 +189,7 @@ def plot_multi_diagnostic_return_levels(data, var,
                                                 block_size= str(interval_day) + "D")
         df_emp_rl.loc[:,'return_periods'] = df_emp_rl.loc[:,'return_periods']/(365.24/interval_day)
         df_emp_rl = df_emp_rl.loc[df_emp_rl.loc[:,'return_periods'] >= 1.0,:]
-                                                
+
     # Initialize plot and fill in empirical return levels
     fig, ax = plt.subplots()
 
@@ -200,8 +200,12 @@ def plot_multi_diagnostic_return_levels(data, var,
         # period greater than 1 year
         df_emp_rl_cut = df_emp_rl[df_emp_rl['return_periods'] >= 1]
 
+        #ax.scatter(df_emp_rl_cut['return_levels'],
+        #           df_emp_rl_cut['return_periods'],
+        #           marker="o", s=20, lw=1,
+        #           facecolor="k", edgecolor="w", zorder=20)
         ax.scatter(df_emp_rl_cut['return_levels'],
-                   df_emp_rl_cut['return_periods'],
+                   df_emp_rl_cut['prob_non_exceedance'],
                    marker="o", s=20, lw=1,
                    facecolor="k", edgecolor="w", zorder=20)
 
@@ -235,8 +239,11 @@ def plot_multi_diagnostic_return_levels(data, var,
 
         # Plot (return levels, return periods) lines corresponding
         if yaxis == 'prob':
+            #ax.plot(df_model_rl_tmp_cut['return_levels'],
+            #        df_model_rl_tmp_cut.index,
+            #        label=dist)
             ax.plot(df_model_rl_tmp_cut['return_levels'],
-                    df_model_rl_tmp_cut.index,
+                    df_model_rl_tmp_cut['prob_non_exceedance'],
                     label=dist)
 
         elif yaxis == 'rp':
@@ -251,19 +258,23 @@ def plot_multi_diagnostic_return_levels(data, var,
     # of non-exceedance or return period
     if yaxis == 'prob':
         ax.set_ylabel("Probability of non-exceedance")
-        list_yticks = [1/0.9, 2, 10]\
-                    + [10**i for i in range(2, 5) if max(periods) > 10**i]
-        if max(periods) > max(list_yticks):
-            list_yticks = list_yticks + [max(periods)]
-
+        #list_yticks = [1/0.9, 2, 10]\
+        #            + [10**i for i in range(2, 5) if max(periods) > 10**i]
+        #if max(periods) > max(list_yticks):
+        #    list_yticks = list_yticks + [max(periods)]
+        #ax.set_yticks(list_yticks)
+        #ax.set_yticklabels([round(1-1/rp,5) for rp in list_yticks])
+        list_yticks = [10**-1,10**(np.log10(0.5)),10**0]
+        list_ytickslabels = ['0.1','0.5','1']
         ax.set_yticks(list_yticks)
-        ax.set_yticklabels([round(1-1/rp,5) for rp in list_yticks])
+        ax.set_yticklabels(list_ytickslabels)
+        ax.set_ylim(0.05,1.02)
 
     elif yaxis == 'rp':
-        ax.set_ylabel("Return period")
+        ax.set_ylabel("Return period [yr]")
 
-    ax.set_xlabel("Return levels")
-    plt.legend()
+    ax.set_xlabel("Return levels "+var)
+    plt.legend(loc='lower right')
     plt.tight_layout()
 
     # Save the plot if a path is given
diff --git a/metocean_stats/plots/general.py b/metocean_stats/plots/general.py
index ee2b360..46aec96 100644
--- a/metocean_stats/plots/general.py
+++ b/metocean_stats/plots/general.py
@@ -379,3 +379,34 @@ def plot_tidal_levels(data, var='tide',start_time=None , end_time=None ,output_f
 
     return fig
 
+
+def plot_nb_hours_below_threshold(df,var='hs',thr_arr=(np.arange(0.05,20.05,0.05)).tolist(),output_file='plot_name.png'):
+    # Inputs
+    # 1) ds
+    # 2) var: a string
+    # 3) thr_arr: list of thresholds (should be a lot for smooth curve)
+    # 4) String with filename without extension
+    nbhr_arr=nb_hours_below_threshold(df,var,thr_arr)
+    years=df.index.year.to_numpy()
+    years_unique=np.unique(years)
+    yr1=int(years_unique[0])
+    yr2=int(years_unique[-1])
+    del df,years,years_unique
+    from matplotlib.patches import Patch
+    fig, ax = plt.subplots(figsize=(10, 8))
+    ax.fill_between(thr_arr,np.min(nbhr_arr,axis=1),np.max(nbhr_arr,axis=1),color='lightgray')
+    ax.plot(thr_arr,np.mean(nbhr_arr,axis=1),linewidth=2,color='k')
+    ax.grid(axis='both', color='gray',linestyle='dashed')
+    ax.set_ylim(0,9000)
+    ax.set_xlim(0,20)
+    ax.set_xlabel(var+' threshold [m]',fontsize=20)
+    ax.set_ylabel('Number of hours',fontsize=20)
+    ax.tick_params(axis='both', which='major', labelsize=18)
+    ax.set_title('Mean number of hours per year with '+var+' < threshold\n'+str(yr1)+'-'+str(yr2),fontsize=20)
+    legend_elements = [Patch(facecolor='lightgray', edgecolor=None, label='Min-Max range')]
+    ax.legend(handles=legend_elements, loc='lower right', prop={'size': 20})
+    plt.savefig(output_file, dpi=150, bbox_inches='tight')
+    plt.close()
+    del ax
+    return fig
+
diff --git a/metocean_stats/stats/aux_funcs.py b/metocean_stats/stats/aux_funcs.py
index 304b4c2..95a9120 100644
--- a/metocean_stats/stats/aux_funcs.py
+++ b/metocean_stats/stats/aux_funcs.py
@@ -496,3 +496,12 @@ def detrend_ts(df, column_name='tide'):
     import scipy.signal as signal
     df[column_name] = signal.detrend(df[column_name])
     return df
+
+
+def degminsec_to_decimals(degrees,minutes,seconds):
+    if degrees<=0:
+        loc_decimals=degrees-(minutes/60)-(seconds/3600)
+    else:
+        loc_decimals=degrees+(minutes/60)+(seconds/3600)
+    return loc_decimals
+
diff --git a/metocean_stats/stats/general.py b/metocean_stats/stats/general.py
index 533ac54..757c576 100644
--- a/metocean_stats/stats/general.py
+++ b/metocean_stats/stats/general.py
@@ -434,7 +434,7 @@ def calculate_monthly_weather_window(data: pd.DataFrame, var: str,threshold=5, w
     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],month=month ,threshold=threshold,op_duration=window_size,timestep=3)
+        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:
@@ -515,18 +515,17 @@ def calculate_weather_window(data: pd.DataFrame, var: str,threshold=5, window_si
     #breakpoint()
     return mean, p10, p50,p90
 
-def weather_window_length(time_series,month,threshold,op_duration,timestep):
-    print(month)
-    # time_series: input timeseries (numpy array)
+def weather_window_length(time_series,threshold,op_duration,timestep,month=None):
+    # time_series: input timeseries
     # threshold over which operation is possible (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 an array with all weather windows duration in hours
     month_ts = time_series.index.month
     ts_mask=np.where(time_series<threshold,1,0)
     od=int(op_duration/timestep)
     ts=np.zeros((len(ts_mask)-od+1))
-    mon=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()
@@ -553,13 +552,41 @@ def weather_window_length(time_series,month,threshold,op_duration,timestep):
     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
-    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)
+    if month==None:
+        mean = np.mean(wt1)
+        p10 = np.percentile(wt1,10)
+        p50 = np.percentile(wt1,50)
+        p90 = np.percentile(wt1,90)
+    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)
     return mean, p10, p50, p90
 
 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)
-    return  surge_min, surge_max
\ No newline at end of file
+    return  surge_min, surge_max
+
+
+def nb_hours_below_threshold(df,var,thr_arr):
+    thr_arr=np.array(thr_arr)
+    delta_t=(df.index.to_series().diff().dropna().dt.total_seconds()/3600).mean()
+    years=df.index.year.to_numpy()
+    years_unique=np.unique(years)
+    # Calculate the number of hours with hs below a threshold
+    nbhr_arr=np.zeros((len(thr_arr),len(years_unique)))
+    j=0
+    for yr in years_unique:
+        df1=df[df.index.year == yr]
+        df1=df1[var]
+        i=0
+        for t in thr_arr:
+            nbhr_arr[i,j]=len(df1[(df1 < t)])*delta_t
+            i=i+1
+        j=j+1
+        del df1
+    del i,j
+    del years,years_unique,delta_t
+    return nbhr_arr
diff --git a/metocean_stats/tables/general.py b/metocean_stats/tables/general.py
index 6b565b3..267dc02 100644
--- a/metocean_stats/tables/general.py
+++ b/metocean_stats/tables/general.py
@@ -11,7 +11,7 @@
 from ..stats.general import *
 
 
-def scatter_diagram(data: pd.DataFrame, var1: str, step_var1: float, var2: str, step_var2: float, output_file):
+def scatter_diagram(data: pd.DataFrame, var1: str, step_var1: float, var2: str, step_var2: float,number_of_events=False, output_file='scatter_diagram.csv'):
     """
     The function is written by dung-manh-nguyen and KonstantinChri.
     Plot scatter diagram (heatmap) of two variables (e.g, var1='hs', var2='tp')
@@ -21,6 +21,11 @@ def scatter_diagram(data: pd.DataFrame, var1: str, step_var1: float, var2: str,
      """
 
     sd = calculate_scatter(data, var1, step_var1, var2, step_var2)
+    start_date = data.index.min()
+    end_date = data.index.max() 
+    # Calculate the difference in years
+    number_of_years = (end_date - start_date).days / 365.25
+    dt = (data.index.to_series().diff().dropna().dt.total_seconds() / 3600).mean()
 
     # Convert to percentage
     tbl = sd.values
@@ -53,6 +58,13 @@ def scatter_diagram(data: pd.DataFrame, var1: str, step_var1: float, var2: str,
     #cols.insert(0,var1+' / '+var2 )
     rows = rows[::-1]
     tbl = tbl[::-1,:]
+    
+    if number_of_events==True:
+        #number of events 
+        tbl = np.flip(sd.values,axis=0)
+    else:
+        pass
+
     dfout = pd.DataFrame(data=np.round(tbl,2), index=rows, columns=cols)
     if output_file.split('.')[-1]=='csv':
         dfout.to_csv(output_file,index_label=var1+'/'+var2)
@@ -586,3 +598,63 @@ def table_max_min_water_level(data: pd.DataFrame, var_total_water_level: str,var
     df.to_csv(output_file, index=True)
     
     return df
+
+
+def table_nb_hours_below_threshold(df,var='hs',threshold=[0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,9,10,12.5,15,17.5,20],output_file='table_hours_below_thresh.csv'):
+    # Inputs
+    # 1) df
+    # 2) var: a string
+    # 3) thr_arr: list of thresholds (should be a lot for smooth curve)
+    # 4) thresholds_chosen: list of thresholds to be included in the table
+    # 5) output_file: String with filename without extension
+    thr_arr=(np.arange(0.05,20.05,0.05)).tolist()
+    nbhr_arr=nb_hours_below_threshold(df,var,thr_arr)
+    # Create file
+    threshold=np.array(threshold)
+    arr1=np.zeros((len(threshold),3))
+    rows=[]
+    for j in range(len(threshold)):
+        # Needs to add/subtract 0.001 because of precision problem
+        t=np.where((thr_arr>=threshold[j]-0.001) & (thr_arr<=threshold[j]+0.001))[0]
+        arr1[j,0]=np.min(nbhr_arr[t,:])
+        arr1[j,1]=np.round(np.mean(nbhr_arr[t,:]),0)
+        arr1[j,2]=np.max(nbhr_arr[t,:])
+        del t
+        rows.append('<='+str(threshold[j]))
+    # Convert numpy array to dataframe
+    df_out = pd.DataFrame()
+    df_out[var] = ['<'+str(threshold[j]) for j in range(len(threshold))]
+    #df_out = pd.DataFrame(data=np.round(arr1,0), index=rows, columns=cols)
+    df_out['Minimum'] = np.round(arr1[:,0],0).astype('int')
+    df_out['Mean'] = np.round(arr1[:,1],0).astype('int')
+    df_out['Maximum'] = np.round(arr1[:,2],0).astype('int')
+    # Save to csv format
+    df_out.to_csv(output_file, index=False)
+    return df_out
+
+def table_weather_window_thresholds(ds,var,threshold,op_duration=[6,12,24,48],output_file='table_ww.csv'):
+    # ds: dataframe
+    # var: string with var name
+    # threshold: list of thresholds appropriate for the variable
+    # op_duration: list, default is 6h, 12, 24h, and 48h
+    # output_file: string with file name
+    ds = ds[var]
+    timestep = (ds.index.to_series().diff().dropna().dt.total_seconds()/3600).mean()
+    arr_o=np.zeros((len(op_duration),len(threshold)))
+    t1=0
+    for t in threshold:
+        op1=0
+        for op in op_duration:
+            arr_o[op1,t1],p1,p2,p3 = weather_window_length(ds,threshold=t,op_duration=op,timestep=timestep,month=None)
+            op1=op1+1
+            del p1,p2,p3
+        t1=t1+1
+    # Convert numpy array to dataframe
+    cols=[var+'<'+str(th)+'m' for th in threshold]
+    df_tmp1 = pd.DataFrame(data=op_duration, columns=['Operation duration [h]'])
+    df_tmp2 = pd.DataFrame(data=np.round(arr_o,2), columns=cols)
+    df_out = pd.concat([df_tmp1, df_tmp2], axis=1, sort=True)
+    del df_tmp1,df_tmp2
+    # Save to csv format
+    df_out.to_csv(output_file, index=False)
+    return df_out


Test patch:
diff --git a/tests/test_maps.py b/tests/test_maps.py
index e55f7da..cd8c9bc 100644
--- a/tests/test_maps.py
+++ b/tests/test_maps.py
@@ -45,3 +45,18 @@
 #     else:
 #         raise FileNotFoundError(f"{output_file} was not created")
 #     return
+
+#def test_plot_mean_air_temperature_map():
+#     output_file = 'test_temp.png'
+#     fig = maps.plot_mean_air_temperature_map(product='NORA3', title='Mean 2-m air temperature 1991-2020 (NORA3)', set_extent=[-25, -10, 60.5, 68], unit='degC', output_file=output_file)
+#     if fig.dpi == 100.0:
+#         pass
+#     else:
+#         raise ValueError("FigValue is not correct")
+#    
+#     # Check if file is created
+#     if os.path.isfile(output_file):
+#         os.remove(output_file)
+#     else:
+#         raise FileNotFoundError(f"{output_file} was not created")
+#     return
\ No newline at end of file
diff --git a/tests/test_plots.py b/tests/test_plots.py
index c927ad9..d92e088 100644
--- a/tests/test_plots.py
+++ b/tests/test_plots.py
@@ -284,4 +284,12 @@ def test_plot_monthly_return_periods_cur_pot(ds=ds_ocean):
 #def test_predict_ts_LSTM(ds=ds):
 #    ml_stats.predict_ts(ts_origin=ds.data,var_origin=['hs','tp','Pdir'],ts_train=ds.data.loc['2000-01-01':'2000-01-10'],var_train=['hs'], model='LSTM')
 
-  
\ No newline at end of file
+def test_plot_nb_hours_below_threshold(ds=ds):
+    output_file = 'test_plot_nb_hr_below_t.png'
+    fig = plots.plot_nb_hours_below_threshold(ds,var='HS',thr_arr=(np.arange(0.05,20.05,0.05)).tolist(),output_file=output_file)
+    if os.path.exists(output_file):
+        os.remove(output_file)
+    if fig.dpi == 100.0:
+        pass
+    else:
+        raise ValueError("FigValue is not correct")
\ No newline at end of file
diff --git a/tests/test_tables.py b/tests/test_tables.py
index 18f0bcd..fb0c716 100644
--- a/tests/test_tables.py
+++ b/tests/test_tables.py
@@ -400,5 +400,27 @@ def test_table_max_min_water_level(ds=ds):
         os.remove(output_file)
     if df.shape == (2, 4):
         pass
+    else:
+        raise ValueError("Shape is not correct")
+    
+
+def test_table_nb_hours_below_threshold(ds=ds):
+    output_file = 'table_nb_hr_below_t.csv'
+    df = tables.table_nb_hours_below_threshold(ds,var='HS',threshold=[1,2,3,4,5,6,7,8,9,10],output_file=output_file)
+    if os.path.exists(output_file):
+        os.remove(output_file)
+    if df.shape == (10, 4):
+        pass
+    else:
+        raise ValueError("Shape is not correct")
+    
+
+def test_table_weather_window_thresholds(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)
+    if os.path.exists(output_file):
+        os.remove(output_file)
+    if df.shape == (4, 4):
+        pass
     else:
         raise ValueError("Shape is not correct")
\ No newline at end of file