Dev -> Main

.gitignore

@@ -2,4 +2,5 @@ venv/
 build/
 *.egg-info/
 .vscode/
-*.eggs
+*.eggs
+*__pycache__*

docs/index.rst

@@ -288,6 +288,53 @@ Joint Distribution Hs-Tp Plot
 .. image:: files/Hs.Tp.joint.distribution.png
    :width: 500
 
+Joint 2D distribution contours
+------------------------------
+
+.. code-block:: python
+
+   plots.plot_joint_2D_contour(
+       df, 
+       var1='HS', 
+       var2='TP', 
+       return_periods=[1,10,100,1000], 
+   )
+
+.. image:: files/contours.png
+   :width: 500
+
+Joint 3D distribution contour
+------------------------------
+
+.. code-block:: python
+
+   plots.plot_joint_3D_contour(
+       df, 
+       var1 = 'W10',
+       var2 = 'HS', 
+       var3 = 'TP', 
+       return_period=100
+   )
+
+.. image:: files/3D_contour.png
+   :width: 500
+
+Joint 3D distribution contour cross-sections
+------------------------------
+
+.. code-block:: python
+
+   plots.plot_joint_3D_contour_slices(
+      df,
+      var1 = 'W10',
+      var2 = 'HS',
+      var3 = 'TP',
+      slice_values = [5,10,15,20,25]
+   )
+
+.. image:: files/3D_contour_slices.png
+   :width: 500
+
 Monthly Joint Distribution Hs-Tp Parameter Table
 ------------------------------------------------
 

examples/example_joint_model.py

@@ -0,0 +1,65 @@
+import os
+from metocean_stats.stats.aux_funcs import readNora10File
+from metocean_stats.CMA import JointProbabilityModel,predefined
+import matplotlib.pyplot as plt
+plt.rcParams["axes.grid"] = True
+
+data = readNora10File(os.path.join(os.path.dirname(__file__),'../tests/data/NORA_test.txt'))
+
+# Initialize and fit 2D model of Hs (wave height) and U (wind speed)
+model = JointProbabilityModel(predefined.get_DNVGL_Hs_U)
+model.fit(data,"HS","W10") # fit to HS and W10
+
+# Create the 2D contours
+ax=model.plot_contours(periods=[1,10,100,1000],method="IFORM")
+model.plot_data_scatter(ax,s=2)
+model.plot_dependent_percentiles(ax)
+model.plot_legend(ax)
+model.reset_labels()
+
+# Fitted distributions and histograms per interval
+plt.rcParams["figure.figsize"] = (20,20)
+fig,ax=model.plot_histograms_of_interval_distributions()
+fig[1].subplots_adjust(hspace=0.3)
+
+# Dependency functions of U on Hs
+plt.rcParams["figure.figsize"] = (15,7)
+fig,ax = plt.subplots(1,2)
+model.plot_dependence_functions(ax)
+
+# Analytical vs empirical quantiles
+fig,ax = plt.subplots(1,2)
+model.plot_marginal_quantiles(ax)
+
+# Initialize 3D model of wind, wave height and wave period
+model = JointProbabilityModel(predefined.get_LiGaoMoan_U_hs_tp)
+model.fit(data,"W10","HS","TP")
+
+# Full 3D contour obtained with IFORM
+model.plot_3D_contour(return_period=100,n_samples=720)
+
+# 2D slices of the 3D contour at a range of wind levels
+plt.rcParams["figure.figsize"] = (20,12)
+ax=model.plot_3D_contour_slices(subplots=False,slice_values=[i for i in range(1,40)])
+model.plot_legend(ax)
+model.reset_labels()
+
+# Isodensity contour - defined by a marginal wind return value
+model.plot_3D_isodensity_contour(marginal_value=30)
+
+# 2D slices of isodensity contours, at a range of probability densities.
+plt.rcParams["figure.figsize"] = (20,7)
+model.reset_labels()
+ax=model.plot_3D_isodensity_contour_slice(
+    density_levels=[1e-10,1e-9,1e-8,1e-7,1e-6,1e-5,1e-4,1e-3,2e-3,4e-3,6e-3],slice_value=10)
+ax.set_xlim([0,40])
+ax.set_ylim([0,10])
+model.plot_legend(ax)
+model.reset_labels()
+
+# Get a dataframe describing the joint model.
+params = model.parameters()
+print(params)
+
+# Show figures
+plt.show()

metocean_stats/CMA/contours.py

@@ -41,7 +41,8 @@
     "IFORMContour",
     "ISORMContour",
     "get_contour",
-    "sort_contour"
+    "sort_contour",
+    "split_contour",
 ]
 
 def get_contour(
@@ -154,6 +155,55 @@ def sort_contour(contour):
         if nearest_idx == 0:
             return contour[ordered_indices]
 
+def split_contour(contour,split_dim):
+    """
+    Where "contour" is a Nx2 array, 
+    where each point in the array is placed to the previous, 
+    divide into left hand side and right hand side on the second coordinate (e.g., Tp),
+    and sort from low to high on the first coordinate (e.g., Hs).
+    """
+    range_dim = 0 if split_dim else 1
+
+    bot_idx = np.argmin(contour[:,range_dim])
+    top_idx = np.argmax(contour[:,range_dim])
+    i = 0
+    state = "none"
+    lhs = []
+    rhs = []
+    while True:
+        if i == len(contour): i = 0
+
+        # We start counting at the bottom of the contour
+        if state == "none" and i == bot_idx:
+            if contour[i-1,split_dim] < contour[i,split_dim]:
+                state = "rhs"
+            else:
+                state = "lhs"
+
+        # If we reached top, switch side
+        elif i == top_idx:
+            if state == "lhs": 
+                state = "rhs"
+            elif state == "rhs":
+                state = "lhs"
+
+        # If we reached bottom again, we stop
+        elif i == bot_idx: break
+
+        # Adding values depending on state
+        if state == "lhs": lhs.append(contour[i])
+        if state == "rhs": rhs.append(contour[i])
+
+        # Increment
+        i += 1
+
+    lhs = np.array(lhs)
+    lhs = lhs[np.argsort(lhs[:,range_dim])]
+    rhs = np.array(rhs)
+    rhs = rhs[np.argsort(rhs[:,range_dim])]
+
+    return lhs,rhs
+
 def get_random_unit_sphere_points(n_dim,n_samples):
     """
     Generates equally distributed points on the sphere's surface.
@@ -297,6 +347,7 @@ def _compute(
             unit_sphere_points = get_random_unit_sphere_points(n_dim=self.model.n_dim,n_samples=self.n_points)
         elif self.point_distribution == "gridded":
             unit_sphere_points = get_3D_gridded_unit_sphere_points(primary_dim=0,n_polar=self.n_points,n_azimuthal=self.n_points)
+
         sphere_points = beta * unit_sphere_points
 
         # Get probabilities for coordinates