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

Issue:
Increase fontsizes in plot_multi_diagnostic_return_levels_uncertainty



Discussion / hints:

Fix patch:
diff --git a/docs/files/Hs.Tp.joint.distribution.multi.binned.var3.png b/docs/files/Hs.Tp.joint.distribution.multi.binned.var3.png
new file mode 100644
index 0000000..ee78ff7
Binary files /dev/null and b/docs/files/Hs.Tp.joint.distribution.multi.binned.var3.png differ
diff --git a/docs/index.rst b/docs/index.rst
index e428af9..ca35362 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -340,6 +340,27 @@ Directional Joint Distribution Hs-Tp Parameter Table
    :header-rows: 1
    :file: files/dir_Hs_Tp_joint_param.csv
 
+
+Joint Distribution Hs-Tp Mulit Plot with binned Var3
+------------------------------
+
+.. code-block:: python
+
+   plots.plots.plot_multi_joint_distribution_Hs_Tp_var3(
+      df,
+      var_hs='HS',
+      var_tp='TP',
+      var3='W10',
+      var3_units='m/s',
+      periods=[100],
+      var3_bin=5,
+      threshold_min=100,
+      output_file='Hs.Tp.joint.distribution.multi.binned.var3.png'
+      )
+
+.. image:: files/Hs.Tp.joint.distribution.multi.binned.var3.png
+   :width: 500
+
 Monthly Weather Window Plot
 ---------------------------
 
diff --git a/examples/example_NORA.py b/examples/example_NORA.py
index e4a25cb..28abccb 100644
--- a/examples/example_NORA.py
+++ b/examples/example_NORA.py
@@ -11,64 +11,65 @@
     folder.mkdir(parents=True)
 
 # Map:
-maps.plot_points_on_map(lon=[3.35,3.10], lat=[60.40,60.90],label=['NORA3','NORKYST800'],bathymetry='NORA3',output_file=folder / 'map.png')
-maps.plot_extreme_wave_map(return_period=100, product='NORA3', title='100-yr return values Hs (NORA3)', set_extent = [0,30,52,73],output_file=folder / 'wave_100yrs.png')
-maps.plot_extreme_wind_map(return_period=100, product='NORA3',z=10, title='100-yr return values Wind at 100 m (NORA3)', set_extent = [0,30,52,73], output_file=folder / 'wind_100yrs.png')
+maps.plot_points_on_map(lon=[3.35,3.10], lat=[60.40,60.90],label=['NORA3','NORKYST800'],bathymetry='NORA3',output_file='map.png')
+maps.plot_extreme_wave_map(return_period=100, product='NORA3', title='100-yr return values Hs (NORA3)', set_extent = [0,30,52,73],output_file='wave_100yrs.png')
+maps.plot_extreme_wind_map(return_period=100, product='NORA3',z=10, title='100-yr return values Wind at 100 m (NORA3)', set_extent = [0,30,52,73], output_file='wind_100yrs.png')
 
 
 # Wind:
-plots.var_rose(df,var_dir='D10',var='W10',method='overall',max_perc=40,decimal_places=1, units='m/s',output_file=folder / 'wind_omni.png')
-plots.var_rose(df,var_dir='D10',var='W10',method='monthly',max_perc=40,decimal_places=1, units='m/s',output_file=folder / 'wind_monthly.png')
-plots.plot_directional_stats(df,var='W10',step_var=0.1,var_dir='D10',title = '', output_file=folder / 'directional_stats.png')
-plots.table_directional_non_exceedance(df,var='W10',step_var=2,var_dir='D10',output_file=folder / 'table_directional_non_exceedance.csv')
-plots.plot_monthly_stats(df,var='W10',title = 'Wind Speed at 10 m [m/s]', output_file=folder / 'monthly_stats.png')
-tables.table_monthly_non_exceedance(df,var='W10',step_var=2,output_file=folder / 'table_monthly_non_exceedance.csv')
-plots.plot_prob_non_exceedance_fitted_3p_weibull(df,var='W10',output_file=folder / 'prob_non_exceedance_fitted_3p_weibull.png')
-plots.plot_monthly_return_periods(df,var='W10',periods=[1, 10, 100],distribution='Weibull3P_MOM', units='m/s',output_file=folder / 'W10_monthly_extremes.png')
-plots.plot_directional_return_periods(df,var='W10',var_dir='D10',periods=[1, 10, 100],distribution='Weibull3P_MOM', units='m/s',adjustment='NORSOK',output_file=folder / 'W10_dir_extremes_Weibull_norsok.png')
-plots.plot_monthly_weather_window(df,var='W10',threshold=10, window_size=12,output_file=folder / 'NORA10_monthly_weather_window4_12_plot.png')
-plots.plot_scatter(df,var1='W10',var2='W100',var1_units='m/s',var2_units='m/s', title='Scatter',regression_line=True,qqplot=False,density=True,output_file=folder / 'scatter_plot.png')
-plots.plot_multi_diagnostic_return_levels(df, var='HS', periods=[10, 100],threshold=None,output_file=folder / 'plot_diagnostic_return_levels.png')
+plots.var_rose(df,var_dir='D10',var='W10',method='overall',max_perc=40,decimal_places=1, units='m/s',output_file='wind_omni.png')
+plots.var_rose(df,var_dir='D10',var='W10',method='monthly',max_perc=40,decimal_places=1, units='m/s',output_file='wind_monthly.png')
+plots.plot_directional_stats(df,var='W10',step_var=0.1,var_dir='D10',title = '', output_file='directional_stats.png')
+plots.table_directional_non_exceedance(df,var='W10',step_var=2,var_dir='D10',output_file='table_directional_non_exceedance.csv')
+plots.plot_monthly_stats(df,var='W10',title = 'Wind Speed at 10 m [m/s]', output_file='monthly_stats.png')
+tables.table_monthly_non_exceedance(df,var='W10',step_var=2,output_file='table_monthly_non_exceedance.csv')
+plots.plot_prob_non_exceedance_fitted_3p_weibull(df,var='W10',output_file='prob_non_exceedance_fitted_3p_weibull.png')
+plots.plot_monthly_return_periods(df,var='W10',periods=[1, 10, 100],distribution='Weibull3P_MOM', units='m/s',output_file='W10_monthly_extremes.png')
+plots.plot_directional_return_periods(df,var='W10',var_dir='D10',periods=[1, 10, 100],distribution='Weibull3P_MOM', units='m/s',adjustment='NORSOK',output_file='W10_dir_extremes_Weibull_norsok.png')
+plots.plot_monthly_weather_window(df,var='W10',threshold=10, window_size=12,output_file='NORA10_monthly_weather_window4_12_plot.png')
+plots.plot_scatter(df,var1='W10',var2='W100',var1_units='m/s',var2_units='m/s', title='Scatter',regression_line=True,qqplot=False,density=True,output_file='scatter_plot.png')
+plots.plot_multi_diagnostic_return_levels(df, var='HS', periods=[10, 100],threshold=None,output_file='plot_diagnostic_return_levels.png')
 
 # Waves:
-plots.plot_prob_non_exceedance_fitted_3p_weibull(df,var='HS',output_file=folder / 'prob_non_exceedance_fitted_3p_weibull.png')
-tables.scatter_diagram(df, var1='HS', step_var1=1, var2='TP', step_var2=1, output_file= str(folder / 'Hs_Tp_scatter.csv'))
-tables.table_var_sorted_by_hs(df, var='TP', var_hs='HS', output_file=folder / 'Tp_sorted_by_Hs.csv')
-tables.table_monthly_non_exceedance(df,var='HS',step_var=0.5,output_file=folder / 'Hs_table_monthly_non_exceedance.csv')
-plots.plot_monthly_stats(df,var='HS',show=['Maximum','P99','Mean'], title = 'Hs[m]', output_file=folder / 'Hs_monthly_stats.png')
-tables.table_directional_non_exceedance(df,var='HS',step_var=0.5,var_dir='DIRM',output_file=folder / 'table_directional_non_exceedance.csv')
-plots.plot_directional_stats(df,var='HS',step_var=0.5, var_dir='DIRM', title = '$H_s$[m]', output_file=folder / 'directional_stats.png')
-plots.plot_joint_distribution_Hs_Tp(df,var_hs='HS',var_tp='TP',periods=[1,10,100,1000], title='Hs-Tp joint distribution',output_file=folder / 'Hs.Tp.joint.distribution.png',density_plot=True)
-tables.table_monthly_joint_distribution_Hs_Tp_param(df,var_hs='HS',var_tp='TP',periods=[1,10,100,10000],output_file=folder / 'monthly_Hs_Tp_joint_param.csv')
-tables.table_directional_joint_distribution_Hs_Tp_param(df,var_hs='HS',var_tp='TP',var_dir='DIRM',periods=[1,10,100],output_file=folder / 'dir_Hs_Tp_joint_param.csv')
-plots.plot_monthly_weather_window(df,var='HS',threshold=4, window_size=12,output_file=folder / 'NORA10_monthly_weather_window4_12_plot.png')
-tables.table_monthly_return_periods(df,var='HS',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM', units='m',output_file=folder / 'HS_monthly_extremes_Weibull.csv')
-tables.table_directional_return_periods(df,var='HS',periods=[1, 10, 100, 10000], units='m',var_dir = 'DIRM',distribution='Weibull3P_MOM', adjustment='NORSOK' ,output_file=folder / 'directional_extremes_weibull.csv')
-plots.plot_monthly_return_periods(df,var='HS',periods=[1, 10, 100],distribution='Weibull3P_MOM', units='m',output_file=folder / 'HS_monthly_extremes.png')
-plots.plot_directional_return_periods(df,var='HS',var_dir='DIRM',periods=[1, 10, 100, 10000 ],distribution='GUM', units='m',output_file=folder / 'dir_extremes_GUM.png')
-plots.plot_directional_return_periods(df,var='HS',var_dir='DIRM',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM', units='m',adjustment='NORSOK',output_file=folder / 'dir_extremes_Weibull_norsok.png')
-tables.table_monthly_joint_distribution_Hs_Tp_return_values(df,var_hs='HS',var_tp='TP',periods=[1,10,100,10000],output_file=folder / 'monthly_Hs_Tp_joint_return_values.csv')
-tables.table_directional_joint_distribution_Hs_Tp_return_values(df,var_hs='HS',var_tp='TP',var_dir='DIRM',periods=[1,10,100,1000],adjustment='NORSOK',output_file=str(folder / 'directional_Hs_Tp_joint_return_values.csv'))
-tables.table_Hs_Tpl_Tph_return_values(df,var_hs='HS',var_tp='TP',periods=[1,10,100,10000],output_file=folder / 'hs_tpl_tph_return_values.csv')
-plots.plot_tp_for_given_hs(df, 'HS', 'TP',output_file=folder / 'tp_for_given_hs.png')
-tables.table_tp_for_given_hs(df, 'HS', 'TP',max_hs=20,output_file=folder / 'tp_for_given_hs.csv')
-tables.table_tp_for_rv_hs(df, var_hs='HS', var_tp='TP',periods=[1,10,100,10000],output_file=folder / 'tp_for_rv_hs.csv')
-tables.table_wave_induced_current(df, var_hs='HS',var_tp='TP',depth=200,ref_depth=200, spectrum = 'JONSWAP',output_file=folder / 'JONSWAP_wave_induced_current_depth200.csv')
-tables.table_wave_induced_current(df, var_hs='HS',var_tp='TP',depth=200,ref_depth=200, spectrum = 'TORSEHAUGEN',output_file=folder / 'TORSEHAUGEN_wave_induced_current_depth200.csv')
-tables.table_hs_for_given_wind(df, 'HS','W10', bin_width=2, max_wind=42, output_file= str(folder / 'table_perc_hs_for_wind.csv'))
-plots.plot_hs_for_given_wind(df, 'HS', 'W10',output_file=folder / 'hs_for_given_wind.png')
-tables.table_hs_for_rv_wind(df, var_wind='W10', var_hs='HS',periods=[1,10,100,10000],output_file=folder / 'hs_for_rv_wind.csv')
-tables.table_Hmax_crest_return_periods(df,var_hs='HS', var_tp='TP', depth=200, periods=[1, 10, 100,10000],sea_state='long-crested',output_file=folder / 'table_Hmax_crest_rp.csv')
-tables.table_directional_Hmax_return_periods(df,var_hs='HS', var_tp = 'TP',var_dir='DIRM', periods=[10, 100,10000],adjustment='NORSOK', output_file=folder / 'table_dir_Hmax_return_values.csv')
+plots.plot_prob_non_exceedance_fitted_3p_weibull(df,var='HS',output_file='prob_non_exceedance_fitted_3p_weibull.png')
+tables.scatter_diagram(df, var1='HS', step_var1=1, var2='TP', step_var2=1, output_file= str('Hs_Tp_scatter.csv'))
+tables.table_var_sorted_by_hs(df, var='TP', var_hs='HS', output_file='Tp_sorted_by_Hs.csv')
+tables.table_monthly_non_exceedance(df,var='HS',step_var=0.5,output_file='Hs_table_monthly_non_exceedance.csv')
+plots.plot_monthly_stats(df,var='HS',show=['Maximum','P99','Mean'], title = 'Hs[m]', output_file='Hs_monthly_stats.png')
+tables.table_directional_non_exceedance(df,var='HS',step_var=0.5,var_dir='DIRM',output_file='table_directional_non_exceedance.csv')
+plots.plot_directional_stats(df,var='HS',step_var=0.5, var_dir='DIRM', title = '$H_s$[m]', output_file='directional_stats.png')
+plots.plot_joint_distribution_Hs_Tp(df,var_hs='HS',var_tp='TP',periods=[1,10,100,1000], title='Hs-Tp joint distribution',output_file='Hs.Tp.joint.distribution.png',density_plot=True)
+tables.table_monthly_joint_distribution_Hs_Tp_param(df,var_hs='HS',var_tp='TP',periods=[1,10,100,10000],output_file='monthly_Hs_Tp_joint_param.csv')
+tables.table_directional_joint_distribution_Hs_Tp_param(df,var_hs='HS',var_tp='TP',var_dir='DIRM',periods=[1,10,100],output_file='dir_Hs_Tp_joint_param.csv')
+plots.plot_monthly_weather_window(df,var='HS',threshold=4, window_size=12,output_file='NORA10_monthly_weather_window4_12_plot.png')
+tables.table_monthly_return_periods(df,var='HS',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM', units='m',output_file='HS_monthly_extremes_Weibull.csv')
+tables.table_directional_return_periods(df,var='HS',periods=[1, 10, 100, 10000], units='m',var_dir = 'DIRM',distribution='Weibull3P_MOM', adjustment='NORSOK' ,output_file='directional_extremes_weibull.csv')
+plots.plot_monthly_return_periods(df,var='HS',periods=[1, 10, 100],distribution='Weibull3P_MOM', units='m',output_file='HS_monthly_extremes.png')
+plots.plot_directional_return_periods(df,var='HS',var_dir='DIRM',periods=[1, 10, 100, 10000 ],distribution='GUM', units='m',output_file='dir_extremes_GUM.png')
+plots.plot_directional_return_periods(df,var='HS',var_dir='DIRM',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM', units='m',adjustment='NORSOK',output_file='dir_extremes_Weibull_norsok.png')
+tables.table_monthly_joint_distribution_Hs_Tp_return_values(df,var_hs='HS',var_tp='TP',periods=[1,10,100,10000],output_file='monthly_Hs_Tp_joint_return_values.csv')
+tables.table_directional_joint_distribution_Hs_Tp_return_values(df,var_hs='HS',var_tp='TP',var_dir='DIRM',periods=[1,10,100,1000],adjustment='NORSOK',output_file=str('directional_Hs_Tp_joint_return_values.csv'))
+tables.table_Hs_Tpl_Tph_return_values(df,var_hs='HS',var_tp='TP',periods=[1,10,100,10000],output_file='hs_tpl_tph_return_values.csv')
+plots.plot_tp_for_given_hs(df, 'HS', 'TP',output_file='tp_for_given_hs.png')
+tables.table_tp_for_given_hs(df, 'HS', 'TP',max_hs=20,output_file='tp_for_given_hs.csv')
+tables.table_tp_for_rv_hs(df, var_hs='HS', var_tp='TP',periods=[1,10,100,10000],output_file='tp_for_rv_hs.csv')
+tables.table_wave_induced_current(df, var_hs='HS',var_tp='TP',depth=200,ref_depth=200, spectrum = 'JONSWAP',output_file='JONSWAP_wave_induced_current_depth200.csv')
+tables.table_wave_induced_current(df, var_hs='HS',var_tp='TP',depth=200,ref_depth=200, spectrum = 'TORSEHAUGEN',output_file='TORSEHAUGEN_wave_induced_current_depth200.csv')
+tables.table_hs_for_given_wind(df, 'HS','W10', bin_width=2, max_wind=42, output_file= str('table_perc_hs_for_wind.csv'))
+plots.plot_hs_for_given_wind(df, 'HS', 'W10',output_file='hs_for_given_wind.png')
+tables.table_hs_for_rv_wind(df, var_wind='W10', var_hs='HS',periods=[1,10,100,10000],output_file='hs_for_rv_wind.csv')
+tables.table_Hmax_crest_return_periods(df,var_hs='HS', var_tp='TP', depth=200, periods=[1, 10, 100,10000],sea_state='long-crested',output_file='table_Hmax_crest_rp.csv')
+tables.table_directional_Hmax_return_periods(df,var_hs='HS', var_tp = 'TP',var_dir='DIRM', periods=[10, 100,10000],adjustment='NORSOK', output_file='table_dir_Hmax_return_values.csv')
+plots.plot_multi_joint_distribution_Hs_Tp_var3(df,var_hs='HS',var_tp='TP',var3='W10',var3_units='m/s',periods=[100],var3_bin=5,threshold_min=100,output_file='Hs.Tp.joint.distribution.multi.binned.var3.png')
 
 
 # Air Temperature:
-plots.plot_monthly_return_periods(df,var='T2m',periods=[1, 10, 100],distribution='GUM_L',method='minimum', units='°C',output_file=folder / 'T2m_monthly_extremes_neg.png')
-tables.table_monthly_return_periods(df,var='T2m',periods=[1, 10, 100],distribution='GUM_L', method='minimum' ,units='°C',output_file=folder / 'T2m_monthly_extremes_neg.csv')
-plots.plot_monthly_return_periods(df,var='T2m',periods=[1, 10, 100],distribution='GUM', method='maximum', units='°C',output_file=folder / 'T2m_monthly_extremes_pos.png')
-tables.table_monthly_return_periods(df,var='T2m',periods=[1, 10, 100],distribution='GUM', method='maximum' ,units='°C',output_file=folder / 'T2m_monthly_extremes_pos.csv')
-plots.plot_monthly_stats(df,var='T2m',show=['Minimum','Mean','Maximum'], title = 'T2m', output_file=folder / 'T2m_monthly_stats.png')
-tables.table_monthly_non_exceedance(df,var='T2m',step_var=0.5,output_file=folder / 'T2m_table_monthly_non_exceedance.csv')
+plots.plot_monthly_return_periods(df,var='T2m',periods=[1, 10, 100],distribution='GUM_L',method='minimum', units='°C',output_file='T2m_monthly_extremes_neg.png')
+tables.table_monthly_return_periods(df,var='T2m',periods=[1, 10, 100],distribution='GUM_L', method='minimum' ,units='°C',output_file='T2m_monthly_extremes_neg.csv')
+plots.plot_monthly_return_periods(df,var='T2m',periods=[1, 10, 100],distribution='GUM', method='maximum', units='°C',output_file='T2m_monthly_extremes_pos.png')
+tables.table_monthly_return_periods(df,var='T2m',periods=[1, 10, 100],distribution='GUM', method='maximum' ,units='°C',output_file='T2m_monthly_extremes_pos.csv')
+plots.plot_monthly_stats(df,var='T2m',show=['Minimum','Mean','Maximum'], title = 'T2m', output_file='T2m_monthly_stats.png')
+tables.table_monthly_non_exceedance(df,var='T2m',step_var=0.5,output_file='T2m_table_monthly_non_exceedance.csv')
 
 
 # Currents:
@@ -80,63 +81,63 @@
 #ds_all = ds_all.dropna(how='all')
 
 
-#plots.plot_monthly_return_periods(ds_ocean,var='current_speed_0m',periods=[1, 10, 100],distribution='Weibull3P_MOM',method='POT',threshold='P99', units='m/s',output_file=folder / 'csp0m_monthly_extremes.png')
-#plots.var_rose(df, f'current_direction_{depth}',f'current_speed_{depth}',max_perc=30,decimal_places=2, units='m/s', method='monthly', output_file=folder / 'monthly_rose.png')
-#plots.var_rose(df,f'current_direction_{depth}',f'current_speed_{depth}',max_perc=30,decimal_places=2, units='m/s', method='overall', output_file=folder / 'overall_rose.png')
+#plots.plot_monthly_return_periods(ds_ocean,var='current_speed_0m',periods=[1, 10, 100],distribution='Weibull3P_MOM',method='POT',threshold='P99', units='m/s',output_file='csp0m_monthly_extremes.png')
+#plots.var_rose(df, f'current_direction_{depth}',f'current_speed_{depth}',max_perc=30,decimal_places=2, units='m/s', method='monthly', output_file='monthly_rose.png')
+#plots.var_rose(df,f'current_direction_{depth}',f'current_speed_{depth}',max_perc=30,decimal_places=2, units='m/s', method='overall', output_file='overall_rose.png')
 #plots.plot_monthly_stats(df,var1=f'current_speed_{depth}',show=['Mean','P99','Maximum'], title = f'Current[m/s], depth:{depth}', output_file=f'current{depth}_monthly_stats.png')
 #plots.plot_directional_stats(df,var1=f'current_speed_{depth}',var_dir=f'current_direction_{depth}',step_var1=0.05,show=['Mean','P99','Maximum'], title = f'Current[m/s], depth:{depth}', output_file=f'current{depth}_dir_stats.png')
 #tables.table_directional_return_periods(df,var=f'current_speed_{depth}',periods=[1, 10, 100, 10000], units='m/s',var_dir = f'current_direction_{depth}',distribution='Weibull', adjustment='NORSOK' ,output_file=f'directional_extremes_weibull_current{depth}.csv')
 #tables.table_monthly_return_periods(df,var=f'current_speed_{depth}',periods=[1, 10, 100, 10000], units='m/s',distribution='Weibull3P_MOM',method='POT',threshold='P99',output_file=f'monthly_extremes_weibull_current{depth}.csv')
-#df = tables.table_monthly_return_periods(df,var='HS',periods=[1, 10, 100, 10000], units='m',distribution='Weibull3P_MOM',method='POT',threshold='P99',output_file=folder / 'HS_monthly_extremes_Weibull_POT_P99.csv')
-#df = tables.table_monthly_return_periods(df,var='HS',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM', units='m',output_file=folder / 'HS_monthly_extremes_Weibull.csv')
-#df = tables.table_directional_return_periods(df,var='HS',var_dir='DIRM',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM', units='m',output_file=folder / 'HS_dir_extremes_Weibull.csv')
-#df = tables.table_directional_return_periods(df,var='HS',var_dir='DIRM',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM',method='POT',threshold='P99', units='m',output_file=folder / 'HS_dir_extremes_Weibull_POT.csv')
-# plots.plot_directional_return_periods(df,var='HS',var_dir='DIRM',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM', units='m',adjustment='NORSOK',method='POT',threshold='P99',output_file=folder / 'dir_extremes_Weibull_norsok.png')
-# df = tables.table_profile_return_values(df,var=['W10','W50','W80','W100','W150'], z=[10, 50, 80, 100, 150], periods=[1, 10, 100, 10000], output_file=folder / 'RVE_wind_profile.csv')
-#fig = plots.plot_profile_return_values(ds_ocean,var=['current_speed_' + d for d in depth], z=[float(d[:-1]) for d in depth], periods=[1, 10, 100, 10000],reverse_yaxis=True, output_file=folder / 'RVE_current_profile.png')
-#df = tables.table_current_for_given_wind(ds_all, var_curr='current_speed_0m', var_wind='W10', bin_width=2, max_wind=42, output_file=folder / 'table_perc_current_for_wind.csv')
-#plots.plot_current_for_given_wind(ds_all, var_curr='current_speed_0m', var_wind='W10',max_wind=40 ,output_file=folder / 'curr_for_given_wind.png')
-#df = tables.table_current_for_given_hs(ds_all, var_curr='current_speed_0m', var_hs='HS', bin_width=2, max_hs=20, output_file=folder / 'table_perc_current_for_Hs.csv')
+#df = tables.table_monthly_return_periods(df,var='HS',periods=[1, 10, 100, 10000], units='m',distribution='Weibull3P_MOM',method='POT',threshold='P99',output_file='HS_monthly_extremes_Weibull_POT_P99.csv')
+#df = tables.table_monthly_return_periods(df,var='HS',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM', units='m',output_file='HS_monthly_extremes_Weibull.csv')
+#df = tables.table_directional_return_periods(df,var='HS',var_dir='DIRM',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM', units='m',output_file='HS_dir_extremes_Weibull.csv')
+#df = tables.table_directional_return_periods(df,var='HS',var_dir='DIRM',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM',method='POT',threshold='P99', units='m',output_file='HS_dir_extremes_Weibull_POT.csv')
+# plots.plot_directional_return_periods(df,var='HS',var_dir='DIRM',periods=[1, 10, 100, 10000],distribution='Weibull3P_MOM', units='m',adjustment='NORSOK',method='POT',threshold='P99',output_file='dir_extremes_Weibull_norsok.png')
+# df = tables.table_profile_return_values(df,var=['W10','W50','W80','W100','W150'], z=[10, 50, 80, 100, 150], periods=[1, 10, 100, 10000], output_file='RVE_wind_profile.csv')
+#fig = plots.plot_profile_return_values(ds_ocean,var=['current_speed_' + d for d in depth], z=[float(d[:-1]) for d in depth], periods=[1, 10, 100, 10000],reverse_yaxis=True, output_file='RVE_current_profile.png')
+#df = tables.table_current_for_given_wind(ds_all, var_curr='current_speed_0m', var_wind='W10', bin_width=2, max_wind=42, output_file='table_perc_current_for_wind.csv')
+#plots.plot_current_for_given_wind(ds_all, var_curr='current_speed_0m', var_wind='W10',max_wind=40 ,output_file='curr_for_given_wind.png')
+#df = tables.table_current_for_given_hs(ds_all, var_curr='current_speed_0m', var_hs='HS', bin_width=2, max_hs=20, output_file='table_perc_current_for_Hs.csv')
 #ds['current_speed_0m'] = 0.05*ds['W10']
-#df = tables.table_current_for_given_wind(df, var_curr='current_speed_0m', var_wind='W10', bin_width=2, max_wind=42, output_file=folder / 'table_perc_current_for_wind.csv')
-#df = tables.table_current_for_given_hs(df, var_curr='current_speed_0m', var_hs='HS', bin_width=2, max_hs=20, output_file=folder / 'table_perc_current_for_Hs.csv')
+#df = tables.table_current_for_given_wind(df, var_curr='current_speed_0m', var_wind='W10', bin_width=2, max_wind=42, output_file='table_perc_current_for_wind.csv')
+#df = tables.table_current_for_given_hs(df, var_curr='current_speed_0m', var_hs='HS', bin_width=2, max_hs=20, output_file='table_perc_current_for_Hs.csv')
 
-#plots.plot_current_for_given_hs(ds_all, var_curr='current_speed_0m', var_hs='HS', max_hs=20, output_file=folder / 'curr_for_given_hs.png')
-#df = tables.table_extreme_current_profile_rv(ds_ocean, var=['current_speed_' + d for d in depth], z=[float(d[:-1]) for d in depth], periods=[1,100,1000],percentile=95, output_file=folder / 'table_extreme_current_profile_rv.png')
-#df = tables.table_profile_stats(ds_ocean, var=['current_speed_' + d for d in depth], z=[float(d[:-1]) for d in depth], var_dir=['current_direction_' + d for d in depth], output_file=folder / 'table_profile_stats.csv')
-#fig = plots.plot_profile_stats(ds_ocean,var=['current_speed_' + d for d in depth], z=[float(d[:-1]) for d in depth],reverse_yaxis=True, output_file=folder / 'stats_current_profile.png')
-#df = tables.table_current_for_rv_wind(ds_all, var_curr='current_speed_0m', var_wind='W10',periods=[1,10,100,10000],output_file=folder / 'Uc_for_rv_wind.csv')
-#df = tables.table_current_for_rv_hs(ds_all, var_curr='current_speed_0m', var_hs='HS',periods=[1,10,100,10000],output_file=folder / 'Uc_for_rv_hs.csv')
+#plots.plot_current_for_given_hs(ds_all, var_curr='current_speed_0m', var_hs='HS', max_hs=20, output_file='curr_for_given_hs.png')
+#df = tables.table_extreme_current_profile_rv(ds_ocean, var=['current_speed_' + d for d in depth], z=[float(d[:-1]) for d in depth], periods=[1,100,1000],percentile=95, output_file='table_extreme_current_profile_rv.png')
+#df = tables.table_profile_stats(ds_ocean, var=['current_speed_' + d for d in depth], z=[float(d[:-1]) for d in depth], var_dir=['current_direction_' + d for d in depth], output_file='table_profile_stats.csv')
+#fig = plots.plot_profile_stats(ds_ocean,var=['current_speed_' + d for d in depth], z=[float(d[:-1]) for d in depth],reverse_yaxis=True, output_file='stats_current_profile.png')
+#df = tables.table_current_for_rv_wind(ds_all, var_curr='current_speed_0m', var_wind='W10',periods=[1,10,100,10000],output_file='Uc_for_rv_wind.csv')
+#df = tables.table_current_for_rv_hs(ds_all, var_curr='current_speed_0m', var_hs='HS',periods=[1,10,100,10000],output_file='Uc_for_rv_hs.csv')
 
 
 # Sea temperature
-# df = tables.table_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='mean', output_file=folder / 'table_mean_temp_profile_monthly_stats.png')
-# df = tables.table_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='std.dev', output_file=folder / 'table_std_temp_profile_monthly_stats.png')
-# df = tables.table_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='minimum', output_file=folder / 'table_min_temp_profile_monthly_stats.png')
-# df = tables.table_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='maximum', output_file=folder / 'table_max_temp_profile_monthly_stats.png')
+# df = tables.table_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='mean', output_file='table_mean_temp_profile_monthly_stats.png')
+# df = tables.table_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='std.dev', output_file='table_std_temp_profile_monthly_stats.png')
+# df = tables.table_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='minimum', output_file='table_min_temp_profile_monthly_stats.png')
+# df = tables.table_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='maximum', output_file='table_max_temp_profile_monthly_stats.png')
 
-# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='mean',title='Mean Sea Temperature [°C]', output_file=folder / 'plot_mean_temp_profile_monthly_stats.png')
-# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='std.dev',title='St.Dev Sea Temperature [°C]', output_file=folder / 'plot_std_temp_profile_monthly_stats.png')
-# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='minimum',title='Min. Sea Temperature [°C]', output_file=folder / 'plot_min_temp_profile_monthly_stats.png')
-# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='maximum',title='Max. Sea Temperature [°C]', output_file=folder / 'plot_max_temp_profile_monthly_stats.png')
+# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='mean',title='Mean Sea Temperature [°C]', output_file='plot_mean_temp_profile_monthly_stats.png')
+# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='std.dev',title='St.Dev Sea Temperature [°C]', output_file='plot_std_temp_profile_monthly_stats.png')
+# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='minimum',title='Min. Sea Temperature [°C]', output_file='plot_min_temp_profile_monthly_stats.png')
+# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['temp_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='maximum',title='Max. Sea Temperature [°C]', output_file='plot_max_temp_profile_monthly_stats.png')
 
 # # Sainity:
-# df = tables.table_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='mean', output_file=folder / 'table_mean_sal_profile_monthly_stats.png')
-# df = tables.table_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='std.dev', output_file=folder / 'table_std_sal_profile_monthly_stats.png')
-# df = tables.table_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='minimum', output_file=folder / 'table_min_sal_profile_monthly_stats.png')
-# df = tables.table_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='maximum', output_file=folder / 'table_max_sal_profile_monthly_stats.png')
+# df = tables.table_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='mean', output_file='table_mean_sal_profile_monthly_stats.png')
+# df = tables.table_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='std.dev', output_file='table_std_sal_profile_monthly_stats.png')
+# df = tables.table_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='minimum', output_file='table_min_sal_profile_monthly_stats.png')
+# df = tables.table_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='maximum', output_file='table_max_sal_profile_monthly_stats.png')
 
-# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='mean',title='Mean Salinity [PSU]', output_file=folder / 'plot_mean_sal_profile_monthly_stats.png')
-# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='std.dev',title='St.Dev Salinity [PSU]', output_file=folder / 'plot_std_sal_profile_monthly_stats.png')
-# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='minimum',title='Min. Salinity [PSU]', output_file=folder / 'plot_min_sal_profile_monthly_stats.png')
-# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='maximum',title='Max. Salinity [PSU]', output_file=folder / 'plot_max_sal_profile_monthly_stats.png')
+# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='mean',title='Mean Salinity [PSU]', output_file='plot_mean_sal_profile_monthly_stats.png')
+# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='std.dev',title='St.Dev Salinity [PSU]', output_file='plot_std_sal_profile_monthly_stats.png')
+# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='minimum',title='Min. Salinity [PSU]', output_file='plot_min_sal_profile_monthly_stats.png')
+# fig = plots.plot_profile_monthly_stats(ds_ocean, var=['salt_' + d for d in depth], z=[float(d[:-1]) for d in depth], method='maximum',title='Max. Salinity [PSU]', output_file='plot_max_sal_profile_monthly_stats.png')
 
 
 # Water levels:
 #ds_tide = pd.read_csv('../tests/data/GTSM_test.csv',comment='#',index_col=0, parse_dates=True)
-#df = tables.table_tidal_levels(ds_tide, var='tide', output_file=folder / 'tidal_levels.csv')
-#fig = plots.plot_tidal_levels(ds_tide, var='tide',start_time='2010-01-01',end_time='2010-03-30', output_file=folder / 'tidal_levels.png')
-#df = tables.table_storm_surge_for_given_hs(ds_all, var_surge='zeta_0m', var_hs='HS', bin_width=1, max_hs=20, output_file=folder / 'table_perc_surge_for_Hs.csv')
-#fig = plots.plot_storm_surge_for_given_hs(df,var_surge='zeta_0m', var_hs='HS', max_hs=20, output_file=folder / 'surge_for_given_hs.png')
-#df = tables.table_extreme_total_water_level(df, var_hs='HS',var_tp='TP',var_surge='zeta_0m', var_tide='tide', periods=[100,10000], output_file=folder / 'table_extreme_total_water_level.csv')
-#df = tables.table_storm_surge_for_rv_hs(df, var_hs='HS',var_tp='TP',var_surge='zeta_0m', var_tide='tide', periods=[1,10,100,10000],depth=200, output_file=folder / 'table_storm_surge_for_rv_hs.csv')
+#df = tables.table_tidal_levels(ds_tide, var='tide', output_file='tidal_levels.csv')
+#fig = plots.plot_tidal_levels(ds_tide, var='tide',start_time='2010-01-01',end_time='2010-03-30', output_file='tidal_levels.png')
+#df = tables.table_storm_surge_for_given_hs(ds_all, var_surge='zeta_0m', var_hs='HS', bin_width=1, max_hs=20, output_file='table_perc_surge_for_Hs.csv')
+#fig = plots.plot_storm_surge_for_given_hs(df,var_surge='zeta_0m', var_hs='HS', max_hs=20, output_file='surge_for_given_hs.png')
+#df = tables.table_extreme_total_water_level(df, var_hs='HS',var_tp='TP',var_surge='zeta_0m', var_tide='tide', periods=[100,10000], output_file='table_extreme_total_water_level.csv')
+#df = tables.table_storm_surge_for_rv_hs(df, var_hs='HS',var_tp='TP',var_surge='zeta_0m', var_tide='tide', periods=[1,10,100,10000],depth=200, output_file='table_storm_surge_for_rv_hs.csv')
diff --git a/metocean_stats/plots/extreme.py b/metocean_stats/plots/extreme.py
index f0786a5..dc938cb 100644
--- a/metocean_stats/plots/extreme.py
+++ b/metocean_stats/plots/extreme.py
@@ -430,8 +430,9 @@ def plot_multi_diagnostic_return_levels_uncertainty(data, var,
         plt.legend(loc='lower right')
     elif yaxis == 'rp':
         ax.set_xscale('log')
-        ax.set_xlabel("Return period [yr]")
-        ax.set_ylabel("Return levels "+var)
+        ax.set_xlabel("Return period [yr]", fontsize=22)
+        ax.set_ylabel("Return levels "+var, fontsize=22)
+        ax.tick_params(axis='both', which='major', labelsize=20)
         plt.legend(loc='upper left')
 
     plt.tight_layout()
@@ -603,6 +604,64 @@ def plot_RVE_ALL(dataframe,var='hs',periods=np.array([1,10,100,1000]),distributi
        
     return 
 
+################
+
+def plot_multi_joint_distribution_Hs_Tp_var3(data,var_hs='hs',var_tp='tp',var3='W10',var3_units='m/s',periods=[100],var3_bin=5,threshold_min=100,output_file='Hs.Tp.joint.distribution.multi.binned.var3.png'):  
+    import matplotlib.pyplot as plt
+    from cycler import cycler
+
+    """
+    Plot joint distribution of Hs-Tp for a given return period for eached binned var3 data
+    Input:
+        var3: e.g., 'W10' 
+        var3_units: e.g., 'm/s'
+        var3_bin: sets the bin size, e.g., 5 
+        threshold_min: provide is the minimum number of datapoints in the dataset, e.g., 100
+    """
+    var3_name = var3 
+    df = data.dropna()
+    max_var3 = max(df[var3])
+    window = np.arange(0, max_var3 + var3_bin, var3_bin)
+    t3_ = []; h3_= []
+
+    fig, ax = plt.subplots(figsize=(8,6))
+    custom_cycler = cycler(color=['b', 'g', 'r', 'c', 'm', 'y', 'k'])
+    ax.set_prop_cycle(custom_cycler)
+
+    for i in range(len(window)-1):
+        var_3 = df[var3].where((df[var3] > window[i]) & (df[var3] <= window[i+1])).dropna() #w10 = (wind1,wind2]
+        if len(var_3)< threshold_min:
+            continue
+        hs = df[var_hs].where(var_3.notnull()).dropna()
+        tp = df[var_tp].where(var_3.notnull()).dropna()
+        dff = pd.DataFrame({'HS': hs, 'TP': tp, var3_name: var_3})
+        a1, a2, a3, b1, b2, b3, pdf_Hs, h, t3,h3,X,hs_tpl_tph = joint_distribution_Hs_Tp(data=dff,var_hs=var_hs,var_tp=var_tp,periods=periods)
+        t3_.append(t3)
+        h3_.append(h3)
+        linestyle = '-' if i % 2 == 0 else '--'
+        labels = str(np.round(window[i],2))+'-'+str(np.round(window[i+1],2))
+        plt.plot(t3[0],h3[0], linestyle=linestyle, label=var3_name+'$\\in$' + labels+' ['+var3_units+']')  
+
+    #include the whole dataset also:
+    a1, a2, a3, b1, b2, b3, pdf_Hs, h, t3,h3,X,hs_tpl_tph = joint_distribution_Hs_Tp(data=df,var_hs=var_hs,var_tp=var_tp,periods=periods)
+    labels = str(np.round(window[0],2))+'-'+str(np.round(window[-1],2))
+    plt.plot(t3[0],h3[0],label=var3_name+'$_{all}$'+'$\\in$ ' + labels+ ' ['+var3_units+']') 
+    plt.xlabel('Tp - Peak Period [s]')
+    plt.suptitle('Return period = '+ str(periods[0])+'-year')
+    plt.ylabel('Hs - Significant Wave Height [m]')
+    plt.grid()
+    plt.legend()
+    #plt.xlim([0,np.ceil(max(max(t3_[0])))]) #When the dataset is split 
+    #plt.ylim([0,np.ceil(max(max(h3_[-1])))]) #When the dataset is split
+    plt.xlim([0,np.ceil(max(max(t3)))]) #When the whole dataset is included 
+    plt.ylim([0,np.ceil(max(max(h3)))]) #When the whole dataset is included    
+    #plt.show()
+    plt.savefig(output_file,dpi=100,facecolor='white',bbox_inches='tight')
+
+    return fig
+
+##################
+
 def plot_joint_distribution_Hs_Tp(data,var_hs='hs',var_tp='tp',periods=[1,10,100,10000], title='Hs-Tp joint distribution',output_file='Hs.Tp.joint.distribution.png',density_plot=False):
     a1, a2, a3, b1, b2, b3, pdf_Hs, h, t3,h3,X,hs_tpl_tph = joint_distribution_Hs_Tp(data=data,var_hs=var_hs,var_tp=var_tp,periods=periods)
     df = data
@@ -640,8 +699,8 @@ def plot_joint_distribution_Hs_Tp(data,var_hs='hs',var_tp='tp',periods=[1,10,100
 
     for i in range(len(periods)):
         plt.plot(t3[i],h3[i],label=str(X[i])+'-year')  
-        
-   
+
+
     plt.plot(t_steepness,h_steepness,'k--',label='steepness')
     
     plt.plot(percentile50[0],percentile50[1],'g',label='Tp-mean',linewidth=5)
diff --git a/metocean_stats/stats/extreme.py b/metocean_stats/stats/extreme.py
index 156ed87..cc3b7e4 100644
--- a/metocean_stats/stats/extreme.py
+++ b/metocean_stats/stats/extreme.py
@@ -617,16 +617,18 @@ def joint_distribution_Hs_Tp(data,var_hs='hs',var_tp='tp',periods=[1,10,100,1000
     alpha = mu
     sigma = std
     
-    h = np.linspace(start=0.01, stop=30, num=1500)
+    #h = np.linspace(start=0.01, stop=30, num=1500)
+    h = np.linspace(start=0.01, stop=np.ceil(max(df[var_hs].values)*np.sqrt(np.pi)), num=1500)
     pdf_Hs1 = h*0
     pdf_Hs2 = h*0
     
     if 0 < mu < 5 : 
-    	pdf_Hs1 = 1/(np.sqrt(2*np.pi)*alpha*h)*np.exp(-(np.log(h)-sigma)**2/(2*alpha**2))
+      #Based on Moan et al. (2005), "Uncertainty of wave-induced response of marine structures due to long-term variation of extratropical wave conditions":
+    	pdf_Hs1 = 1/(np.sqrt(2*np.pi)*sigma*h)*np.exp(-(np.log(h)-alpha)**2/(2*sigma**2))
     else:
-    	param = stats.lognorm.fit(df.hs.values,) # shape, loc, scale
-    	pdf_lognorm = stats.lognorm.pdf(h, param[0], loc=param[1], scale=param[2])
-    	pdf_Hs1 = pdf_lognorm
+        param = stats.lognorm.fit(df.hs.values,) # shape, loc, scale
+        pdf_lognorm = stats.lognorm.pdf(h, param[0], loc=param[1], scale=param[2])
+        pdf_Hs1 = pdf_lognorm
     
     param = Weibull_method_of_moment(df.hs.values) #stats.weibull_min.fit(df.hs.values) # shape, loc, scale
     pdf_Hs2 = stats.weibull_min.pdf(h, param[0], loc=param[1], scale=param[2])
@@ -708,7 +710,7 @@ def joint_distribution_Hs_Tp(data,var_hs='hs',var_tp='tp',periods=[1,10,100,1000
     start = 1
     x = mean_hs[start:]
     y = variance_lnTp[start:]
-    parameters, covariance = curve_fit(Gauss4, x, y)
+    parameters, covariance = curve_fit(Gauss4, x, y)#,maxfev=10000)
     b1 = 0.005
     b2 = parameters[0]
     b3 = parameters[1]
@@ -729,7 +731,6 @@ def joint_distribution_Hs_Tp(data,var_hs='hs',var_tp='tp',periods=[1,10,100,1000
         pdf_Hs_Tp[i,:] = pdf_Hs[i]*f_Hs_Tp[i,:]
     
     interval = ((df.index[-1]-df.index[0]).days + 1)*24/df.shape[0] # in hours 
-
     t3 = []
     h3 = []
     X = []


Test patch:
diff --git a/tests/test_plots.py b/tests/test_plots.py
index 888e541..99503ac 100644
--- a/tests/test_plots.py
+++ b/tests/test_plots.py
@@ -304,4 +304,16 @@ def test_plot_multi_diagnostic_with_uncertainty(ds=ds):
     else:
         raise ValueError("FigValue is not correct")
     
-    
\ No newline at end of file
+
+def test_plot_multi_joint_distribution_Hs_Tp_var3(ds=ds):
+    output_file = 'test_mulit_joint_distribution_Hs_Tp_var3.png'
+    fig = plots.plot_multi_joint_distribution_Hs_Tp_var3(ds,var_hs='HS',var_tp='TP',var3='W10',var3_units='m/s',periods=[100],var3_bin=10,threshold_min=100,output_file=output_file)
+    if os.path.exists(output_file):
+        os.remove(output_file)
+#    if int(fig.axes[0].lines[0].get_xdata()[0]) == 3:
+    if fig.dpi == 100.0:
+        pass
+    else:
+        raise ValueError("FigValue is not correct")
+
+