Repository: SIO-ODF/ctdcal (language: python)

Issue:
Implement TS plots
Not a param vs param plot - TS plots that include density surfaces are useful for identifying mistrips or highlight SSSCC-related sensor drift in a way different than simple residual plots.

See the 2022 osnap branch for an example.


Discussion / hints:
```
# TODO: more plots! what does ODV have?
# section plots (w/ PyDiva)
# parameter-parameter plots
```

Old TODO marker - add functionality for parameter-parameter to turn off density lines

Fix patch:
diff --git a/ctdcal/ctd_plots.py b/ctdcal/ctd_plots.py
index f924923f..9ce1f01b 100755
--- a/ctdcal/ctd_plots.py
+++ b/ctdcal/ctd_plots.py
@@ -4,6 +4,7 @@
 import logging
 from pathlib import Path
 
+import gsw
 import matplotlib.pyplot as plt
 import matplotlib.ticker as ticker
 import numpy as np
@@ -222,3 +223,71 @@ def _intermediate_residual_plot(
 
     # save to path or return axis (primarily for testing)
     return _save_fig(ax, f_out)
+
+
+def param_vs_param(param1, label1, param2, label2, f_out=None, stn=None, tsT=None, tsS=None):
+    """
+    Scatter plots two CTD parameters against each other with support for TS plotting.
+
+    Support for bottle or time data.
+
+    Most TS plots have temperature on the Y axis and salinity on the X axis.
+    
+    Parameters
+    ----------
+    param1 : array-like of float
+        The first column/array for plotting on the X axis
+    label1 : string
+        Label for the X axis, attributed to param1
+    param2 : array-like of float
+        The second column/array for plotting on the Y-axis
+    label2 : string
+        Label for the Y axis, attributed to param2
+    f_out : String, optional
+        Path and filename to save TS plot
+    stn : array-like, optional
+        An array of station numbers or labels for 
+    tsT : array-like of float, optional
+        Temperature array for use in density gradient calculations
+    tsS : array-like of float, optional
+        Salinity array for use in density gradient calculations
+    """
+
+    fig = plt.figure(figsize=(7, 6))
+    ax = fig.add_subplot()
+
+    if stn is not None:
+        stn = np.array(stn)  # Ensure stn is a NumPy array
+        sc = ax.scatter(param1, param2, c=stn.astype(int), marker="+")
+        sc.set_clim(vmin=stn.astype(int).min(), vmax=stn.astype(int).max()) #   Color mapping based on station number
+        cbar = plt.colorbar(sc, format=ticker.FormatStrFormatter("%.0f"))
+        cbar.set_label("Station/Cast Label")
+    else:
+        sc = ax.scatter(param1, param2, color='black', marker="+")
+
+    ax.set_xlabel(label1)
+    ax.set_ylabel(label2)
+
+    #   TS plotting - define density contours
+    if tsT is not None:
+        #   Demarcate temperature, salinty, and limits for density
+        smin = tsS.min() - (0.01 * tsS.min())
+        smax = tsS.max() + (0.01 * tsS.max())
+        tmin = tsT.min() - (0.1 * tsT.max())
+        tmax = tsT.max() + (0.1 * tsT.max())
+        #   Gridding for a contour
+        xdim = round((smax - smin) / 0.1 + 1.0)
+        ydim = round((tmax - tmin) + 1.0)
+        #   Creating density vectors to fill the grid with densities
+        dens = np.zeros((ydim, xdim))
+        si = np.linspace(1, xdim - 1, xdim) * 0.1 + smin
+        ti = np.linspace(1, ydim - 1, ydim) + tmin
+        for j in range(0, int(ydim)):
+            for i in range(0, int(xdim)):
+                dens[j, i] = gsw.rho(si[i], ti[j], 0)
+        dens = dens - 1000  #   Convert to sigma-theta
+        #   Density contours formatting
+        cs = plt.contour(si, ti, dens, linestyles="dashed", colors="k")
+        plt.clabel(cs, fontsize=12, inline=1, fmt="%0.1f")
+
+    return _save_fig(ax, f_out) #   Save the figure if f_out is given else return the scatter axis


Test patch:
diff --git a/ctdcal/tests/test_ctd_plots.py b/ctdcal/tests/test_ctd_plots.py
index e6b45ad6..a38b9dc4 100644
--- a/ctdcal/tests/test_ctd_plots.py
+++ b/ctdcal/tests/test_ctd_plots.py
@@ -181,3 +181,49 @@ def test_intermediate_residual_plot(tmp_path):
         assert not file_path.exists()
         ctd_plots.residual_vs_pressure(x, x, x, stn=x, f_out=file_path)
         assert file_path.exists()   # Verify that the file has been created
+
+
+@pytest.fixture
+def example_data():
+    return {
+        'param1': np.array([33, 33.5, 34, 35, 35.5]),
+        'param2': np.array([0.1, 10, 15, 20, 20]),
+        'label1': "Salinity",
+        'label2': "Temperature",
+        'stn': np.array(["1", "2", "3", "4", "5"]),
+        'tsS': np.array([33, 33.5, 34, 35, 35.5]),
+        'tsT': np.array([0.1, 10, 15, 20, 20]),
+    }
+
+def test_param_vs_param(tmp_path, example_data):
+    for ext in [".jpg", ".png", ".pdf"]:
+        output_path = tmp_path / "figures"
+        output_path.mkdir(parents=True, exist_ok=True)
+        file_path = output_path / f"not_TS{ext}"
+
+        #   Check that the plotted figure generates and is saved as expected
+        assert not file_path.exists()
+        ax = ctd_plots.param_vs_param(example_data["param1"], example_data["label1"], 
+                                      example_data["param2"], example_data["label2"])   #   Should return figure axis
+        assert np.allclose(ax.get_xticks(), [32.5, 33. , 33.5, 34. , 34.5, 35. , 35.5, 36. ])
+        assert np.allclose(ax.get_yticks(), [-2.5,  0. ,  2.5,  5. ,  7.5, 10. , 12.5, 15. , 17.5, 20. , 22.5])
+        assert ax.get_xlabel() == example_data["label1"]
+        assert ax.get_ylabel() == example_data["label2"]
+        ctd_plots.param_vs_param(example_data["param1"], example_data["label1"], 
+                                 example_data["param2"], example_data["label2"], f_out=file_path)
+        assert file_path.exists()
+
+        #   Now check that the TS contour plotting subroutine works
+        file_path = output_path / f"TS{ext}"
+        assert not file_path.exists()
+        ax = ctd_plots.param_vs_param(example_data["param1"], example_data["label1"], 
+                                      example_data["param2"], example_data["label2"], 
+                                      stn=example_data["stn"], tsT=example_data["tsT"], tsS=example_data["tsS"])
+        assert np.allclose(ax.get_xticks(), [32.5, 33. , 33.5, 34. , 34.5, 35. , 35.5, 36. ])  #   Should be unchanged
+        assert len(ax.collections[-1].get_paths()) == 9 #   Should generate 9 contours
+        assert len(ax.collections[-1].labelTexts) == 7  #   Should have text on 7 of the contour lines
+        assert ax.collections[-1].labelTexts[4].get_text() == "26.4"    #   The density calculations should have produced specific text on the 5th contour
+        ctd_plots.param_vs_param(example_data["param1"], example_data["label1"], 
+                                 example_data["param2"], example_data["label2"], 
+                                 f_out=file_path, stn=example_data["stn"], tsT=example_data["tsT"], tsS=example_data["tsS"])
+        assert file_path.exists()