Merge branch 'main' into tmp

environment.yml

@@ -3,6 +3,7 @@ channels:
   - conda-forge
 dependencies:
   - python=>3.9
+  - xarray
   - matplotlib>=3.1
   - numpy>=1.17
   - pandas

examples/example_NORA3_stats.py

@@ -7,7 +7,7 @@
 # Import data from thredds.met.no and save it as csv
 #ds.import_data(save_csv=True)
 # Load data from local file
-ds.load_data('tests/data/'+ds.datafile)
+ds.load_data('/home/birgitterf/dev/github/metocean-stats/tests/data/'+ds.datafile)
 
 # Generate Statistics
 general_stats.scatter_diagram(data=ds.data, var1='hs', step_var1=1, var2='tp', step_var2=1, output_file='scatter_hs_tp_diagram.png')
@@ -24,4 +24,5 @@
 rl_am = extreme_stats.return_levels_annual_max(data=ds.data, var='hs', periods=[20,50,100,1000],method='GEV',output_file='return_levels_GEV.png')
 
 # Profile Statistics
-mean_prof = 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='wind_profile.png')
+mean_prof = 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=[5,95], output_file='wind_profile.png')
+alfa = profile_stats.profile_shear(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],z=[20,50], perc=[5,95], output_file='wind_profile_shear.png')

metocean_stats/stats/profile_stats.py

@@ -1,6 +1,7 @@
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
 
 def mean_profile(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='mean_wind_profile.png'):  
     """
@@ -16,20 +17,67 @@ def mean_profile(data, vars = ['wind_speed_10m','wind_speed_20m','wind_speed_50m
     df = data[vars]
     nparr=df.to_numpy()
     mean_values = np.mean(nparr,axis=0) 
+    #breakpoint()
 
     if output_file==False:
         pass
     else:
         fig, ax = plt.subplots(figsize=(4,8))
-        plt.fill_betweenx(np.array(height_levels),np.percentile(nparr,perc[0],axis=0),np.percentile(nparr,perc[1],axis=0),color='lightgray')
-        plt.plot(mean_values,np.array(height_levels),'k-',linewidth=2.5)
+        plt.fill_betweenx(np.array(height_levels),np.percentile(nparr,perc[0],axis=0),np.percentile(nparr,perc[1],axis=0),color='lightgray',label=f'{perc[0]}-{perc[1]}%')
+        plt.plot(mean_values,np.array(height_levels),'k-',linewidth=2.5,label='Mean wind speed')
         plt.ylim(np.min(height_levels),np.max(height_levels))
         plt.xlabel('wind speed [m/s]')
         plt.ylabel('height [m]')
         plt.grid()
         plt.title(output_file.split('.')[0],fontsize=18)
         plt.tight_layout()
+        plt.legend()
         plt.savefig(output_file)
         plt.close()
     return mean_values
 
+def profile_shear(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], z=[20,250], perc = [25,75], output_file='shear_distribution.png'):  
+    """ This function is written by birgitterf.
+    It calculates and plot the shear between two height levels as a histogram 
+    with the median value in black and and lower/higher percentile in gray.
+    data = dataframe
+    vars = list of variables in dataframe e.g., ['wind_speed_10m','wind_speed_20m','wind_speed_50m']
+    height_levels = [10,20,50] 
+    z = interval to calculate the shear over e.g. [10,50]
+    perc = lower, higher percentile e.g., [25,75] 
+    output_file = 'shear_histogram.png' or False for no 
+    return shear values alfa as numpy array.
+    """
+    df = data[vars]
+    nparr = df.to_numpy()
+    H = np.array(height_levels)
+
+    i_upper = np.arange(len(H))[H==z[1]]
+    i_lower = np.arange(len(H))[H==z[0]]
+    if len(i_lower)==1:
+        i_lower = i_lower[0]
+        i_upper = i_upper[0]
+    else:
+        raise Exception
+
+    a = np.log(nparr[:,i_upper]/nparr[:,i_lower])
+    b = np.log(z[1]/z[0])
+    alfa = a/b
+
+    if output_file==False:
+        pass
+    else:
+        fig, ax = plt.subplots(figsize=(8,4))
+        N, bins, range = plt.hist(alfa,50)
+        locs, _ = plt.yticks() 
+        print(locs)
+        plt.yticks(locs,np.round(100*(locs/len(alfa)),1))
+        plt.fill_betweenx([0,np.max(N)],np.percentile(alfa,perc[0]),np.percentile(alfa,perc[1]),alpha=0.5,color='lightgray',label=f'{perc[0]}-{perc[1]}%')
+        plt.plot(np.ones(2)*np.median(alfa),[0,np.max(N)],color='black',label='median')
+        #ax.yaxis.set_major_formatter(mtick.PercentFormatter(xmax=1))
+        plt.xlabel(f'Shear between {z[0]} m and {z[1]} m')
+        plt.ylabel('Density [%]')
+        plt.legend()
+        plt.savefig(output_file)
+        plt.close()
+    return alfa

tests/test_stats.py

@@ -5,6 +5,7 @@
 # Define TimeSeries-object for NORA3
 ds = ts.TimeSeries(lon=1.32, lat=53.324,start_time='2000-01-01', end_time='2010-12-31' , product='NORA3_wind_wave')
 # Import data from thredds.met.no and save it as csv
+#ds.load_data('tests/data/'+ds.datafile)
 ds.load_data('tests/data/'+ds.datafile)
 
 def test_scatter_diagram(ds=ds):
@@ -42,10 +43,11 @@ def test_return_levels_idm(ds=ds):
 
 def test_mean_profile(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)
-
-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')
+
+def test_profile_shear(ds=ds):
+    profile_stats.profile_shear(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], z=[20,250], perc = [25,75], output_file=False)
+
+
 
 def test_diagnostic_return_level_plot_multi(ds=ds):
     extreme_stats.diagnostic_return_level_plot_multi(data=ds.data, 

version.py

@@ -1,4 +1,4 @@
-__version__ = "0.1.5"
+__version__ = "0.1.6"
 
 
 def git_describe():