Add files via upload

pyproject.toml

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

src/IHSetDean/IHSetDean.py

@@ -1,7 +1,8 @@
 import numpy as np
 from scipy.interpolate import interp1d
 import xarray as xr
-from IHSetUtils import wMOORE, ADEAN
+import pandas as pd
+from IHSetUtils import wMOORE, ADEAN, Hs12Calc, depthOfClosure
 
 class cal_Dean(object):
     """
@@ -11,37 +12,61 @@ class cal_Dean(object):
     
     This class reads input datasets, performs its calibration.
     """
-    def __init__(self, path):
-        self.path = path
-
-        cfg = xr.open_dataset(path+'config.nc')
-        ens = xr.open_dataset(path+'ens.nc')
-
-        self.Dmin = cfg['Dmin'].values
-        self.Dmax = cfg['Dmax'].values
-        self.dCal = cfg['dCal'].values
-        self.D50 = ens['D50'].values
-        self.dp = ens['d'].values
-        self.zp = ens['z'].values
-
+    def __init__(self, path_prof, path_wav, Switch_Calibrate, Switch_Cal_DoC, **kwargs):
+        self.path_prof = path_prof
+        self.path_wav = path_wav
+        prof = pd.read_csv(path_prof)
+        self.D50 = kwargs['D50']                  # D50 is not necesarry if Switch_Calibrate = 1
+        self.MSL = -kwargs['MSL']       
+        self.xm = np.linspace(kwargs['Xm'][0], kwargs['Xm'][1], 1000).reshape(-1, 1)
+
+        self.Switch_Calibrate = Switch_Calibrate
+        # Calibrate Dean parameter using profile data [0: no without obs (using D50); 1: no with obs (using D50); 2: yes (using Obs)]
+        if Switch_Calibrate == 1 or Switch_Calibrate == 2:
+            self.xp = prof.iloc[:, 0]
+            self.zp = prof.iloc[:, 1]
+            self.zp = abs(self.zp)
+            xp_inx = self.xp[(self.zp >= self.MSL)]
+            self.xp = self.xp - min(xp_inx)
+
+        if Switch_Calibrate == 2:
+            self.Zmin = kwargs['Zmin']
+            self.Zmax = kwargs['Zmax']
+
+        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 calibrate(self):
-        self.zp = self.zp - self.zp[0]
-        self.d = self.dp - self.dp[0]
-
-        # Profile with equidistant points
-        dp = np.linspace(self.Dmin, self.Dmax, 500).reshape(-1, 1)
-        interp_func = interp1d(self.d, self.zp, kind="linear", fill_value="extrapolate")
-        zp = interp_func(dp)
-        zp = zp[1:]
-        dp = dp[1:]
 
-        if self.dCal == 1:    
+        if self.Switch_Calibrate == 2:
+            self.xpp = self.xp[(self.zp >= self.MSL)]
+            self.zpp = self.zp[(self.zp >= self.MSL)]
+
+            if self.Zmin < np.min(self.zpp):
+                self.Zmin = np.min(self.zpp)
+            if self.Zmax > np.max(self.zpp):
+                self.Zmax = np.max(self.zpp)
+            zp = np.linspace(self.Zmin, self.Zmax, 100).reshape(-1, 1)
+            interp_func = interp1d(self.zpp, self.xpp, kind="linear", fill_value="extrapolate")
+            xp = interp_func(zp)
+            zp = -zp[1:] + self.MSL
+            xp = xp[1:]
+
             ws = None
             if self.D50 is not None:
                 ws = wMOORE(self.D50)
 
             Y = np.log(-zp)
-            Y2 = 2 / 3 * np.log(dp)
+            Y2 = 2 / 3 * np.log(xp)
 
             fc = np.arange(-20, 20, 0.0001)
             Y2_grid, fc_grid = np.meshgrid(Y2, fc, indexing="ij")
@@ -55,16 +80,19 @@ def RMSEq(Y, Y2t):
 
             self.A = np.exp(fc[I])
             self.kk = np.exp(fc[I] - np.log(ws**0.44)) if ws is not None else None
-            self.dp_value = dp
+
+            # self.zp = self.zp - self.MSL
 
-        if self.dCal == 0:
+        if self.Switch_Calibrate == 0 or self.Switch_Calibrate == 1:
             D50 = self.D50 / 1000
             self.A = ADEAN(D50)
 
         return self
 
 def Dean(self):
-
-    self.hm = -self.A * self.dp ** (2 / 3)
+    self.zm = self.A * self.xm ** (2 / 3) + self.MSL
+    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

src/IHSetDean/tests/test_dean.py

@@ -1,16 +1,9 @@
 from IHSetDean import IHSetDean
-import xarray as xr
 import os
 import matplotlib.pyplot as plt
 
-config = xr.Dataset(coords={'dCal': 0,                # Calibrate using the wave data : 1 / Not calibrate : 0 (using D50)
-                            'Dmin': 0,                # Calibrate the minimum observed distnace value
-                            'Dmax': 500,               # Calibrate the maximum observed distance value
-                            })
-
 wrkDir = os.getcwd()
-config.to_netcdf(wrkDir+'/data/config.nc', engine='netcdf4')
-model = IHSetDean.cal_Dean(wrkDir+'/data/')
+model = IHSetDean.cal_Dean(wrkDir+'/data/prof.csv',wrkDir+'/data/wav.nc', 2, 1, Xm = [0, 500], Zmin = 0.0, Zmax = 1.0, D50 = 0.5, MSL = 0.0)
 self = IHSetDean.Dean(model.calibrate())
 
 plt.rcParams.update({'font.family': 'serif'})
@@ -20,9 +13,42 @@
         'weight': 'bold',
         'size': 8}
 
-hk = []
-hk.append(plt.plot(self.dp, self.zp, '--k')[0])
-hk.append(plt.plot(self.dp, self.hm, 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='Dean 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.MSL-1,5+0.5])
+        plt.plot([xLim[-1],xLim[-1]-5,xLim[-1]+5,xLim[-1]],
+                 [self.MSL,self.MSL-0.25,self.MSL-0.25,self.MSL], 'b', linewidth=2)
+        plt.text(xLim[-1], self.MSL-0.4,'MSL', 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='Dean 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.MSL,self.MSL-0.25,self.MSL-0.25,self.MSL], 'b', linewidth=2)
+        plt.text(xm_DoC_fill[-1], self.MSL-0.4,'MSL', fontdict=font, horizontalalignment='center', verticalalignment='center')
+        plt.xlim([self.xm[0]-20,xm_DoC_fill[-1]+20])
+        plt.ylim([self.MSL-1,self.DoC+0.5])
+
+if self.Switch_Calibrate == 2:
+        plt.fill([min(self.xm), min(self.xm), max(self.xm), max(self.xm)],
+                 [self.Zmin, self.Zmax, self.Zmax, self.Zmin], color=[0.5, 0.5, 0.5], alpha=0.25)
+if self.Switch_Calibrate == 1 or self.Switch_Calibrate == 2:
+        plt.plot(self.xp, self.zp, '--k', linewidth=2, label='Observed profile')[0]
+plt.plot([min(self.xm), max(self.xm)], [self.MSL, self.MSL], '--b', linewidth=2)[0]
+
+plt.xlabel('Offshore distance [m]', fontdict=font)
 plt.ylabel('Water depth [m]', fontdict=font)
-plt.show()
+
+plt.grid(True)
+plt.gca().invert_yaxis()
+plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.1), ncol=2)
+plt.show()