From cc149c7e34d683caa5e44708e7f3d2ec0943765d Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Thu, 11 Dec 2025 12:39:20 -0800
Subject: [PATCH 01/20] Update instance attributes from parsed arguments.

---
 src/data/create_products.py | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index fdf0806c..4c869f3e 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -572,9 +572,18 @@ def process_command_line(self):
         )
 
         self.args = parser.parse_args()
-        self.logger.setLevel(self._log_levels[self.verbose])
         self.commandline = " ".join(sys.argv)
 
+        # Update instance attributes from parsed arguments
+        self.auv_name = self.args.auv_name
+        self.mission = self.args.mission
+        self.base_path = self.args.base_path
+        self.start_esecs = self.args.start_esecs
+        self.local = self.args.local
+        self.verbose = self.args.verbose
+
+        self.logger.setLevel(self._log_levels[self.args.verbose])
+
 
 if __name__ == "__main__":
     cp = CreateProducts()

From 0593333d997b741c4f44a81428b6ed9eec6bc8be Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Thu, 11 Dec 2025 14:47:16 -0800
Subject: [PATCH 02/20] WIP: Generate classic Dorado mission Quick Look plots.

---
 src/data/create_products.py | 246 +++++++++++++++++++++++++++++++-----
 1 file changed, 212 insertions(+), 34 deletions(-)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index 4c869f3e..ef16d81c 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -17,6 +17,7 @@
 from pathlib import Path
 
 import cmocean
+import gsw
 import matplotlib  # noqa: ICN001
 import matplotlib.pyplot as plt
 import numpy as np
@@ -109,6 +110,7 @@ def __init__(  # noqa: PLR0913
         "sea_water_temperature": "thermal",
         "sea_water_salinity": "haline",
         "sea_water_sigma_t": "dense",
+        "sea_water_density": "dense",
         "mass_concentration_of_chlorophyll_in_sea_water": "algae",
         "mass_concentration_of_oxygen_in_sea_water": "oxy",
         "downwelling_photosynthetic_photon_flux_in_sea_water": "solar",
@@ -143,6 +145,49 @@ def _open_ds(self):
             self.logger.debug("%s not available yet", dap_url)
             self.ds = xr.open_dataset(local_nc)
 
+    def _compute_density(self, best_ctd: str = "ctd1") -> None:
+        """Compute sigma-t density from temperature and salinity using EOS-80.
+
+        Args:
+            best_ctd: The CTD instrument to use for temperature and salinity
+        """
+        if "density" in self.ds:
+            self.logger.debug("density already exists in dataset")
+            return
+
+        temp_var = f"{best_ctd}_temperature"
+        sal_var = f"{best_ctd}_salinity"
+
+        if temp_var not in self.ds or sal_var not in self.ds:
+            self.logger.warning(
+                "Cannot compute density: %s or %s not in dataset",
+                temp_var,
+                sal_var,
+            )
+            return
+
+        # Get temperature, salinity, and pressure (approximated from depth)
+        temp = self.ds[temp_var].to_numpy()
+        sal = self.ds[sal_var].to_numpy()
+
+        # Compute sigma-t using gsw (uses EOS-80 formulation)
+        # sigma-t is density - 1000 kg/m³
+        density = gsw.density.sigma0(sal, temp)
+
+        # Add to dataset
+        self.ds["density"] = xr.DataArray(
+            density,
+            dims=self.ds[temp_var].dims,
+            coords=self.ds[temp_var].coords,
+            attrs={
+                "long_name": "Sigma-t (density anomaly)",
+                "standard_name": "sea_water_density",
+                "units": "kg/m^3",
+                "comment": f"Computed from {temp_var} and {sal_var} using gsw.density.sigma0",
+            },
+        )
+        self.logger.info("Computed density (sigma-t) from %s and %s", temp_var, sal_var)
+
     def _grid_dims(self) -> tuple:
         # From Matlab code in plot_sections.m:
         # auvnav positions are too fine for distance calculations, they resolve
@@ -242,9 +287,28 @@ def _plot_var(  # noqa: PLR0913
         scale: str = "linear",
         num_colors: int = 256,
     ):
+        # Handle both 1D and 2D axis arrays
+        curr_ax = ax[row] if col == 0 and hasattr(ax, "ndim") and ax.ndim == 1 else ax[row, col]
+
         var_to_plot = (
             np.log10(self.ds[var].to_numpy()) if scale == "log" else self.ds[var].to_numpy()
         )
+
+        # Check if variable has any valid (non-NaN) data
+        valid_data = var_to_plot[~np.isnan(var_to_plot)]
+        if len(valid_data) == 0:
+            self.logger.warning("%s has no valid data, skipping plot", var)
+            curr_ax.text(
+                0.5,
+                0.5,
+                f"No valid data for {var}",
+                ha="center",
+                va="center",
+                transform=curr_ax.transAxes,
+            )
+            curr_ax.set_ylabel("Depth (m)")
+            return
+
         scafac = max(idist) / max(iz)
         gridded_var = griddata(
             (distnav.to_numpy() / 1000.0 / scafac, self.ds.cf["depth"].to_numpy()),
@@ -265,8 +329,28 @@ def _plot_var(  # noqa: PLR0913
             # Likely a cmocean colormap
             cmap = getattr(cmocean.cm, color_map_name)
 
-        v2_5 = np.percentile(var_to_plot[~np.isnan(var_to_plot)], 2.5)
-        v97_5 = np.percentile(var_to_plot[~np.isnan(var_to_plot)], 97.5)
+        v2_5 = np.percentile(valid_data, 2.5)
+        v97_5 = np.percentile(valid_data, 97.5)
+
+        # Additional check: ensure percentiles are valid numbers
+        if np.isnan(v2_5) or np.isnan(v97_5) or v2_5 == v97_5:
+            self.logger.warning(
+                "%s has invalid range (v2.5=%.2f, v97.5=%.2f), skipping plot",
+                var,
+                v2_5,
+                v97_5,
+            )
+            curr_ax.text(
+                0.5,
+                0.5,
+                f"Invalid data range for {var}",
+                ha="center",
+                va="center",
+                transform=curr_ax.transAxes,
+            )
+            curr_ax.set_ylabel("Depth (m)")
+            return
+
         norm = matplotlib.colors.BoundaryNorm(
             np.linspace(v2_5, v97_5, num_colors),
             num_colors,
@@ -279,8 +363,8 @@ def _plot_var(  # noqa: PLR0913
             v2_5,
             v97_5,
         )
-        ax[row, col].set_ylim(max(iz), min(iz))
-        cntrf = ax[row, col].contourf(
+        curr_ax.set_ylim(max(iz), min(iz))
+        cntrf = curr_ax.contourf(
             idist / 1000.0,
             iz,
             gridded_var,
@@ -293,39 +377,42 @@ def _plot_var(  # noqa: PLR0913
         xb = np.append(
             profile_bottoms.get_index("dist").to_numpy(),
             [
-                ax[row, col].get_xlim()[1],
-                ax[row, col].get_xlim()[1],
-                ax[row, col].get_xlim()[0],
-                ax[row, col].get_xlim()[0],
+                curr_ax.get_xlim()[1],
+                curr_ax.get_xlim()[1],
+                curr_ax.get_xlim()[0],
+                curr_ax.get_xlim()[0],
             ],
         )
         yb = np.append(
             profile_bottoms.to_numpy(),
             [
                 profile_bottoms.to_numpy()[-1],
-                ax[row][col].get_ylim()[0],
-                ax[row, col].get_ylim()[0],
+                curr_ax.get_ylim()[0],
+                curr_ax.get_ylim()[0],
                 profile_bottoms.to_numpy()[0],
             ],
         )
-        ax[row][col].fill(list(reversed(xb)), list(reversed(yb)), "w")
+        curr_ax.fill(list(reversed(xb)), list(reversed(yb)), "w")
 
-        ax[row, col].set_ylabel("Depth (m)")
-        cb = fig.colorbar(cntrf, ax=ax[row, col])
+        curr_ax.set_ylabel("Depth (m)")
+        cb = fig.colorbar(cntrf, ax=curr_ax)
         cb.locator = matplotlib.ticker.LinearLocator(numticks=3)
         cb.minorticks_off()
         cb.update_ticks()
         cb.ax.set_yticklabels([f"{x:.1f}" for x in cb.get_ticks()])
-        if scale == "log":
-            cb.set_label(
-                f"{self.ds[var].attrs['long_name']}\n[log10({self.ds[var].attrs['units']})]",
-                fontsize=7,
-            )
+
+        # Get long_name and units with fallbacks
+        long_name = self.ds[var].attrs.get("long_name", var)
+        units = self.ds[var].attrs.get("units", "")
+
+        if scale == "log" and units:
+            cb.set_label(f"{long_name}\n[log10({units})]", fontsize=7)
+        elif scale == "log":
+            cb.set_label(f"{long_name}\n[log10]", fontsize=7)
+        elif units:
+            cb.set_label(f"{long_name} [{units}]", fontsize=9)
         else:
-            cb.set_label(
-                f"{self.ds[var].attrs['long_name']} [{self.ds[var].attrs['units']}]",
-                fontsize=9,
-            )
+            cb.set_label(long_name, fontsize=9)
 
     def plot_2column(self) -> str:
         """Create 2column plot similar to plot_sections.m and stoqs/utils/Viz/plotting.py
@@ -338,21 +425,25 @@ def plot_2column(self) -> str:
         scfac = max(idist) / max(iz)  # noqa: F841
 
         fig, ax = plt.subplots(nrows=5, ncols=2, figsize=(18, 10))
-        fig.tight_layout()
+        fig.tight_layout(rect=[0, 0.03, 1, 0.96])
 
         best_ctd = self._get_best_ctd()
+
+        # Compute density (sigma-t) if not already present
+        self._compute_density(best_ctd)
+
         profile_bottoms = self._profile_bottoms(distnav)
         row = 0
         col = 1
         for var, scale in (
-            ("ctd1_oxygen_mll", "linear"),
             ("density", "linear"),
-            ("hs2_bb420", "linear"),
             (f"{best_ctd}_temperature", "linear"),
-            ("hs2_bb700", "linear"),
             (f"{best_ctd}_salinity", "linear"),
-            ("hs2_fl700", "linear"),
             ("isus_nitrate", "linear"),
+            ("ctd1_oxygen_mll", "linear"),
+            ("hs2_bbp420", "linear"),
+            ("hs2_bb700", "linear"),
+            ("hs2_fl700", "linear"),
             ("biolume_avg_biolume", "log"),
         ):
             self.logger.info("Plotting %s...", var)
@@ -375,6 +466,9 @@ def plot_2column(self) -> str:
                 )
             if row != 4:  # noqa: PLR2004
                 ax[row, col].get_xaxis().set_visible(False)
+            else:
+                # Add x-axis label only to bottom row
+                ax[row, col].set_xlabel("Distance (km)")
 
             if col == 1:
                 row += 1
@@ -382,19 +476,103 @@ def plot_2column(self) -> str:
             else:
                 col = 1
 
+        # Add title to the figure
+        title = f"{self.auv_name} {self.mission}"
+        if "title" in self.ds.attrs:
+            title = self.ds.attrs["title"]
+        fig.suptitle(title, fontsize=12, fontweight="bold")
+
         # Save plot to file
-        images_dir = Path(BASE_PATH, self.auv_name, MISSIONIMAGES)
+        images_dir = Path(BASE_PATH, self.auv_name, MISSIONIMAGES, self.mission)
         Path(images_dir).mkdir(parents=True, exist_ok=True)
 
-        plt.savefig(
-            Path(
-                images_dir,
-                f"{self.auv_name}_{self.mission}_{FREQ}_2column.png",
-            ),
+        output_file = Path(
+            images_dir,
+            f"{self.auv_name}_{self.mission}_{FREQ}_2column.png",
         )
+        plt.savefig(output_file, dpi=100, bbox_inches="tight")
+        plt.show()
+        plt.close(fig)
+
+        self.logger.info("Saved 2column plot to %s", output_file)
+        return str(output_file)
 
     def plot_biolume(self) -> str:
-        "Create biolume plot"
+        """Create bioluminescence plot showing raw signal and proxy variables"""
+        self._open_ds()
+
+        # Check if biolume variables exist
+        biolume_vars = [v for v in self.ds.variables if v.startswith("biolume_")]
+        if not biolume_vars:
+            self.logger.warning("No biolume variables found in dataset")
+            return None
+
+        idist, iz, distnav = self._grid_dims()
+        profile_bottoms = self._profile_bottoms(distnav)
+
+        # Create figure with subplots for biolume variables
+        num_plots = min(len(biolume_vars), 6)  # Limit to 6 most important variables
+        fig, ax = plt.subplots(nrows=num_plots, ncols=1, figsize=(18, 12))
+        if num_plots == 1:
+            ax = [ax]
+        fig.tight_layout(rect=[0, 0.03, 1, 0.96])
+
+        # Priority order for biolume variables to plot
+        priority_vars = [
+            "biolume_avg_biolume",
+            "biolume_bg_biolume",
+            "biolume_nbflash_high",
+            "biolume_nbflash_low",
+            "biolume_proxy_diatoms",
+            "biolume_proxy_adinos",
+        ]
+
+        vars_to_plot = []
+        for pvar in priority_vars:
+            if pvar in self.ds:
+                scale = "log" if "avg_biolume" in pvar or "bg_biolume" in pvar else "linear"
+                vars_to_plot.append((pvar, scale))
+            if len(vars_to_plot) >= num_plots:
+                break
+
+        for i, (var, scale) in enumerate(vars_to_plot):
+            self.logger.info("Plotting %s...", var)
+            self._plot_var(
+                var,
+                idist,
+                iz,
+                distnav,
+                fig,
+                ax,
+                i,
+                0,
+                profile_bottoms,
+                scale=scale,
+            )
+            if i != num_plots - 1:
+                ax[i].get_xaxis().set_visible(False)
+            else:
+                ax[i].set_xlabel("Distance (km)")
+
+        # Add title to the figure
+        title = f"{self.auv_name} {self.mission} - Bioluminescence"
+        if "title" in self.ds.attrs:
+            title = f"{self.ds.attrs['title']} - Bioluminescence"
+        fig.suptitle(title, fontsize=12, fontweight="bold")
+
+        # Save plot to file
+        images_dir = Path(BASE_PATH, self.auv_name, MISSIONIMAGES, self.mission)
+        Path(images_dir).mkdir(parents=True, exist_ok=True)
+
+        output_file = Path(
+            images_dir,
+            f"{self.auv_name}_{self.mission}_{FREQ}_biolume.png",
+        )
+        plt.savefig(output_file, dpi=100, bbox_inches="tight")
+        plt.close(fig)
+
+        self.logger.info("Saved biolume plot to %s", output_file)
+        return str(output_file)
 
     def _get_best_ctd(self) -> str:
         """Determine best CTD to use for ODV lookup table based on metadata"""

From 3de60f488dd37ce9110876e79bf365fa86630885 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Sun, 14 Dec 2025 14:36:46 -0800
Subject: [PATCH 03/20] WIP: Add map plot, deal with no_data, adjust witespace,
 ...

---
 src/data/create_products.py | 251 +++++++++++++++++++++++++++++-------
 1 file changed, 202 insertions(+), 49 deletions(-)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index ef16d81c..f51143ff 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -21,6 +21,7 @@
 import matplotlib  # noqa: ICN001
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 import pyproj
 import xarray as xr
 
@@ -180,7 +181,7 @@ def _compute_density(self, best_ctd: str = "ctd1") -> None:
             dims=self.ds[temp_var].dims,
             coords=self.ds[temp_var].coords,
             attrs={
-                "long_name": "Sigma-t (density anomaly)",
+                "long_name": "Sigma-t",
                 "standard_name": "sea_water_density",
                 "units": "kg/m^3",
                 "comment": f"Computed from {temp_var} and {sal_var} using gsw.density.sigma0",
@@ -245,8 +246,9 @@ def _grid_dims(self) -> tuple:
             distnav.to_numpy().max(),
             3 * self.ds["profile_number"].to_numpy()[-1],
         )
-        # Vertical gridded to .5 m
-        iz = np.arange(2.0, self.ds.cf["depth"].max(), 0.5)
+        # Vertical gridded to .5 m, rounded down to nearest 50m
+        max_depth = np.floor(self.ds.cf["depth"].max() / 50) * 50
+        iz = np.arange(2.0, max_depth, 0.5)
         if not iz.any():
             self.logger.warning(
                 "Gridding vertical for a surface only mission: {self.ds.cf['depth'].max() =}",
@@ -273,7 +275,104 @@ def _profile_bottoms(
         window = int(len(distnav) * window_frac)
         return depth_dist.rolling(dist=window).max()
 
-    def _plot_var(  # noqa: PLR0913
+    def _plot_track_map(
+        self, map_ax: matplotlib.axes.Axes, reference_ax: matplotlib.axes.Axes
+    ) -> None:
+        """Plot AUV track map on left side with title and times on right.
+
+        Args:
+            map_ax: The axes object to plot the map on
+            reference_ax: The axes below to align with (for left edge)
+        """
+        # Get lat/lon data
+        lons = self.ds.cf["longitude"].to_numpy()
+        lats = self.ds.cf["latitude"].to_numpy()
+
+        # Get time data for start/end
+        times = self.ds.cf["time"].to_numpy()
+        start_time = pd.to_datetime(times[0]).strftime("%Y-%m-%d %H:%M:%S")
+        end_time = pd.to_datetime(times[-1]).strftime("%Y-%m-%d %H:%M:%S")
+
+        # Get title from netCDF attributes
+        title = self.ds.attrs.get("title", f"{self.auv_name} {self.mission}")
+
+        # Get the position of the reference axes below to align with
+        ref_pos = reference_ax.get_position()
+
+        # Store original position
+        pos = map_ax.get_position()
+
+        # Make the plot square by using equal aspect (do this first)
+        map_ax.set_aspect("equal", adjustable="datalim")
+
+        # Plot the track with depth coloring
+        depths = self.ds.cf["depth"].to_numpy()
+        map_ax.scatter(
+            lons,
+            lats,
+            c=depths,
+            cmap="viridis_r",  # Reversed so shallow is yellow, deep is dark
+            s=1,
+            alpha=0.6,
+        )
+
+        # Add start and end markers
+        map_ax.plot(lons[0], lats[0], "go", markersize=8, label="Start", zorder=5)
+        map_ax.plot(lons[-1], lats[-1], "r^", markersize=8, label="End", zorder=5)
+
+        # Set fixed axis limits for Monterey Bay area
+        map_ax.set_xlim([-122.41, -121.77])
+        map_ax.set_ylim([36.5, 37.0])
+
+        # Now position map aligned with left edge of reference, 50% width
+        # Use a square aspect ratio based on the y-dimension
+        map_height = pos.height
+        aspect_ratio = (37.0 - 36.5) / (122.41 - 121.77)  # data aspect ratio
+        map_width = map_height / aspect_ratio * 0.7  # scale to fit nicely
+
+        map_ax.set_position([ref_pos.x0, pos.y0, map_width, map_height])
+
+        # Remove axes, labels, and ticks but keep the border
+        map_ax.set_xticks([])
+        map_ax.set_yticks([])
+        map_ax.set_xlabel("")
+        map_ax.set_ylabel("")
+
+        # Add a border around the map
+        for spine in map_ax.spines.values():
+            spine.set_edgecolor("black")
+            spine.set_linewidth(1)
+
+        # Add legend
+        map_ax.legend(loc="upper right", fontsize=7, framealpha=0.9)
+
+        # Add text on the right side with title and times
+        # Wrap title for better formatting
+        import textwrap
+
+        wrapped_title = textwrap.fill(title, width=40)
+
+        text_content = f"{wrapped_title}\n\nStart: {start_time}\nEnd:   {end_time}"
+
+        # Get updated position after aspect adjustment
+        updated_pos = map_ax.get_position()
+
+        # Position text in figure coordinates, to the right of the map
+        text_x = updated_pos.x0 + updated_pos.width + 0.01
+        text_y = updated_pos.y0 + updated_pos.height * 0.5
+
+        map_ax.figure.text(
+            text_x,
+            text_y,
+            text_content,
+            fontsize=11,
+            fontweight="bold",
+            family="sans-serif",
+            verticalalignment="center",
+            horizontalalignment="left",
+        )
+
+    def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
         self,
         var: str,
         idist: np.array,
@@ -286,29 +385,55 @@ def _plot_var(  # noqa: PLR0913
         profile_bottoms: xr.DataArray,
         scale: str = "linear",
         num_colors: int = 256,
+        best_ctd: str = "ctd1",
     ):
         # Handle both 1D and 2D axis arrays
         curr_ax = ax[row] if col == 0 and hasattr(ax, "ndim") and ax.ndim == 1 else ax[row, col]
 
-        var_to_plot = (
-            np.log10(self.ds[var].to_numpy()) if scale == "log" else self.ds[var].to_numpy()
-        )
+        # Check if variable exists and has valid data
+        no_data = False
+        if var not in self.ds:
+            self.logger.warning("%s not in dataset", var)
+            no_data = True
+        else:
+            var_to_plot = (
+                np.log10(self.ds[var].to_numpy()) if scale == "log" else self.ds[var].to_numpy()
+            )
+            valid_data = var_to_plot[~np.isnan(var_to_plot)]
+            if len(valid_data) == 0:
+                self.logger.warning("%s has no valid data", var)
+                no_data = True
+
+        # If no data, set up minimal axes and return early
+        if no_data:
+            curr_ax.set_xlim([min(idist) / 1000.0, max(idist) / 1000.0])
+            curr_ax.set_ylim([max(iz), min(iz)])
+
+            # Only show y-label on left column or top of right column
+            if col == 0 or (col == 1 and row == 0):
+                curr_ax.set_ylabel("Depth (m)")
+            else:
+                curr_ax.set_ylabel("")
 
-        # Check if variable has any valid (non-NaN) data
-        valid_data = var_to_plot[~np.isnan(var_to_plot)]
-        if len(valid_data) == 0:
-            self.logger.warning("%s has no valid data, skipping plot", var)
+            # Set y-axis ticks at 0, 50, 100, 150, etc.
+            y_min = max(iz)
+            y_ticks = np.arange(0, int(y_min) + 50, 50)
+            curr_ax.set_yticks(y_ticks)
+
+            # Add "No Data" text
             curr_ax.text(
                 0.5,
                 0.5,
-                f"No valid data for {var}",
+                "No Data",
                 ha="center",
                 va="center",
+                fontsize=14,
+                fontweight="bold",
                 transform=curr_ax.transAxes,
             )
-            curr_ax.set_ylabel("Depth (m)")
             return
 
+        # Normal plotting path - we have valid data
         scafac = max(idist) / max(iz)
         gridded_var = griddata(
             (distnav.to_numpy() / 1000.0 / scafac, self.ds.cf["depth"].to_numpy()),
@@ -317,6 +442,7 @@ def _plot_var(  # noqa: PLR0913
             method="linear",
             rescale=True,
         )
+
         color_map_name = "cividis"
         with contextlib.suppress(KeyError):
             color_map_name = self.cmocean_lookup.get(
@@ -332,23 +458,40 @@ def _plot_var(  # noqa: PLR0913
         v2_5 = np.percentile(valid_data, 2.5)
         v97_5 = np.percentile(valid_data, 97.5)
 
-        # Additional check: ensure percentiles are valid numbers
+        # Check for invalid percentiles
         if np.isnan(v2_5) or np.isnan(v97_5) or v2_5 == v97_5:
             self.logger.warning(
-                "%s has invalid range (v2.5=%.2f, v97.5=%.2f), skipping plot",
+                "%s has invalid range (v2.5=%.2f, v97.5=%.2f)",
                 var,
                 v2_5,
                 v97_5,
             )
+            # Set up minimal axes and return early
+            curr_ax.set_xlim([min(idist) / 1000.0, max(idist) / 1000.0])
+            curr_ax.set_ylim([max(iz), min(iz)])
+
+            # Only show y-label on left column or top of right column
+            if col == 0 or (col == 1 and row == 0):
+                curr_ax.set_ylabel("Depth (m)")
+            else:
+                curr_ax.set_ylabel("")
+
+            # Set y-axis ticks at 0, 50, 100, 150, etc.
+            y_min = max(iz)
+            y_ticks = np.arange(0, int(y_min) + 50, 50)
+            curr_ax.set_yticks(y_ticks)
+
+            # Add "No Data" text
             curr_ax.text(
                 0.5,
                 0.5,
-                f"Invalid data range for {var}",
+                "No Data",
                 ha="center",
                 va="center",
+                fontsize=14,
+                fontweight="bold",
                 transform=curr_ax.transAxes,
             )
-            curr_ax.set_ylabel("Depth (m)")
             return
 
         norm = matplotlib.colors.BoundaryNorm(
@@ -394,7 +537,18 @@ def _plot_var(  # noqa: PLR0913
         )
         curr_ax.fill(list(reversed(xb)), list(reversed(yb)), "w")
 
-        curr_ax.set_ylabel("Depth (m)")
+        # Only show y-label on left column or top of right column
+        if col == 0 or (col == 1 and row == 0):
+            curr_ax.set_ylabel("Depth (m)")
+        else:
+            curr_ax.set_ylabel("")
+
+        # Set y-axis ticks at 0, 50, 100, 150, etc.
+        y_min, y_max = curr_ax.get_ylim()
+        # Since y-axis is inverted (max at bottom), y_min is the deeper value
+        y_ticks = np.arange(0, int(y_min) + 50, 50)
+        curr_ax.set_yticks(y_ticks)
+
         cb = fig.colorbar(cntrf, ax=curr_ax)
         cb.locator = matplotlib.ticker.LinearLocator(numticks=3)
         cb.minorticks_off()
@@ -425,16 +579,19 @@ def plot_2column(self) -> str:
         scfac = max(idist) / max(iz)  # noqa: F841
 
         fig, ax = plt.subplots(nrows=5, ncols=2, figsize=(18, 10))
-        fig.tight_layout(rect=[0, 0.03, 1, 0.96])
+        plt.subplots_adjust(hspace=0.15, wspace=0.1, left=0.05, right=0.98, top=0.96, bottom=0.03)
 
         best_ctd = self._get_best_ctd()
 
         # Compute density (sigma-t) if not already present
         self._compute_density(best_ctd)
 
+        # Create map in top-left subplot (row=0, col=0), aligned with ax[1,0] below
+        self._plot_track_map(ax[0, 0], ax[1, 0])
+
         profile_bottoms = self._profile_bottoms(distnav)
-        row = 0
-        col = 1
+        row = 1  # Start at row 1, col 0 (below the map)
+        col = 0
         for var, scale in (
             ("density", "linear"),
             (f"{best_ctd}_temperature", "linear"),
@@ -442,45 +599,41 @@ def plot_2column(self) -> str:
             ("isus_nitrate", "linear"),
             ("ctd1_oxygen_mll", "linear"),
             ("hs2_bbp420", "linear"),
-            ("hs2_bb700", "linear"),
+            ("hs2_bbp700", "linear"),
             ("hs2_fl700", "linear"),
-            ("biolume_avg_biolume", "log"),
+            ("biolume_avg_biolume", "linear"),
         ):
             self.logger.info("Plotting %s...", var)
             if var not in self.ds:
-                self.logger.warning("%s not in dataset", var)
-                ax[row, col].get_xaxis().set_visible(False)
-                ax[row, col].get_yaxis().set_visible(False)
-            else:
-                self._plot_var(
-                    var,
-                    idist,
-                    iz,
-                    distnav,
-                    fig,
-                    ax,
-                    row,
-                    col,
-                    profile_bottoms,
-                    scale=scale,
-                )
+                self.logger.warning("%s not in dataset, plotting with no data", var)
+
+            self._plot_var(
+                var,
+                idist,
+                iz,
+                distnav,
+                fig,
+                ax,
+                row,
+                col,
+                profile_bottoms,
+                scale=scale,
+                best_ctd=best_ctd,
+            )
             if row != 4:  # noqa: PLR2004
                 ax[row, col].get_xaxis().set_visible(False)
             else:
                 # Add x-axis label only to bottom row
                 ax[row, col].set_xlabel("Distance (km)")
 
-            if col == 1:
-                row += 1
-                col = 0
-            else:
+            # Column-major order: fill down first column, then second column
+            if row == 4 and col == 0:  # noqa: PLR2004
+                # Finished first column, move to top of second column
+                row = 0
                 col = 1
-
-        # Add title to the figure
-        title = f"{self.auv_name} {self.mission}"
-        if "title" in self.ds.attrs:
-            title = self.ds.attrs["title"]
-        fig.suptitle(title, fontsize=12, fontweight="bold")
+            else:
+                # Move down in current column
+                row += 1
 
         # Save plot to file
         images_dir = Path(BASE_PATH, self.auv_name, MISSIONIMAGES, self.mission)

From 25e9694d877fa6812d904f769c3906128ca8a506 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 07:37:41 -0800
Subject: [PATCH 04/20] Add y-axis label for variable with no data.

---
 src/data/create_products.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index f51143ff..0a4a7d1c 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -431,6 +431,38 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
                 fontweight="bold",
                 transform=curr_ax.transAxes,
             )
+
+            # Add a fake colorbar to maintain layout consistency
+            # Create a dummy mappable with a simple gray colormap
+            dummy_data = np.array([[0, 1]])
+            dummy_im = curr_ax.imshow(
+                dummy_data,
+                cmap="gray",
+                aspect="auto",
+                visible=False,
+            )
+            cb = fig.colorbar(dummy_im, ax=curr_ax)
+
+            # Hide the colorbar ticks and tick labels but keep the label visible
+            cb.set_ticks([])
+            cb.ax.set_facecolor("white")
+            cb.outline.set_visible(False)
+
+            # Get variable name for the label
+            if var in self.ds:
+                long_name = self.ds[var].attrs.get("long_name", var)
+                units = self.ds[var].attrs.get("units", "")
+            else:
+                # Extract readable name from variable string (e.g., "isus_nitrate" -> "Nitrate")
+                long_name = var.split("_")[-1].capitalize()
+                units = ""
+
+            # Add label to the colorbar area
+            if units:
+                cb.set_label(f"{long_name} [{units}]", fontsize=9)
+            else:
+                cb.set_label(long_name, fontsize=9)
+
             return
 
         # Normal plotting path - we have valid data

From 5e2b388f5ff323c7ea0aab8fbd5f3956574ed2f4 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 08:55:38 -0800
Subject: [PATCH 05/20] Add standard_names for dissolved oxygen, correct
 spelling of Bioluminescence.

---
 src/data/calibrate.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/data/calibrate.py b/src/data/calibrate.py
index 272ed731..b0e118e0 100755
--- a/src/data/calibrate.py
+++ b/src/data/calibrate.py
@@ -2332,6 +2332,7 @@ def _calibrated_oxygen(  # noqa: PLR0913
         oxygen_mll.attrs = {
             "long_name": "Dissolved Oxygen",
             "units": "ml/l",
+            "standard_name": "volume_fraction_of_oxygen_in_sea_water",
             "comment": mll_comment,
         }
 
@@ -2344,6 +2345,7 @@ def _calibrated_oxygen(  # noqa: PLR0913
         oxygen_umolkg.attrs = {
             "long_name": "Dissolved Oxygen",
             "units": "umol/kg",
+            "standard_name": "moles_of_oxygen_per_unit_mass_in_sea_water",
             "comment": umolkg_comment,
         }
         return oxygen_mll, oxygen_umolkg
@@ -2942,7 +2944,7 @@ def _biolume_process(self, sensor):
             name=f"{sensor}_flow",
         )
         self.combined_nc["biolume_flow"].attrs = {
-            "long_name": "Bioluminesence pump flow rate",
+            "long_name": "Bioluminescence pump flow rate",
             "units": "mL/s",
             "coordinates": f"{sensor}_time {sensor}_depth",
             "comment": f"flow from {source}",
@@ -2960,7 +2962,7 @@ def _biolume_process(self, sensor):
             name=f"{sensor}_avg_biolume",
         )
         self.combined_nc["biolume_avg_biolume"].attrs = {
-            "long_name": "Bioluminesence Average of 60Hz data",
+            "long_name": "Bioluminescence Average of 60Hz data",
             "units": "photons s^-1",
             "coordinates": f"{sensor}_{TIME} {sensor}_depth",
             "comment": f"avg_biolume from {source} {lag_info}",

From 218cab72585f2103e3d558bb03202e6b1f88b756 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 09:35:11 -0800
Subject: [PATCH 06/20] Use correct spelling of "Bioluminescence".

---
 src/data/combine.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/data/combine.py b/src/data/combine.py
index 791e26ba..2bf9cd60 100755
--- a/src/data/combine.py
+++ b/src/data/combine.py
@@ -317,7 +317,7 @@ def _biolume_process(self, sensor):
             name=f"{sensor}_flow",
         )
         self.combined_nc["biolume_flow"].attrs = {
-            "long_name": "Bioluminesence pump flow rate",
+            "long_name": "Bioluminescence pump flow rate",
             "units": "mL/s",
             "coordinates": f"{sensor}_time {sensor}_depth",
             "comment": f"flow from {source}",
@@ -335,7 +335,7 @@ def _biolume_process(self, sensor):
             name=f"{sensor}_avg_biolume",
         )
         self.combined_nc["biolume_avg_biolume"].attrs = {
-            "long_name": "Bioluminesence Average of 60Hz data",
+            "long_name": "Bioluminescence Average of 60Hz data",
             "units": "photons s^-1",
             "coordinates": f"{sensor}_{TIME} {sensor}_depth",
             "comment": f"avg_biolume from {source} {lag_info}",

From 9382ba9373c1983a8fc4264bb2b6145e77bc9f65 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 09:44:10 -0800
Subject: [PATCH 07/20] Improve plot labels and colormap selection

- Add MAX_LONG_NAME_LENGTH constant (40 chars) to control when variable
  names are used instead of long_name attributes in colorbar labels
- Update x-axis labels to "Distance along track (km)" for clarity
- Increase bottom margins from 0.03 to 0.06 to ensure x-axis labels
  are fully visible in both 2-column and biolume plots
- Add variable_colormap_lookup dict for fallback colormap selection
  when standard_name attribute is not available
- Add _get_colormap_name() helper method that tries standard_name
  first, then matches variable name patterns, then defaults to 'cividis'
- Simplify colormap selection logic by using the new helper method

These changes improve plot readability and make colormap selection
more robust for variables without standard_name attributes.
---
 .vscode/launch.json         |  4 +-
 src/data/create_products.py | 77 +++++++++++++++++++++++++++++--------
 2 files changed, 64 insertions(+), 17 deletions(-)

diff --git a/.vscode/launch.json b/.vscode/launch.json
index 1f8dcd2a..c89e6b7c 100644
--- a/.vscode/launch.json
+++ b/.vscode/launch.json
@@ -294,7 +294,7 @@
             //"args": ["-v", "1", "--last_n_days", "30", "--start_year", "2022", "--end_year", "2022", "--resample", "--noinput"]
             //"args": ["-v", "1", "--last_n_days", "30", "--start_year", "2022", "--end_year", "2022", "--archive", "--noinput"]
             //"args": ["-v", "1", "--start_year", "2021", "--end_year", "2022",  "--noinput"]
-            //"args": ["-v", "1", "--mission", "2022.201.00", "--clobber", "--noinput", "--no_cleanup"]
+            "args": ["-v", "1", "--mission", "2022.201.00", "--clobber", "--noinput", "--no_cleanup"]
             //"args": ["-v", "1", "--mission", "2009.084.00",  "--clobber", "--noinput"]
             //"args": ["-v", "1", "--mission", "2022.243.00",  "--calibrate", "--clobber", "--noinput"]
             //"args": ["-v", "1", "--start_year", "2004",  "--end_year", "2004", "--start_yd", "168", "--use_portal", "--clobber", "--noinput"]
@@ -338,7 +338,7 @@
             //"args": ["-v", "2", "--mission", "2004.029.03", "--noinput", "--no_cleanup"],
             //"args": ["-v", "1", "--mission", "2023.192.01", "--noinput", "--no_cleanup"],
             //"args": ["-v", "1", "--mission", "2010.151.04", "--noinput", "--no_cleanup", "--clobber"],
-            "args": ["-v", "1", "--mission", "2025.316.02", "--noinput", "--no_cleanup", "--add_seconds", "619315200"],
+            //"args": ["-v", "1", "--mission", "2025.316.02", "--noinput", "--no_cleanup", "--add_seconds", "619315200"],
 
         },
         {
diff --git a/src/data/create_products.py b/src/data/create_products.py
index 0a4a7d1c..346b42b0 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -8,7 +8,6 @@
 __copyright__ = "Copyright 2023, Monterey Bay Aquarium Research Institute"
 
 import argparse  # noqa: I001
-import contextlib
 import logging
 import os
 import re
@@ -74,6 +73,9 @@ def __init__(  # noqa: PLR0913
         self.verbose = verbose
         self.commandline = commandline
 
+    # Maximum length for long_name before using variable name instead
+    MAX_LONG_NAME_LENGTH = 40
+
     # Column name format required by ODV - will be tab delimited
     ODV_COLUMN_NAMES = [  # noqa: RUF012
         "Cruise",
@@ -122,6 +124,20 @@ def __init__(  # noqa: PLR0913
         "eastward_sea_water_velocity": "balance",
         "northward_wind": "balance",
         "eastward_wind": "balance",
+        "volume_fraction_of_oxygen_in_sea_water": "oxy",
+        "moles_of_oxygen_per_unit_mass_in_sea_water": "oxy",
+    }
+    # Fallback colormap lookup by variable name
+    variable_colormap_lookup = {  # noqa: RUF012
+        "nitrate": "matter",
+        "hs2_bbp420": "matter",
+        "hs2_bbp700": "matter",
+        "hs2_bbp470": "matter",
+        "hs2_bbp676": "matter",
+        "hs2_fl700": "algae",
+        "hs2_fl676": "algae",
+        "ecopuck_chla": "algae",
+        "biolume_avg_biolume": "cividis",
     }
 
     def _open_ds(self):
@@ -275,6 +291,35 @@ def _profile_bottoms(
         window = int(len(distnav) * window_frac)
         return depth_dist.rolling(dist=window).max()
 
+    def _get_colormap_name(self, var: str) -> str:
+        """Get colormap name for a variable.
+
+        Tries in order:
+        1. Lookup by standard_name attribute
+        2. Lookup by matching variable name parts
+        3. Default to 'cividis'
+
+        Args:
+            var: Variable name
+
+        Returns:
+            Colormap name
+        """
+        # Try standard_name first
+        if var in self.ds and "standard_name" in self.ds[var].attrs:
+            standard_name = self.ds[var].attrs["standard_name"]
+            if standard_name in self.cmocean_lookup:
+                return self.cmocean_lookup[standard_name]
+
+        # Fallback: try matching variable name parts
+        var_lower = var.lower()
+        for key, colormap in self.variable_colormap_lookup.items():
+            if key in var_lower:
+                return colormap
+
+        # Default
+        return "cividis"
+
     def _plot_track_map(
         self, map_ax: matplotlib.axes.Axes, reference_ax: matplotlib.axes.Axes
     ) -> None:
@@ -385,7 +430,6 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
         profile_bottoms: xr.DataArray,
         scale: str = "linear",
         num_colors: int = 256,
-        best_ctd: str = "ctd1",
     ):
         # Handle both 1D and 2D axis arrays
         curr_ax = ax[row] if col == 0 and hasattr(ax, "ndim") and ax.ndim == 1 else ax[row, col]
@@ -399,7 +443,8 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
             var_to_plot = (
                 np.log10(self.ds[var].to_numpy()) if scale == "log" else self.ds[var].to_numpy()
             )
-            valid_data = var_to_plot[~np.isnan(var_to_plot)]
+            # Filter out both NaN and infinite values (e.g., log10(0) = -inf)
+            valid_data = var_to_plot[~np.isnan(var_to_plot) & ~np.isinf(var_to_plot)]
             if len(valid_data) == 0:
                 self.logger.warning("%s has no valid data", var)
                 no_data = True
@@ -457,6 +502,10 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
                 long_name = var.split("_")[-1].capitalize()
                 units = ""
 
+            # Use variable name if long_name is too long
+            if len(long_name) > self.MAX_LONG_NAME_LENGTH:
+                long_name = var
+
             # Add label to the colorbar area
             if units:
                 cb.set_label(f"{long_name} [{units}]", fontsize=9)
@@ -475,12 +524,7 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
             rescale=True,
         )
 
-        color_map_name = "cividis"
-        with contextlib.suppress(KeyError):
-            color_map_name = self.cmocean_lookup.get(
-                self.ds[var].attrs["standard_name"],
-                "cividis",
-            )
+        color_map_name = self._get_colormap_name(var)
         try:
             cmap = plt.get_cmap(color_map_name)
         except ValueError:
@@ -591,6 +635,10 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
         long_name = self.ds[var].attrs.get("long_name", var)
         units = self.ds[var].attrs.get("units", "")
 
+        # Use variable name if long_name is too long
+        if len(long_name) > self.MAX_LONG_NAME_LENGTH:
+            long_name = var
+
         if scale == "log" and units:
             cb.set_label(f"{long_name}\n[log10({units})]", fontsize=7)
         elif scale == "log":
@@ -611,7 +659,7 @@ def plot_2column(self) -> str:
         scfac = max(idist) / max(iz)  # noqa: F841
 
         fig, ax = plt.subplots(nrows=5, ncols=2, figsize=(18, 10))
-        plt.subplots_adjust(hspace=0.15, wspace=0.1, left=0.05, right=0.98, top=0.96, bottom=0.03)
+        plt.subplots_adjust(hspace=0.15, wspace=0.1, left=0.05, right=0.98, top=0.96, bottom=0.06)
 
         best_ctd = self._get_best_ctd()
 
@@ -633,7 +681,7 @@ def plot_2column(self) -> str:
             ("hs2_bbp420", "linear"),
             ("hs2_bbp700", "linear"),
             ("hs2_fl700", "linear"),
-            ("biolume_avg_biolume", "linear"),
+            ("biolume_avg_biolume", "log"),
         ):
             self.logger.info("Plotting %s...", var)
             if var not in self.ds:
@@ -650,13 +698,12 @@ def plot_2column(self) -> str:
                 col,
                 profile_bottoms,
                 scale=scale,
-                best_ctd=best_ctd,
             )
             if row != 4:  # noqa: PLR2004
                 ax[row, col].get_xaxis().set_visible(False)
             else:
                 # Add x-axis label only to bottom row
-                ax[row, col].set_xlabel("Distance (km)")
+                ax[row, col].set_xlabel("Distance along track (km)")
 
             # Column-major order: fill down first column, then second column
             if row == 4 and col == 0:  # noqa: PLR2004
@@ -700,7 +747,7 @@ def plot_biolume(self) -> str:
         fig, ax = plt.subplots(nrows=num_plots, ncols=1, figsize=(18, 12))
         if num_plots == 1:
             ax = [ax]
-        fig.tight_layout(rect=[0, 0.03, 1, 0.96])
+        fig.tight_layout(rect=[0, 0.06, 1, 0.96])
 
         # Priority order for biolume variables to plot
         priority_vars = [
@@ -737,7 +784,7 @@ def plot_biolume(self) -> str:
             if i != num_plots - 1:
                 ax[i].get_xaxis().set_visible(False)
             else:
-                ax[i].set_xlabel("Distance (km)")
+                ax[i].set_xlabel("Distance along track (km)")
 
         # Add title to the figure
         title = f"{self.auv_name} {self.mission} - Bioluminescence"

From 70cdf075dbb17e0b3c82f0c0c88581cf341fdf18 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 10:06:55 -0800
Subject: [PATCH 08/20] Make subplots wider.

---
 src/data/create_products.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index 346b42b0..983a5199 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -659,7 +659,7 @@ def plot_2column(self) -> str:
         scfac = max(idist) / max(iz)  # noqa: F841
 
         fig, ax = plt.subplots(nrows=5, ncols=2, figsize=(18, 10))
-        plt.subplots_adjust(hspace=0.15, wspace=0.1, left=0.05, right=0.98, top=0.96, bottom=0.06)
+        plt.subplots_adjust(hspace=0.15, wspace=0.01, left=0.05, right=1.01, top=0.96, bottom=0.06)
 
         best_ctd = self._get_best_ctd()
 
@@ -747,7 +747,7 @@ def plot_biolume(self) -> str:
         fig, ax = plt.subplots(nrows=num_plots, ncols=1, figsize=(18, 12))
         if num_plots == 1:
             ax = [ax]
-        fig.tight_layout(rect=[0, 0.06, 1, 0.96])
+        fig.tight_layout(rect=[0, 0.06, 0.99, 0.96])
 
         # Priority order for biolume variables to plot
         priority_vars = [

From 011033109aaf3dc686c0b0ce2c5a86f8a48e8a3c Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 10:39:55 -0800
Subject: [PATCH 09/20] Use Blues and Reds for hs2 backscatter and exponential
 notation for all hs2_ variables.

---
 src/data/create_products.py | 31 ++++++++++++++++++++++++++-----
 1 file changed, 26 insertions(+), 5 deletions(-)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index 983a5199..06f7519a 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -130,10 +130,10 @@ def __init__(  # noqa: PLR0913
     # Fallback colormap lookup by variable name
     variable_colormap_lookup = {  # noqa: RUF012
         "nitrate": "matter",
-        "hs2_bbp420": "matter",
-        "hs2_bbp700": "matter",
-        "hs2_bbp470": "matter",
-        "hs2_bbp676": "matter",
+        "hs2_bbp420": "Blues",
+        "hs2_bbp700": "Reds",
+        "hs2_bbp470": "Blues",
+        "hs2_bbp676": "Reds",
         "hs2_fl700": "algae",
         "hs2_fl676": "algae",
         "ecopuck_chla": "algae",
@@ -629,7 +629,28 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
         cb.locator = matplotlib.ticker.LinearLocator(numticks=3)
         cb.minorticks_off()
         cb.update_ticks()
-        cb.ax.set_yticklabels([f"{x:.1f}" for x in cb.get_ticks()])
+
+        # Use scientific notation with offset for hs2 variables
+        if var.startswith("hs2_"):
+            # Calculate the order of magnitude
+            tick_values = cb.get_ticks()
+            max_abs = max(abs(tick_values.min()), abs(tick_values.max()))
+            if max_abs > 0:
+                order = int(np.floor(np.log10(max_abs)))
+                scale = 10**order
+                # Set clean tick labels
+                cb.ax.set_yticklabels([f"{x / scale:.2f}" for x in tick_values])
+                # Add offset text
+                cb.ax.text(
+                    1.5,
+                    1.05,
+                    f"10$^{{{order}}}$",
+                    transform=cb.ax.transAxes,
+                    fontsize=9,
+                    verticalalignment="bottom",
+                )
+        else:
+            cb.ax.set_yticklabels([f"{x:.1f}" for x in cb.get_ticks()])
 
         # Get long_name and units with fallbacks
         long_name = self.ds[var].attrs.get("long_name", var)

From 114528b4d0a1ea444eaff6fe012e7a76b339460a Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 10:57:58 -0800
Subject: [PATCH 10/20] Add measurement and Gulper number locations.

---
 src/data/create_products.py | 65 +++++++++++++++++++++++++++++++++++++
 1 file changed, 65 insertions(+)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index 06f7519a..87736ab7 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -291,6 +291,40 @@ def _profile_bottoms(
         window = int(len(distnav) * window_frac)
         return depth_dist.rolling(dist=window).max()
 
+    def _get_gulper_locations(self, distnav: xr.DataArray) -> dict:
+        """Get gulper bottle locations in distance/depth space.
+
+        Returns:
+            Dictionary mapping bottle number to (distance_km, depth_m) tuple
+        """
+        gulper = Gulper()
+        gulper.args = argparse.Namespace()
+        gulper.args.base_path = self.base_path
+        gulper.args.auv_name = self.auv_name
+        gulper.args.mission = self.mission
+        gulper.args.local = self.local
+        gulper.args.verbose = 0  # Suppress gulper logging
+        gulper.args.start_esecs = self.start_esecs
+        gulper.logger.setLevel(logging.WARNING)
+
+        gulper_times = gulper.parse_gulpers()
+        if not gulper_times:
+            return {}
+
+        locations = {}
+        for bottle, esec in gulper_times.items():
+            # Find closest time index
+            time_ns = np.datetime64(int(esec * 1e9), "ns")
+            time_idx = np.abs(self.ds.cf["time"].to_numpy() - time_ns).argmin()
+
+            # Get distance and depth at that time
+            dist_km = distnav.to_numpy()[time_idx] / 1000.0
+            depth_m = self.ds.cf["depth"].to_numpy()[time_idx]
+
+            locations[bottle] = (dist_km, depth_m)
+
+        return locations
+
     def _get_colormap_name(self, var: str) -> str:
         """Get colormap name for a variable.
 
@@ -430,6 +464,7 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
         profile_bottoms: xr.DataArray,
         scale: str = "linear",
         num_colors: int = 256,
+        gulper_locations: dict = None,
     ):
         # Handle both 1D and 2D axis arrays
         curr_ax = ax[row] if col == 0 and hasattr(ax, "ndim") and ax.ndim == 1 else ax[row, col]
@@ -613,6 +648,32 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
         )
         curr_ax.fill(list(reversed(xb)), list(reversed(yb)), "w")
 
+        # Add measurement location dots for density plot
+        if var == "density":
+            curr_ax.scatter(
+                distnav.to_numpy() / 1000.0,
+                self.ds.cf["depth"].to_numpy(),
+                s=0.1,
+                c="black",
+                alpha=0.1,
+                zorder=3,
+            )
+
+        # Add gulper bottle locations
+        if gulper_locations:
+            for bottle, (dist, depth) in gulper_locations.items():
+                curr_ax.text(
+                    dist,
+                    depth - 5,
+                    str(bottle),
+                    fontsize=7,
+                    ha="center",
+                    va="top",
+                    color="black",
+                    fontweight="bold",
+                    zorder=5,
+                )
+
         # Only show y-label on left column or top of right column
         if col == 0 or (col == 1 and row == 0):
             curr_ax.set_ylabel("Depth (m)")
@@ -690,6 +751,9 @@ def plot_2column(self) -> str:
         # Create map in top-left subplot (row=0, col=0), aligned with ax[1,0] below
         self._plot_track_map(ax[0, 0], ax[1, 0])
 
+        # Parse gulper locations
+        gulper_locations = self._get_gulper_locations(distnav)
+
         profile_bottoms = self._profile_bottoms(distnav)
         row = 1  # Start at row 1, col 0 (below the map)
         col = 0
@@ -719,6 +783,7 @@ def plot_2column(self) -> str:
                 col,
                 profile_bottoms,
                 scale=scale,
+                gulper_locations=gulper_locations,
             )
             if row != 4:  # noqa: PLR2004
                 ax[row, col].get_xaxis().set_visible(False)

From 3044f9a99c4ab2c7ca4979045a0a1b2de5f4d769 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 11:10:01 -0800
Subject: [PATCH 11/20] Use contextily to add OSM basemap to the mapplot.

---
 pyproject.toml              |   1 +
 src/data/create_products.py |  37 +++++++---
 uv.lock                     | 132 ++++++++++++++++++++++++++++++++++++
 3 files changed, 161 insertions(+), 9 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index cb55b1da..1e164196 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -11,6 +11,7 @@ dependencies = [
     "cf-xarray>=0.10.5",
     "cmocean>=4.0.3",
     "coards>=1.0.5",
+    "contextily>=1.7.0",
     "datashader>=0.18.1",
     "defusedxml>=0.7.1",
     "gitpython>=3.1.44",
diff --git a/src/data/create_products.py b/src/data/create_products.py
index 87736ab7..3577190e 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -16,6 +16,7 @@
 from pathlib import Path
 
 import cmocean
+import contextily as ctx
 import gsw
 import matplotlib  # noqa: ICN001
 import matplotlib.pyplot as plt
@@ -367,6 +368,12 @@ def _plot_track_map(
         lons = self.ds.cf["longitude"].to_numpy()
         lats = self.ds.cf["latitude"].to_numpy()
 
+        # Convert to Web Mercator for contextily
+        import pyproj
+
+        transformer = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
+        x_merc, y_merc = transformer.transform(lons, lats)
+
         # Get time data for start/end
         times = self.ds.cf["time"].to_numpy()
         start_time = pd.to_datetime(times[0]).strftime("%Y-%m-%d %H:%M:%S")
@@ -384,11 +391,11 @@ def _plot_track_map(
         # Make the plot square by using equal aspect (do this first)
         map_ax.set_aspect("equal", adjustable="datalim")
 
-        # Plot the track with depth coloring
+        # Plot the track with depth coloring in Web Mercator coordinates
         depths = self.ds.cf["depth"].to_numpy()
         map_ax.scatter(
-            lons,
-            lats,
+            x_merc,
+            y_merc,
             c=depths,
             cmap="viridis_r",  # Reversed so shallow is yellow, deep is dark
             s=1,
@@ -396,12 +403,24 @@ def _plot_track_map(
         )
 
         # Add start and end markers
-        map_ax.plot(lons[0], lats[0], "go", markersize=8, label="Start", zorder=5)
-        map_ax.plot(lons[-1], lats[-1], "r^", markersize=8, label="End", zorder=5)
-
-        # Set fixed axis limits for Monterey Bay area
-        map_ax.set_xlim([-122.41, -121.77])
-        map_ax.set_ylim([36.5, 37.0])
+        map_ax.plot(x_merc[0], y_merc[0], "go", markersize=8, label="Start", zorder=5)
+        map_ax.plot(x_merc[-1], y_merc[-1], "r^", markersize=8, label="End", zorder=5)
+
+        # Set fixed axis limits for Monterey Bay area (in Web Mercator)
+        lon_bounds = [-122.41, -121.77]
+        lat_bounds = [36.5, 37.0]
+        x_bounds, y_bounds = transformer.transform(lon_bounds, lat_bounds)
+        map_ax.set_xlim(x_bounds)
+        map_ax.set_ylim(y_bounds)
+
+        # Add basemap
+        ctx.add_basemap(
+            map_ax,
+            crs="EPSG:3857",
+            source=ctx.providers.OpenStreetMap.Mapnik,
+            alpha=0.6,
+            zorder=0,
+        )
 
         # Now position map aligned with left edge of reference, 50% width
         # Use a square aspect ratio based on the y-dimension
diff --git a/uv.lock b/uv.lock
index 82bbff19..bbfdc67a 100644
--- a/uv.lock
+++ b/uv.lock
@@ -2,6 +2,15 @@ version = 1
 revision = 3
 requires-python = "==3.12.*"
 
+[[package]]
+name = "affine"
+version = "2.4.0"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/69/98/d2f0bb06385069e799fc7d2870d9e078cfa0fa396dc8a2b81227d0da08b9/affine-2.4.0.tar.gz", hash = "sha256:a24d818d6a836c131976d22f8c27b8d3ca32d0af64c1d8d29deb7bafa4da1eea", size = 17132, upload-time = "2023-01-19T23:44:30.696Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/0b/f7/85273299ab57117850cc0a936c64151171fac4da49bc6fba0dad984a7c5f/affine-2.4.0-py3-none-any.whl", hash = "sha256:8a3df80e2b2378aef598a83c1392efd47967afec4242021a0b06b4c7cbc61a92", size = 15662, upload-time = "2023-01-19T23:44:28.833Z" },
+]
+
 [[package]]
 name = "aiofiles"
 version = "24.1.0"
@@ -172,6 +181,7 @@ dependencies = [
     { name = "cf-xarray" },
     { name = "cmocean" },
     { name = "coards" },
+    { name = "contextily" },
     { name = "datashader" },
     { name = "defusedxml" },
     { name = "gitpython" },
@@ -204,6 +214,7 @@ requires-dist = [
     { name = "cf-xarray", specifier = ">=0.10.5" },
     { name = "cmocean", specifier = ">=4.0.3" },
     { name = "coards", specifier = ">=1.0.5" },
+    { name = "contextily", specifier = ">=1.7.0" },
     { name = "datashader", specifier = ">=0.18.1" },
     { name = "defusedxml", specifier = ">=0.7.1" },
     { name = "gitpython", specifier = ">=3.1.44" },
@@ -380,6 +391,42 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/20/94/c5790835a017658cbfabd07f3bfb549140c3ac458cfc196323996b10095a/charset_normalizer-3.4.2-py3-none-any.whl", hash = "sha256:7f56930ab0abd1c45cd15be65cc741c28b1c9a34876ce8c17a2fa107810c0af0", size = 52626, upload-time = "2025-05-02T08:34:40.053Z" },
 ]
 
+[[package]]
+name = "click"
+version = "8.3.1"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "colorama", marker = "sys_platform == 'win32'" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/3d/fa/656b739db8587d7b5dfa22e22ed02566950fbfbcdc20311993483657a5c0/click-8.3.1.tar.gz", hash = "sha256:12ff4785d337a1bb490bb7e9c2b1ee5da3112e94a8622f26a6c77f5d2fc6842a", size = 295065, upload-time = "2025-11-15T20:45:42.706Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/98/78/01c019cdb5d6498122777c1a43056ebb3ebfeef2076d9d026bfe15583b2b/click-8.3.1-py3-none-any.whl", hash = "sha256:981153a64e25f12d547d3426c367a4857371575ee7ad18df2a6183ab0545b2a6", size = 108274, upload-time = "2025-11-15T20:45:41.139Z" },
+]
+
+[[package]]
+name = "click-plugins"
+version = "1.1.1.2"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "click" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/c3/a4/34847b59150da33690a36da3681d6bbc2ec14ee9a846bc30a6746e5984e4/click_plugins-1.1.1.2.tar.gz", hash = "sha256:d7af3984a99d243c131aa1a828331e7630f4a88a9741fd05c927b204bcf92261", size = 8343, upload-time = "2025-06-25T00:47:37.555Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/3d/9a/2abecb28ae875e39c8cad711eb1186d8d14eab564705325e77e4e6ab9ae5/click_plugins-1.1.1.2-py2.py3-none-any.whl", hash = "sha256:008d65743833ffc1f5417bf0e78e8d2c23aab04d9745ba817bd3e71b0feb6aa6", size = 11051, upload-time = "2025-06-25T00:47:36.731Z" },
+]
+
+[[package]]
+name = "cligj"
+version = "0.7.2"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "click" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/ea/0d/837dbd5d8430fd0f01ed72c4cfb2f548180f4c68c635df84ce87956cff32/cligj-0.7.2.tar.gz", hash = "sha256:a4bc13d623356b373c2c27c53dbd9c68cae5d526270bfa71f6c6fa69669c6b27", size = 9803, upload-time = "2021-05-28T21:23:27.935Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/73/86/43fa9f15c5b9fb6e82620428827cd3c284aa933431405d1bcf5231ae3d3e/cligj-0.7.2-py3-none-any.whl", hash = "sha256:c1ca117dbce1fe20a5809dc96f01e1c2840f6dcc939b3ddbb1111bf330ba82df", size = 7069, upload-time = "2021-05-28T21:23:26.877Z" },
+]
+
 [[package]]
 name = "cmocean"
 version = "4.0.3"
@@ -430,6 +477,25 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/e6/75/49e5bfe642f71f272236b5b2d2691cf915a7283cc0ceda56357b61daa538/comm-0.2.2-py3-none-any.whl", hash = "sha256:e6fb86cb70ff661ee8c9c14e7d36d6de3b4066f1441be4063df9c5009f0a64d3", size = 7180, upload-time = "2024-03-12T16:53:39.226Z" },
 ]
 
+[[package]]
+name = "contextily"
+version = "1.7.0"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "geopy" },
+    { name = "joblib" },
+    { name = "matplotlib" },
+    { name = "mercantile" },
+    { name = "pillow" },
+    { name = "rasterio" },
+    { name = "requests" },
+    { name = "xyzservices" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/cc/f3/4c7fdb1ef92c8e3db1de9741595ae1815ddc8a7e2088539e6107c83e2268/contextily-1.7.0.tar.gz", hash = "sha256:6534faa5702b89b46d0d81b4c538754f2d8b3dd8cc298454b11ccedfa67e73ac", size = 22462157, upload-time = "2025-11-24T19:54:39.199Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/6a/bb/6c824d0da874eef900ecf0a71e0889f0d4624560b0160fd9e333a146ee4f/contextily-1.7.0-py3-none-any.whl", hash = "sha256:38393b8e7dc38580ef1f60211a9ba1c3eb142a439c260509c201987754fa9dba", size = 16849, upload-time = "2025-11-24T19:54:37.018Z" },
+]
+
 [[package]]
 name = "contourpy"
 version = "1.3.2"
@@ -602,6 +668,27 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/13/be/0ebbb283f2d91b72beaee2d07760b2c47dab875c49c286f5591d3d157198/frozenlist-1.6.2-py3-none-any.whl", hash = "sha256:947abfcc8c42a329bbda6df97a4b9c9cdb4e12c85153b3b57b9d2f02aa5877dc", size = 12582, upload-time = "2025-06-03T21:48:03.201Z" },
 ]
 
+[[package]]
+name = "geographiclib"
+version = "2.1"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/df/78/4892343230a9d29faa1364564e525307a37e54ad776ea62c12129dbba704/geographiclib-2.1.tar.gz", hash = "sha256:6a6545e6262d0ed3522e13c515713718797e37ed8c672c31ad7b249f372ef108", size = 37004, upload-time = "2025-08-21T21:34:26Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/31/b3/802576f2ea5dcb48501bb162e4c7b7b3ca5654a42b2c968ef98a797a4c79/geographiclib-2.1-py3-none-any.whl", hash = "sha256:e2a873b9b9e7fc38721ad73d5f4e6c9ed140d428a339970f505c07056997d40b", size = 40740, upload-time = "2025-08-21T21:34:24.955Z" },
+]
+
+[[package]]
+name = "geopy"
+version = "2.4.1"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "geographiclib" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/0e/fd/ef6d53875ceab72c1fad22dbed5ec1ad04eb378c2251a6a8024bad890c3b/geopy-2.4.1.tar.gz", hash = "sha256:50283d8e7ad07d89be5cb027338c6365a32044df3ae2556ad3f52f4840b3d0d1", size = 117625, upload-time = "2023-11-23T21:49:32.734Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/e5/15/cf2a69ade4b194aa524ac75112d5caac37414b20a3a03e6865dfe0bd1539/geopy-2.4.1-py3-none-any.whl", hash = "sha256:ae8b4bc5c1131820f4d75fce9d4aaaca0c85189b3aa5d64c3dcaf5e3b7b882a7", size = 125437, upload-time = "2023-11-23T21:49:30.421Z" },
+]
+
 [[package]]
 name = "gitdb"
 version = "4.0.12"
@@ -870,6 +957,15 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/62/a1/3d680cbfd5f4b8f15abc1d571870c5fc3e594bb582bc3b64ea099db13e56/jinja2-3.1.6-py3-none-any.whl", hash = "sha256:85ece4451f492d0c13c5dd7c13a64681a86afae63a5f347908daf103ce6d2f67", size = 134899, upload-time = "2025-03-05T20:05:00.369Z" },
 ]
 
+[[package]]
+name = "joblib"
+version = "1.5.3"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/41/f2/d34e8b3a08a9cc79a50b2208a93dce981fe615b64d5a4d4abee421d898df/joblib-1.5.3.tar.gz", hash = "sha256:8561a3269e6801106863fd0d6d84bb737be9e7631e33aaed3fb9ce5953688da3", size = 331603, upload-time = "2025-12-15T08:41:46.427Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/7b/91/984aca2ec129e2757d1e4e3c81c3fcda9d0f85b74670a094cc443d9ee949/joblib-1.5.3-py3-none-any.whl", hash = "sha256:5fc3c5039fc5ca8c0276333a188bbd59d6b7ab37fe6632daa76bc7f9ec18e713", size = 309071, upload-time = "2025-12-15T08:41:44.973Z" },
+]
+
 [[package]]
 name = "json5"
 version = "0.12.0"
@@ -1272,6 +1368,18 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/b3/38/89ba8ad64ae25be8de66a6d463314cf1eb366222074cfda9ee839c56a4b4/mdurl-0.1.2-py3-none-any.whl", hash = "sha256:84008a41e51615a49fc9966191ff91509e3c40b939176e643fd50a5c2196b8f8", size = 9979, upload-time = "2022-08-14T12:40:09.779Z" },
 ]
 
+[[package]]
+name = "mercantile"
+version = "1.2.1"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "click" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/d2/c6/87409bcb26fb68c393fa8cf58ba09363aa7298cfb438a0109b5cb14bc98b/mercantile-1.2.1.tar.gz", hash = "sha256:fa3c6db15daffd58454ac198b31887519a19caccee3f9d63d17ae7ff61b3b56b", size = 26352, upload-time = "2021-04-21T14:42:41.096Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/b2/d6/de0cc74f8d36976aeca0dd2e9cbf711882ff8e177495115fd82459afdc4d/mercantile-1.2.1-py3-none-any.whl", hash = "sha256:30f457a73ee88261aab787b7069d85961a5703bb09dc57a170190bc042cd023f", size = 14779, upload-time = "2021-04-21T14:42:39.841Z" },
+]
+
 [[package]]
 name = "mistune"
 version = "3.1.3"
@@ -1947,6 +2055,30 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/e4/29/073779afc3ef6f830b8de95026ef20b2d1ec22d0324d767748d806e57379/pyzmq-26.4.0-cp312-cp312-win_arm64.whl", hash = "sha256:2f23c750e485ce1eb639dbd576d27d168595908aa2d60b149e2d9e34c9df40e0", size = 556342, upload-time = "2025-04-04T12:03:59.218Z" },
 ]
 
+[[package]]
+name = "rasterio"
+version = "1.4.4"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "affine" },
+    { name = "attrs" },
+    { name = "certifi" },
+    { name = "click" },
+    { name = "click-plugins" },
+    { name = "cligj" },
+    { name = "numpy" },
+    { name = "pyparsing" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/ec/fa/fce8dc9f09e5bc6520b6fc1b4ecfa510af9ca06eb42ad7bdff9c9b8989d0/rasterio-1.4.4.tar.gz", hash = "sha256:c95424e2c7f009b8f7df1095d645c52895cd332c0c2e1b4c2e073ea28b930320", size = 445004, upload-time = "2025-12-12T18:01:08.971Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/4c/fc/63d89ddfcb4643730553683ee322566b9b15fe56d026e4c21c4f4f5d9d26/rasterio-1.4.4-cp312-cp312-macosx_14_0_arm64.whl", hash = "sha256:f3c4f0cbd188f893011f2a0a6dc2852b3892799b3a0d79eddf92f2b115ec7ed7", size = 21120715, upload-time = "2025-12-12T17:59:19.35Z" },
+    { url = "https://files.pythonhosted.org/packages/43/70/2c003f76a23dbb078fdee35c8e2ec490d2ad8982f4dc956ba08b56027b87/rasterio-1.4.4-cp312-cp312-macosx_15_0_x86_64.whl", hash = "sha256:6fce26090b9f509eab337228420145947c491a13628965410f25bc3e6e05cf75", size = 25732944, upload-time = "2025-12-12T17:59:22.533Z" },
+    { url = "https://files.pythonhosted.org/packages/f6/cc/4a8e92362c0ff496dd1007c3dcba66e9ededf1a45eca8ad1db302b071c49/rasterio-1.4.4-cp312-cp312-manylinux_2_28_aarch64.whl", hash = "sha256:c1c722da390dc264aeccdc0dc200ca37923875d910ca4cd5bec0fec351bb818e", size = 34295209, upload-time = "2025-12-12T17:59:26.035Z" },
+    { url = "https://files.pythonhosted.org/packages/e6/6d/717d2dec47fbefad33ca0d27bd5f0d543b1d1bc9fcab5ef82a13adaaf38d/rasterio-1.4.4-cp312-cp312-manylinux_2_28_x86_64.whl", hash = "sha256:98b6dfb8282b2a54b9d75c3dc8d2520a69bbc66916c7d43de8e0bbf6e0240ca1", size = 35661866, upload-time = "2025-12-12T17:59:29.928Z" },
+    { url = "https://files.pythonhosted.org/packages/ed/60/ae3351fba2726ec0976974ce2eb030c159edd3363b8771e832b8db571c24/rasterio-1.4.4-cp312-cp312-win_amd64.whl", hash = "sha256:9513f4c7a6d93b45098f8dff2421fa9516604e3bfbf35aa144484a88d36a321f", size = 25682853, upload-time = "2025-12-12T17:59:35.869Z" },
+    { url = "https://files.pythonhosted.org/packages/38/ee/35387296bbacfc5cbbb4273228b1b959793d3ce38b0402a07f11a248420b/rasterio-1.4.4-cp312-cp312-win_arm64.whl", hash = "sha256:60b49a482e0f12f12ce9d2cc3090add02f89f3d422e85f2cffaa9207adb83c04", size = 24043249, upload-time = "2025-12-12T17:59:39.915Z" },
+]
+
 [[package]]
 name = "referencing"
 version = "0.36.2"

From 40840073b4db4f466185639ecbf30f930a8500a8 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 11:32:31 -0800
Subject: [PATCH 12/20] Move map legend into title area.

---
 src/data/create_products.py | 67 +++++++++++++++++++++++++++++++++----
 1 file changed, 60 insertions(+), 7 deletions(-)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index 3577190e..bb01c6de 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -441,33 +441,86 @@ def _plot_track_map(
             spine.set_edgecolor("black")
             spine.set_linewidth(1)
 
-        # Add legend
-        map_ax.legend(loc="upper right", fontsize=7, framealpha=0.9)
-
         # Add text on the right side with title and times
         # Wrap title for better formatting
         import textwrap
 
         wrapped_title = textwrap.fill(title, width=40)
 
-        text_content = f"{wrapped_title}\n\nStart: {start_time}\nEnd:   {end_time}"
-
         # Get updated position after aspect adjustment
         updated_pos = map_ax.get_position()
 
         # Position text in figure coordinates, to the right of the map
         text_x = updated_pos.x0 + updated_pos.width + 0.01
-        text_y = updated_pos.y0 + updated_pos.height * 0.5
+        text_y = updated_pos.y0 + updated_pos.height  # Align with top of map
 
+        # Add title
         map_ax.figure.text(
             text_x,
             text_y,
-            text_content,
+            wrapped_title,
+            fontsize=11,
+            fontweight="bold",
+            family="sans-serif",
+            verticalalignment="top",
+            horizontalalignment="left",
+            color="black",
+        )
+
+        # Calculate approximate title height based on number of lines
+        # Each line is roughly 0.015 in figure coordinates at fontsize 11
+        num_title_lines = wrapped_title.count("\n") + 1
+        title_height_approx = num_title_lines * 0.02
+
+        # Position start/end text below the title with some spacing
+        start_end_y = text_y - title_height_approx - 0.03
+
+        # Add start marker and text
+        map_ax.figure.text(
+            text_x,
+            start_end_y,
+            "● ",
+            fontsize=14,
+            fontweight="bold",
+            family="sans-serif",
+            verticalalignment="center",
+            horizontalalignment="left",
+            color="green",
+        )
+        map_ax.figure.text(
+            text_x + 0.015,
+            start_end_y,
+            f"Start: {start_time}",
+            fontsize=11,
+            fontweight="bold",
+            family="sans-serif",
+            verticalalignment="center",
+            horizontalalignment="left",
+            color="black",
+        )
+
+        # Add end marker and text
+        map_ax.figure.text(
+            text_x,
+            start_end_y - 0.025,
+            "▲ ",
+            fontsize=14,
+            fontweight="bold",
+            family="sans-serif",
+            verticalalignment="center",
+            horizontalalignment="left",
+            color="red",
+        )
+        map_ax.figure.text(
+            text_x + 0.015,
+            start_end_y - 0.025,
+            f"End  : {end_time}",
             fontsize=11,
             fontweight="bold",
             family="sans-serif",
             verticalalignment="center",
             horizontalalignment="left",
+            color="black",
         )
 
     def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915

From 6d8da34d62117d74594cb954f97cc4c0111ec99c Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 11:38:47 -0800
Subject: [PATCH 13/20] Update testing values following metadata changes.

---
 src/data/test_process_dorado.py | 10 +++++-----
 src/data/test_process_i2map.py  |  6 +++---
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/data/test_process_dorado.py b/src/data/test_process_dorado.py
index 88217bbc..2e023ee1 100644
--- a/src/data/test_process_dorado.py
+++ b/src/data/test_process_dorado.py
@@ -31,9 +31,9 @@ def test_process_dorado(complete_dorado_processing):
     # but it will alert us if a code change unexpectedly changes the file size.
     # If code changes are expected to change the file size then we should
     # update the expected size here.
-    EXPECTED_SIZE_GITHUB = 621410
-    EXPECTED_SIZE_ACT = 621408
-    EXPECTED_SIZE_LOCAL = 621458
+    EXPECTED_SIZE_GITHUB = 627285
+    EXPECTED_SIZE_ACT = 627283
+    EXPECTED_SIZE_LOCAL = 627333
     if str(proc.args.base_path).startswith("/home/runner"):
         # The size is different in GitHub Actions, maybe due to different metadata
         assert nc_file.stat().st_size == EXPECTED_SIZE_GITHUB  # noqa: S101
@@ -51,8 +51,8 @@ def test_process_dorado(complete_dorado_processing):
     if check_md5:
         # Check that the MD5 hash has not changed
         EXPECTED_MD5_GITHUB = "fec067d16eb5280a8bc2b6ef132821b8"
-        EXPECTED_MD5_ACT = "316955fd489862ad9ed5b63df8aa7db8"
-        EXPECTED_MD5_LOCAL = "f635cee8760aa0f40bd2070bc0c5fa65"
+        EXPECTED_MD5_ACT = "1cd87ec363d3ca3323cf6df3f31ebe33"
+        EXPECTED_MD5_LOCAL = "bd72eed42be79af3bff7c2657b092416"
         if str(proc.args.base_path).startswith("/home/runner"):
             # The MD5 hash is different in GitHub Actions, maybe due to different metadata
             assert hashlib.md5(open(nc_file, "rb").read()).hexdigest() == EXPECTED_MD5_GITHUB  # noqa:  PTH123, S101, S324, SIM115
diff --git a/src/data/test_process_i2map.py b/src/data/test_process_i2map.py
index dc972932..09f3c827 100644
--- a/src/data/test_process_i2map.py
+++ b/src/data/test_process_i2map.py
@@ -30,9 +30,9 @@ def test_process_i2map(complete_i2map_processing):
     # but it will alert us if a code change unexpectedly changes the file size.
     # If code changes are expected to change the file size then we should
     # update the expected size here.
-    EXPECTED_SIZE_GITHUB = 52758
-    EXPECTED_SIZE_ACT = 52728
-    EXPECTED_SIZE_LOCAL = 52858
+    EXPECTED_SIZE_GITHUB = 63069
+    EXPECTED_SIZE_ACT = 63039
+    EXPECTED_SIZE_LOCAL = 64575
     if str(proc.args.base_path).startswith("/home/runner"):
         # The size is different in GitHub Actions, maybe due to different metadata
         assert nc_file.stat().st_size == EXPECTED_SIZE_GITHUB  # noqa: S101

From d837f179529bf50e3919f2832062b1b05938d51b Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 11:40:48 -0800
Subject: [PATCH 14/20] Use "BuPu" instead of "Blues". It's closer to 420nm.

---
 src/data/create_products.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index bb01c6de..8089bb9b 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -131,9 +131,9 @@ def __init__(  # noqa: PLR0913
     # Fallback colormap lookup by variable name
     variable_colormap_lookup = {  # noqa: RUF012
         "nitrate": "matter",
-        "hs2_bbp420": "Blues",
+        "hs2_bbp420": "BuPu",
         "hs2_bbp700": "Reds",
-        "hs2_bbp470": "Blues",
+        "hs2_bbp470": "BuPu",
         "hs2_bbp676": "Reds",
         "hs2_fl700": "algae",
         "hs2_fl676": "algae",

From beda89021e99d2d64fb8f2f188d551ff8c9e2148 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 11:41:38 -0800
Subject: [PATCH 15/20] Update EXPECTED_MD5_GITHUB.

---
 src/data/test_process_dorado.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/data/test_process_dorado.py b/src/data/test_process_dorado.py
index 2e023ee1..255bef4b 100644
--- a/src/data/test_process_dorado.py
+++ b/src/data/test_process_dorado.py
@@ -50,7 +50,7 @@ def test_process_dorado(complete_dorado_processing):
     check_md5 = True
     if check_md5:
         # Check that the MD5 hash has not changed
-        EXPECTED_MD5_GITHUB = "fec067d16eb5280a8bc2b6ef132821b8"
+        EXPECTED_MD5_GITHUB = "3fe4d2c1680db97289273905a386810c"
         EXPECTED_MD5_ACT = "1cd87ec363d3ca3323cf6df3f31ebe33"
         EXPECTED_MD5_LOCAL = "bd72eed42be79af3bff7c2657b092416"
         if str(proc.args.base_path).startswith("/home/runner"):

From 58395542c510cd9056915dad16d6349980006bb6 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 14:51:33 -0800
Subject: [PATCH 16/20] Add installation of GMT library as prerequisite for
 pygmt.

MacOS assumes installed via mac ports.
---
 docker/Dockerfile |  7 ++++++-
 pyproject.toml    |  1 +
 uv.lock           | 17 +++++++++++++++++
 3 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/docker/Dockerfile b/docker/Dockerfile
index 50363723..15af84f6 100644
--- a/docker/Dockerfile
+++ b/docker/Dockerfile
@@ -17,7 +17,12 @@ RUN apt-get update \
         curl \
         build-essential \
         git \
-        rsync
+        rsync \
+        # GMT6 dependencies
+        gmt \
+        gmt-dcw \
+        gmt-gshhg \
+    && rm -rf /var/lib/apt/lists/*
 
 WORKDIR $PYSETUP_PATH
 
diff --git a/pyproject.toml b/pyproject.toml
index 1e164196..90bd7e19 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -23,6 +23,7 @@ dependencies = [
     "numpy>=2.2.6",
     "pooch>=1.8.2",
     "pyarrow>=20.0.0",
+    "pygmt>=0.17.0",
     "pyproj>=3.7.1",
     "pysolar>=0.13",
     "rolling>=0.5.0",
diff --git a/uv.lock b/uv.lock
index bbfdc67a..0bbf808a 100644
--- a/uv.lock
+++ b/uv.lock
@@ -193,6 +193,7 @@ dependencies = [
     { name = "numpy" },
     { name = "pooch" },
     { name = "pyarrow" },
+    { name = "pygmt" },
     { name = "pyproj" },
     { name = "pysolar" },
     { name = "rolling" },
@@ -226,6 +227,7 @@ requires-dist = [
     { name = "numpy", specifier = ">=2.2.6" },
     { name = "pooch", specifier = ">=1.8.2" },
     { name = "pyarrow", specifier = ">=20.0.0" },
+    { name = "pygmt", specifier = ">=0.17.0" },
     { name = "pyproj", specifier = ">=3.7.1" },
     { name = "pysolar", specifier = ">=0.13" },
     { name = "rolling", specifier = ">=0.5.0" },
@@ -1899,6 +1901,21 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/8a/0b/9fcc47d19c48b59121088dd6da2488a49d5f72dacf8262e2790a1d2c7d15/pygments-2.19.1-py3-none-any.whl", hash = "sha256:9ea1544ad55cecf4b8242fab6dd35a93bbce657034b0611ee383099054ab6d8c", size = 1225293, upload-time = "2025-01-06T17:26:25.553Z" },
 ]
 
+[[package]]
+name = "pygmt"
+version = "0.17.0"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "numpy" },
+    { name = "packaging" },
+    { name = "pandas" },
+    { name = "xarray" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/c0/6e/a973f55dbe085713a0f705561fc4214593e11a2e4d25248f6502511bec7b/pygmt-0.17.0.tar.gz", hash = "sha256:06d7581bc35641db318c418aba02163bc519ba267a5a60b5cc419eee40390503", size = 228512, upload-time = "2025-10-02T22:51:40.781Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/11/e5/e249f42f47ca6ac7f88cad10d539334b04d77a76dc9d11b58b6a72b8b702/pygmt-0.17.0-py3-none-any.whl", hash = "sha256:846f93aef999be8f9d6de806e1aeebddc70ed59c33757c16bbb236678d0f3b6c", size = 320982, upload-time = "2025-10-02T22:51:39.077Z" },
+]
+
 [[package]]
 name = "pyparsing"
 version = "3.2.3"

From 6d2bd6c375bd828d3584214c52b5be99a2e041fd Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 14:52:29 -0800
Subject: [PATCH 17/20] Add _get_bathymetry() and draw filled contour of
 seafloor bottom on section plots.

---
 src/data/create_products.py | 151 +++++++++++++++++++++++++++++++++++-
 1 file changed, 148 insertions(+), 3 deletions(-)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index 8089bb9b..2fa281dc 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -22,6 +22,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
+import pooch
 import pyproj
 import xarray as xr
 
@@ -31,6 +32,41 @@
 from resample import AUVCTD_OPENDAP_BASE, FREQ
 from scipy.interpolate import griddata
 
+# Optional import for bathymetry data
+# Set GMT library path for macOS systems with MacPorts or Homebrew installations
+if sys.platform == "darwin":
+    # Try common installation paths, including MacPorts GMT6 subdirectory structure
+    gmt_search_paths = [
+        "/opt/local/lib/gmt6/lib",  # MacPorts GMT6
+        "/opt/local/lib",  # MacPorts older versions
+        "/opt/homebrew/lib",  # Homebrew ARM
+        "/usr/local/lib",  # Homebrew Intel
+    ]
+    for gmt_path in gmt_search_paths:
+        gmt_path_obj = Path(gmt_path)
+        if gmt_path_obj.exists() and any(gmt_path_obj.glob("libgmt.*")):
+            os.environ.setdefault("GMT_LIBRARY_PATH", gmt_path)
+            break
+
+try:
+    import pygmt
+
+    PYGMT_AVAILABLE = True
+
+    # Download Monterey Bay bathymetry grid using pooch for local caching
+    MONTEREY_BAY_GRID_URL = "https://stoqs.mbari.org/terrain/Monterey25.grd"
+    MONTEREY_BAY_GRID = pooch.retrieve(
+        url=MONTEREY_BAY_GRID_URL,
+        known_hash=None,  # Could add SHA256 hash for verification
+        fname="Monterey25.grd",
+        path=pooch.os_cache("auv-python"),
+        progressbar=False,
+    )
+except Exception:  # noqa: BLE001
+    # Any failure (ImportError, GMTCLibNotFoundError, etc.) - just continue without pygmt
+    PYGMT_AVAILABLE = False
+    MONTEREY_BAY_GRID = None
+
 # Define BASE_PATH for backward compatibility
 BASE_PATH = DEFAULT_BASE_PATH
 
@@ -274,6 +310,67 @@ def _grid_dims(self) -> tuple:
 
         return idist, iz, distnav
 
+    def _get_bathymetry(self, lons: np.ndarray, lats: np.ndarray) -> np.ndarray:
+        """Get bathymetry depths at given lon/lat positions using pygmt.
+
+        Args:
+            lons: Array of longitude values
+            lats: Array of latitude values
+
+        Returns:
+            Array of bathymetry depths (positive down) in meters, or None if pygmt unavailable
+        """
+        if not PYGMT_AVAILABLE:
+            self.logger.debug("pygmt not available, cannot retrieve bathymetry")
+            return None
+
+        # Use local Monterey Bay grid if available and coordinates are in range
+        # Otherwise fall back to global grids
+        points = pd.DataFrame({"lon": lons, "lat": lats})
+
+        # Check if coordinates are within Monterey Bay region
+        MB_LON_RANGE = (-122.5, -121.5)
+        MB_LAT_RANGE = (36.0, 37.5)
+        in_mb_region = (
+            (lons >= MB_LON_RANGE[0]).all()
+            and (lons <= MB_LON_RANGE[1]).all()
+            and (lats >= MB_LAT_RANGE[0]).all()
+            and (lats <= MB_LAT_RANGE[1]).all()
+        )
+
+        if in_mb_region and MONTEREY_BAY_GRID:
+            self.logger.info("Using local Monterey Bay bathymetry grid")
+            result = pygmt.grdtrack(
+                points=points,
+                grid=MONTEREY_BAY_GRID,
+                newcolname="depth",
+            )
+            # Convert to positive depths (meters below sea surface)
+            bathymetry = -result["depth"].to_numpy()
+            self.logger.info(
+                "Retrieved bathymetry data from Monterey Bay grid (min: %.1f m, max: %.1f m)",
+                bathymetry.min(),
+                bathymetry.max(),
+            )
+            return bathymetry
+
+        # Fall back to global grids
+        # Try GEBCO first (higher resolution), fall back to ETOPO1
+        result = pygmt.grdtrack(
+            points=points,
+            grid="@earth_relief_15s",  # 15 arc-second resolution (~450m)
+            newcolname="depth",
+        )
+
+        # Convert to positive depths (meters below sea surface)
+        bathymetry = -result["depth"].to_numpy()
+        self.logger.info(
+            "Retrieved bathymetry data using pygmt (min: %.1f m, max: %.1f m)",
+            bathymetry.min(),
+            bathymetry.max(),
+        )
+        return bathymetry
+
     def _profile_bottoms(
         self,
         distnav: xr.DataArray,
@@ -537,6 +634,7 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
         scale: str = "linear",
         num_colors: int = 256,
         gulper_locations: dict = None,
+        bottom_depths: np.array = None,
     ):
         # Handle both 1D and 2D axis arrays
         curr_ax = ax[row] if col == 0 and hasattr(ax, "ndim") and ax.ndim == 1 else ax[row, col]
@@ -699,6 +797,7 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
             extend="both",
             levels=np.linspace(v2_5, v97_5, num_colors),
         )
+
         # Blank out the countoured data below the bottom of the profiles
         xb = np.append(
             profile_bottoms.get_index("dist").to_numpy(),
@@ -720,6 +819,22 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
         )
         curr_ax.fill(list(reversed(xb)), list(reversed(yb)), "w")
 
+        # Add bathymetry as grey filled contour if available
+        if bottom_depths is not None:
+            # Pair bottom_depths with distance along track
+            dist_km = distnav.to_numpy() / 1000.0
+            # Create polygon for seafloor: from left to right along bottom,
+            # then back along axis bottom
+            xb_bathy = np.append(
+                dist_km,
+                [curr_ax.get_xlim()[1], curr_ax.get_xlim()[1], curr_ax.get_xlim()[0], dist_km[0]],
+            )
+            yb_bathy = np.append(
+                bottom_depths,
+                [bottom_depths[-1], curr_ax.get_ylim()[0], curr_ax.get_ylim()[0], bottom_depths[0]],
+            )
+            curr_ax.fill(xb_bathy, yb_bathy, color="#CCCCCC", zorder=1, alpha=0.8)
+
         # Add measurement location dots for density plot
         if var == "density":
             curr_ax.scatter(
@@ -827,6 +942,12 @@ def plot_2column(self) -> str:
         gulper_locations = self._get_gulper_locations(distnav)
 
         profile_bottoms = self._profile_bottoms(distnav)
+
+        bottom_depths = self._get_bathymetry(
+            self.ds.cf["longitude"].to_numpy(),
+            self.ds.cf["latitude"].to_numpy(),
+        )
+
         row = 1  # Start at row 1, col 0 (below the map)
         col = 0
         for var, scale in (
@@ -856,6 +977,7 @@ def plot_2column(self) -> str:
                 profile_bottoms,
                 scale=scale,
                 gulper_locations=gulper_locations,
+                bottom_depths=bottom_depths,
             )
             if row != 4:  # noqa: PLR2004
                 ax[row, col].get_xaxis().set_visible(False)
@@ -1021,6 +1143,27 @@ def gulper_odv(self, sec_bnds: int = 1) -> str:  # noqa: C901, PLR0912, PLR0915
             odv_column_names[24] = "bbp676 [m^{-1}]"
 
         best_ctd = self._get_best_ctd()
+
+        # Get bathymetry data for all gulper locations if available
+        bathymetry_dict = {}
+        if PYGMT_AVAILABLE:
+            gulper_lons = []
+            gulper_lats = []
+            for esec in gulper_times.values():
+                gulper_data = self.ds.sel(
+                    time=slice(
+                        np.datetime64(int((esec - sec_bnds) * 1e9), "ns"),
+                        np.datetime64(int((esec + sec_bnds) * 1e9), "ns"),
+                    ),
+                )
+                gulper_lons.append(gulper_data.cf["longitude"].to_numpy().mean())
+                gulper_lats.append(gulper_data.cf["latitude"].to_numpy().mean())
+
+            bathymetry = self._get_bathymetry(np.array(gulper_lons), np.array(gulper_lats))
+            if bathymetry is not None:
+                for i, bottle in enumerate(gulper_times.keys()):
+                    bathymetry_dict[bottle] = bathymetry[i]
+
         with gulper_odv_filename.open("w") as f:
             f.write("\t".join(odv_column_names) + "\n")
             for bottle, esec in gulper_times.items():
@@ -1056,9 +1199,11 @@ def gulper_odv(self, sec_bnds: int = 1) -> str:  # noqa: C901, PLR0912, PLR0915
                     elif name == "Lat (degrees_north)":
                         f.write(f"{gulper_data.cf['latitude'].to_numpy().mean():8.5f}")
                     elif name == "Bot. Depth [m]":
-                        f.write(
-                            f"{float(1000):8.1f}",
-                        )  # TODO: add proper bottom depth values
+                        # Use pygmt bathymetry if available, otherwise default to 1000m
+                        if bottle in bathymetry_dict:
+                            f.write(f"{bathymetry_dict[bottle]:8.1f}")
+                        else:
+                            f.write(f"{float(1000):8.1f}")
                     elif name == "Bottle Number [count]":
                         f.write(f"{bottle}")
                     elif name == "QF":

From 569571964b6c73fb1fd2a36b56955268e3c2566d Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 15:26:44 -0800
Subject: [PATCH 18/20] Update EXPECTED_MD5_GITHUB.

---
 src/data/test_process_dorado.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/data/test_process_dorado.py b/src/data/test_process_dorado.py
index 255bef4b..d1ff84a4 100644
--- a/src/data/test_process_dorado.py
+++ b/src/data/test_process_dorado.py
@@ -50,7 +50,7 @@ def test_process_dorado(complete_dorado_processing):
     check_md5 = True
     if check_md5:
         # Check that the MD5 hash has not changed
-        EXPECTED_MD5_GITHUB = "3fe4d2c1680db97289273905a386810c"
+        EXPECTED_MD5_GITHUB = "5dacd8536dbfb1ffe2f2c988d47a7257"
         EXPECTED_MD5_ACT = "1cd87ec363d3ca3323cf6df3f31ebe33"
         EXPECTED_MD5_LOCAL = "bd72eed42be79af3bff7c2657b092416"
         if str(proc.args.base_path).startswith("/home/runner"):

From 59cd7164284cdb0050660e13d32f739fbaf151d4 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 15:27:47 -0800
Subject: [PATCH 19/20] Add color bar (jet) for profile_number on maplplot, add
 {best_ctd} text over section plots.

---
 src/data/create_products.py | 56 +++++++++++++++++++++++++++++++------
 1 file changed, 47 insertions(+), 9 deletions(-)

diff --git a/src/data/create_products.py b/src/data/create_products.py
index 2fa281dc..2d1798d4 100755
--- a/src/data/create_products.py
+++ b/src/data/create_products.py
@@ -452,7 +452,7 @@ def _get_colormap_name(self, var: str) -> str:
         # Default
         return "cividis"
 
-    def _plot_track_map(
+    def _plot_track_map(  # noqa: PLR0915
         self, map_ax: matplotlib.axes.Axes, reference_ax: matplotlib.axes.Axes
     ) -> None:
         """Plot AUV track map on left side with title and times on right.
@@ -488,13 +488,13 @@ def _plot_track_map(
         # Make the plot square by using equal aspect (do this first)
         map_ax.set_aspect("equal", adjustable="datalim")
 
-        # Plot the track with depth coloring in Web Mercator coordinates
-        depths = self.ds.cf["depth"].to_numpy()
-        map_ax.scatter(
+        # Plot the track with profile_number coloring in Web Mercator coordinates
+        profile_numbers = self.ds["profile_number"].to_numpy()
+        scatter = map_ax.scatter(
             x_merc,
             y_merc,
-            c=depths,
-            cmap="viridis_r",  # Reversed so shallow is yellow, deep is dark
+            c=profile_numbers,
+            cmap="jet",
             s=1,
             alpha=0.6,
         )
@@ -527,6 +527,22 @@ def _plot_track_map(
 
         map_ax.set_position([ref_pos.x0, pos.y0, map_width, map_height])
 
+        # Add colorbar for profile numbers - create manually positioned axes
+        # to avoid affecting map position
+        # Position colorbar to the right of the map
+        cbar_width = 0.01
+        cbar_pad = 0.005
+        cbar_ax = map_ax.figure.add_axes(
+            [ref_pos.x0 + map_width + cbar_pad, pos.y0, cbar_width, map_height]
+        )
+        cbar = map_ax.figure.colorbar(
+            scatter,
+            cax=cbar_ax,
+            orientation="vertical",
+        )
+        cbar.set_label("Profile Number", fontsize=9)
+        cbar.ax.tick_params(labelsize=8)
+
         # Remove axes, labels, and ticks but keep the border
         map_ax.set_xticks([])
         map_ax.set_yticks([])
@@ -547,8 +563,9 @@ def _plot_track_map(
         # Get updated position after aspect adjustment
         updated_pos = map_ax.get_position()
 
-        # Position text in figure coordinates, to the right of the map
-        text_x = updated_pos.x0 + updated_pos.width + 0.01
+        # Position text in figure coordinates, to the right of the colorbar
+        # Account for colorbar width, padding, and space for colorbar label
+        text_x = ref_pos.x0 + map_width + cbar_pad + cbar_width + 0.03
         text_y = updated_pos.y0 + updated_pos.height  # Align with top of map
 
         # Add title
@@ -635,6 +652,7 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
         num_colors: int = 256,
         gulper_locations: dict = None,
         bottom_depths: np.array = None,
+        best_ctd: str = None,
     ):
         # Handle both 1D and 2D axis arrays
         curr_ax = ax[row] if col == 0 and hasattr(ax, "ndim") and ax.ndim == 1 else ax[row, col]
@@ -841,7 +859,7 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
                 distnav.to_numpy() / 1000.0,
                 self.ds.cf["depth"].to_numpy(),
                 s=0.1,
-                c="black",
+                c="white",
                 alpha=0.1,
                 zorder=3,
             )
@@ -917,6 +935,25 @@ def _plot_var(  # noqa: C901, PLR0912, PLR0913, PLR0915
         else:
             cb.set_label(long_name, fontsize=9)
 
+        # Add CTD label for density, temperature, and salinity plots
+        if best_ctd and (var == "density" or "_temperature" in var or "_salinity" in var):
+            # Position above the plot area (outside y-axis limits)
+            # Since y-axis is inverted (depth), min(iz) is at top, so go even less (shallower)
+            y_pos = (
+                min(iz) - (max(iz) - min(iz)) * 0.025
+            )  # 2.5% above the top (half the whitespace)
+            x_pos = curr_ax.get_xlim()[0]  # Left edge of plot
+            curr_ax.text(
+                x_pos,
+                y_pos,
+                best_ctd,
+                fontsize=8,
+                fontweight="bold",
+                verticalalignment="bottom",
+                horizontalalignment="left",
+                clip_on=False,
+            )
+
     def plot_2column(self) -> str:
         """Create 2column plot similar to plot_sections.m and stoqs/utils/Viz/plotting.py
         Construct a 2D grid of distance and depth and for each parameter grid the data
@@ -978,6 +1015,7 @@ def plot_2column(self) -> str:
                 scale=scale,
                 gulper_locations=gulper_locations,
                 bottom_depths=bottom_depths,
+                best_ctd=best_ctd,
             )
             if row != 4:  # noqa: PLR2004
                 ax[row, col].get_xaxis().set_visible(False)

From 6fc1ef5d9ffc4ca3c24cc69ad807078cf4c1f414 Mon Sep 17 00:00:00 2001
From: Mike McCann <MBARIMike@gmail.com>
Date: Mon, 15 Dec 2025 15:29:49 -0800
Subject: [PATCH 20/20] Go back to previous EXPECTED_MD5_GITHUB.

---
 src/data/test_process_dorado.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/data/test_process_dorado.py b/src/data/test_process_dorado.py
index d1ff84a4..255bef4b 100644
--- a/src/data/test_process_dorado.py
+++ b/src/data/test_process_dorado.py
@@ -50,7 +50,7 @@ def test_process_dorado(complete_dorado_processing):
     check_md5 = True
     if check_md5:
         # Check that the MD5 hash has not changed
-        EXPECTED_MD5_GITHUB = "5dacd8536dbfb1ffe2f2c988d47a7257"
+        EXPECTED_MD5_GITHUB = "3fe4d2c1680db97289273905a386810c"
         EXPECTED_MD5_ACT = "1cd87ec363d3ca3323cf6df3f31ebe33"
         EXPECTED_MD5_LOCAL = "bd72eed42be79af3bff7c2657b092416"
         if str(proc.args.base_path).startswith("/home/runner"):