Add files via upload

pyproject.toml

@@ -23,7 +23,7 @@ classifiers = [
 ]
 dynamic = ["version", "description"]
 
-dependencies = ["numpy", "pytest >=7", "scipy", "shapely", "xarray", "matplotlib"]
+dependencies = ["numpy", "pytest >=7", "scipy", "shapely", "xarray", "matplotlib","pandas"]
 
 
 [project.optional-dependencies]

src/IHSetBernabeu/IHSetBernabeu.py

@@ -2,44 +2,53 @@
 from scipy.interpolate import interp1d
 from shapely.geometry import LineString
 import xarray as xr
-from IHSetUtils import wMOORE
+import pandas as pd
+from IHSetUtils import wMOORE, Hs12Calc, depthOfClosure
 
 class cal_Bernabeu(object):
     """
     cal_Bernabeu
     
     Configuration to calibrate and run the Bernabeu profile.
     
-    This class reads input datasets, performs its calibration.
+    This class reads input datasets, calculates its parameters.
     """
-    def __init__(self, path):
-        self.path = path
-
-        # cfg = xr.open_dataset(path+'config.nc')
-        ens = xr.open_dataset(path+'ens.nc')
-        wav = xr.open_dataset(path+'wav.nc')
-
-        self.D50 = ens['D50'].values
-        self.dp = ens['d'].values
-        self.zp = ens['z'].values
-        self.CM = ens['CM_95'].values
-        self.Hs = wav['Hs50'].values
-        self.Tp = wav['Tp50'].values
-
-    def calibrate(self):
-        self.zp = self.zp - self.zp[0]
-        self.dd = self.dp - self.dp[0]
+    def __init__(self, path_prof, path_wav, Switch_Obs, Switch_Cal_DoC, **kwargs):
+        self.path_prof = path_prof
+        self.path_wav = path_wav
+        prof = pd.read_csv(path_prof)
 
-        # Profile with equidistant points
-        dp = np.linspace(0, self.dp[-1], 500).reshape(-1, 1)
-
-        interp_func = interp1d(self.dd, self.zp, kind="linear", fill_value="extrapolate")
-        zp = interp_func(dp)
-        zp = zp[1:]
-        dp = dp[1:]
+        self.Hs50 = kwargs['Hs50']
+        self.Tp50 = kwargs['Tp50']
+        self.D50 = kwargs['D50']
+        self.HTL = -kwargs['HTL']    
+        self.LTL = -kwargs['LTL']
+        self.CM = abs(self.HTL - self.LTL)
+        self.xm = np.linspace(kwargs['Xm'][0], kwargs['Xm'][1], 1000).reshape(-1, 1)
 
+        self.Switch_Obs = Switch_Obs              # Do you have profile data? (0: no; 1: yes)
+        if Switch_Obs == 1:
+            self.xp = prof.iloc[:, 0]
+            self.zp = prof.iloc[:, 1]
+            self.zp = abs(self.zp)
+            xp_inx = self.xp[(self.zp >= self.HTL)]
+            self.xp = self.xp - min(xp_inx)
+
+        self.Switch_Cal_DoC = Switch_Cal_DoC
+        if Switch_Cal_DoC == 1:                   # Calculate Depth of Closure if you have wave data [0: no; 1: yes]
+            wav = xr.open_dataset(path_wav)
+            Hs = wav['Hs'].values
+            Hs = Hs.reshape(-1, 1)
+            Tp = wav['Tp'].values
+            Tp = Tp.reshape(-1, 1)
+
+            H12,T12 = Hs12Calc(Hs,Tp)
+            self.DoC = depthOfClosure(H12,T12)
+            self.DoC = self.DoC[0]
+
+    def params(self):        
         ws = wMOORE(self.D50)
-        gamma = self.Hs / (ws * self.Tp)
+        gamma = self.Hs50 / (ws * self.Tp50)
 
         self.Ar = 0.21 - 0.02 * gamma
         self.B = 0.89 * np.exp(-1.24 * gamma)
@@ -51,8 +60,8 @@ def calibrate(self):
         self.c = self.C**(-1.5)
         self.d = self.D / self.C**(1.5)
 
-        self.ha = 3 * self.Hs
-        self.hr = 1.1 * self.Hs
+        self.ha = 3 * self.Hs50
+        self.hr = 1.1 * self.Hs50
 
         return self       
 
@@ -91,8 +100,9 @@ def Bernabeu(self):
     X = np.concatenate([xrot, xaso])
     hm = np.concatenate([hrot, haso])
 
-    self.hmi = interp1d(X, hm, kind='linear', fill_value='extrapolate')(self.dp)  
-
-    # err = RMSEq(self.zp, self.hmi)
+    self.zm = -interp1d(X, hm, kind='linear', fill_value='extrapolate')(self.xm) + self.HTL 
+    if self.Switch_Cal_DoC == 1:
+        self.xm_DoC = np.mean(self.xm[(self.zm <= self.DoC + 0.05) & (self.zm >= self.DoC - 0.05)])
+        self.zm_DoC = np.mean(self.zm[(self.zm <= self.DoC + 0.05) & (self.zm >= self.DoC - 0.05)])
 
-    return self   
+    return self  

src/IHSetBernabeu/tests/test_bernabeu.py

@@ -3,8 +3,8 @@
 import matplotlib.pyplot as plt
 
 wrkDir = os.getcwd()
-model = IHSetBernabeu.cal_Bernabeu(wrkDir+'/data/')
-self = IHSetBernabeu.Bernabeu(model.calibrate())
+model = IHSetBernabeu.cal_Bernabeu(wrkDir+'/data/prof.csv',wrkDir+'/data/wav.nc', 1, 1, Xm = [0, 750], Hs50 = 0.4466, Tp50 = 10.3337, D50 = 0.5, HTL = 0.0, LTL = -3.1264)
+self = IHSetBernabeu.Bernabeu(model.params())
 
 plt.rcParams.update({'font.family': 'serif'})
 plt.rcParams.update({'font.size': 7})
@@ -13,9 +13,46 @@
         'weight': 'bold',
         'size': 8}
 
-hk = []
-hk.append(plt.plot(self.dp, self.zp, '--k')[0])
-hk.append(plt.plot(self.dp, self.hmi, linewidth=2)[0])
-plt.xlabel('Distance [m]', fontdict=font)
+if self.Switch_Cal_DoC == 0:
+        plt.plot(self.xm, self.zm, '-', color=[0.8, 0.8, 0], linewidth=2, label='Bernabeu profile')[0]
+        plt.fill(self.xm.tolist() + [min(self.xm), min(self.xm)],
+                 self.zm.tolist() + [max(self.zm), min(self.zm)], color='yellow', alpha=0.5)
+        xLim = self.xm[self.zm <= 5]
+        plt.xlim([self.xm[0]-20,xLim[-1]+20])
+        plt.ylim([self.HTL-1,5+0.5])
+        plt.plot([xLim[-1],xLim[-1]-5,xLim[-1]+5,xLim[-1]],
+                 [self.HTL,self.HTL-0.25,self.HTL-0.25,self.HTL], 'b', linewidth=2)
+        plt.text(xLim[-1], self.HTL-0.4,'HTL', fontdict=font, horizontalalignment='center', verticalalignment='center')
+        plt.plot([xLim[-1],xLim[-1]-5,xLim[-1]+5,xLim[-1]],
+                 [self.LTL,self.LTL-0.25,self.LTL-0.25,self.LTL], 'b', linewidth=2)
+        plt.text(xLim[-1], self.LTL-0.4,'LTL', fontdict=font, horizontalalignment='center', verticalalignment='center')
+
+if self.Switch_Cal_DoC == 1:
+        plt.plot(self.xm, self.zm, '-', color=[0.8, 0.8, 0], linewidth=2, label='Bernabeu profile')[0]
+        xm_DoC_fill = self.xm[self.zm <= self.DoC]
+        zm_DoC_fill = self.zm[self.zm <= self.DoC]
+        plt.fill(xm_DoC_fill.tolist() + [min(xm_DoC_fill), min(xm_DoC_fill)],
+                 zm_DoC_fill.tolist() + [max(zm_DoC_fill), min(zm_DoC_fill)], color='yellow', alpha=0.5)
+        plt.plot(self.xm_DoC, self.zm_DoC, 'ro', markersize=8)
+        plt.text(self.xm_DoC, self.zm_DoC-0.2,'h*', fontdict=font, horizontalalignment='center', verticalalignment='center')
+        plt.plot([xm_DoC_fill[-1],xm_DoC_fill[-1]-5,xm_DoC_fill[-1]+5,xm_DoC_fill[-1]],
+                 [self.HTL,self.HTL-0.25,self.HTL-0.25,self.HTL], 'b', linewidth=2)
+        plt.text(xm_DoC_fill[-1], self.HTL-0.4,'HTL', fontdict=font, horizontalalignment='center', verticalalignment='center')
+        plt.plot([xm_DoC_fill[-1],xm_DoC_fill[-1]-5,xm_DoC_fill[-1]+5,xm_DoC_fill[-1]],
+                 [self.LTL,self.LTL-0.25,self.LTL-0.25,self.LTL], 'c', linewidth=2)
+        plt.text(xm_DoC_fill[-1], self.LTL-0.4,'LTL', fontdict=font, horizontalalignment='center', verticalalignment='center')
+        plt.xlim([self.xm[0]-20,xm_DoC_fill[-1]+20])
+        plt.ylim([self.HTL-1,self.DoC+0.5])
+
+if self.Switch_Obs == 1:
+        plt.plot(self.xp, self.zp, '--k', linewidth=2, label='Observed profile')[0]
+plt.plot([min(self.xm), max(self.xm)], [self.HTL, self.HTL], '--b', linewidth=2)[0]
+plt.plot([min(self.xm), max(self.xm)], [self.LTL, self.LTL], '--c', linewidth=2)[0]
+
+plt.xlabel('Offshore distance [m]', fontdict=font)
 plt.ylabel('Water depth [m]', fontdict=font)
+
+plt.grid(True)
+plt.gca().invert_yaxis()
+plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.1), ncol=2)
 plt.show()