correct Dorado HS2 bbp calculation

[messiem]

bbp should be particulate backscatter (unit m-1). The method used to remove water backscatter was corrected in this commit but the new calculation needs to be double-checked against coincident Ecopuck data.

[MBARIMike]

As a first step to correct the Dorado HS2 bbp calculation here's a comparison of the correction "proposed" in this commit additional variables with the "_fixed" suffix added. Here are the (temporary) changes made to the code:

(auv-python-py3.10) ➜  auv-python git:(main) ✗ git diff f419b0d 51d4862 src/data/hs2_proc.py
diff --git a/src/data/hs2_proc.py b/src/data/hs2_proc.py
index 00fd28d..456784d 100644
--- a/src/data/hs2_proc.py
+++ b/src/data/hs2_proc.py
@@ -159,13 +159,13 @@ def hs2_calc_bb(orig_nc, cals):

         b_b_corr = sigma * b_b_uncorr
         # Need to test subtracting b_bw here instead of after multiplying by sigma
-        # b_bp_corr = sigma * (b_b_uncorr - b_bw)
+        b_bp_corr = sigma * (b_b_uncorr - b_bw)

         setattr(hs2, f"bb{wavelength}", b_b_corr)
         # Legacy code that subtracts b_bw after multiplying by sigma
         setattr(hs2, f"bbp{wavelength}", b_b_corr - b_bw)
         # This is likely the correct way to do it, with b_bw subtracted before multiplying by sigma
-        # setattr(hs2, f"bbp{wavelength}", b_bp_corr)
+        setattr(hs2, f"bbp{wavelength}_fixed", b_bp_corr)

     # -% 'hs2.fl700_uncorr = (hs2.Snorm3.*50)./((1 + str2num(CAL.Ch(3).TempCoeff).*(hs2.Temp-str2num(CAL.General.CalTemp))).*hs2.Gain3.*str2num(CAL.Ch(3).RNominal));'  # noqa: E501
     denom = (

(auv-python-py3.10) ➜  auv-python git:(main) ✗ git diff 2b2462c 51d4862 src/data/calibrate.py
diff --git a/src/data/calibrate.py b/src/data/calibrate.py
index a1e8648..328bdaf 100755
--- a/src/data/calibrate.py
+++ b/src/data/calibrate.py
@@ -2271,6 +2271,19 @@ class Calibrate_NetCDF:
                 "units": "m-1",
                 "comment": (f"Computed by hs2_calc_bb() from data in {source}"),
             }
+        if hasattr(hs2, "bbp420_fixed"):
+            blue_bs_fixed = xr.DataArray(
+                hs2.bbp420_fixed.to_numpy(),
+                coords=[hs2.bbp420_fixed.get_index("time")],
+                dims={"hs2_time"},
+                name="hs2_bbp420_fixed",
+            )
+            blue_bs_fixed.attrs = {
+                "long_name": "Particulate backscattering coefficient at 420 nm",
+                "coordinates": coord_str,
+                "units": "m-1",
+                "comment": (f"Computed by hs2_calc_bb() from data in {source}"),
+            }
         if hasattr(hs2, "bbp470"):
             blue_bs = xr.DataArray(
                 hs2.bbp470.to_numpy(),
@@ -2312,6 +2325,19 @@ class Calibrate_NetCDF:
                 "units": "m-1",
                 "comment": (f"Computed by hs2_calc_bb() from data in {source}"),
             }
+        if hasattr(hs2, "bbp700_fixed"):
+            red_bs_fixed = xr.DataArray(
+                hs2.bbp700_fixed.to_numpy(),
+                coords=[hs2.bbp700_fixed.get_index("time")],
+                dims={"hs2_time"},
+                name="hs2_bbp700_fixed",
+            )
+            red_bs_fixed.attrs = {
+                "long_name": "Particulate backscattering coefficient at 700 nm",
+                "coordinates": coord_str,
+                "units": "m-1",
+                "comment": (f"Computed by hs2_calc_bb() from data in {source}"),
+            }

         # Fluorescence
         if hasattr(hs2, "fl676"):
@@ -2384,13 +2410,15 @@ class Calibrate_NetCDF:
                 pend = int(self.args.plot.split("first")[1])
             df_plot = pd.DataFrame(index=blue_bs.get_index("hs2_time")[pbeg:pend])
             df_plot["blue_bs"] = blue_bs[pbeg:pend]
+            df_plot["blue_bs_fixed"] = blue_bs_fixed[pbeg:pend]
             df_plot["red_bs"] = red_bs[pbeg:pend]
-            df_plot["fl"] = fl[pbeg:pend]
+            df_plot["red_bs_fixed"] = red_bs_fixed[pbeg:pend]
+            ## df_plot["fl"] = fl[pbeg:pend]
             title = (
                 f"First {pend} points from"
                 f" {self.args.mission}/{self.sinfo[sensor]['data_filename']}"
             )
-            ax = df_plot.plot(title=title, figsize=(18, 6))
+            ax = df_plot.plot(title=title, figsize=(18, 6), ylim=(-0.003, 0.004))
             ax.grid("on")
             self.logger.debug("Pausing with plot entitled: %s. Close window to continue.", title)
             plt.show()
@@ -2398,12 +2426,16 @@ class Calibrate_NetCDF:
         # Save blue, red, & fl to combined_nc, alsoe
         if hasattr(hs2, "bbp420"):
             self.combined_nc["hs2_bbp420"] = blue_bs
+        if hasattr(hs2, "bbp420_fixed"):
+            self.combined_nc["hs2_bbp420_fixed"] = blue_bs_fixed
         if hasattr(hs2, "bbp470"):
             self.combined_nc["hs2_bbp470"] = blue_bs
         if hasattr(hs2, "bbp676"):
             self.combined_nc["hs2_bbp676"] = red_bs
         if hasattr(hs2, "bbp700"):
             self.combined_nc["hs2_bbp700"] = red_bs
+        if hasattr(hs2, "bbp700_fixed"):
+            self.combined_nc["hs2_bbp700_fixed"] = red_bs_fixed
         if hasattr(hs2, "fl676"):
             self.combined_nc["hs2_fl676"] = fl
         if hasattr(hs2, "fl700"):

And here is a plot of 10 seconds of data comparing the original and fixed values of bbp:

This plot was made by executing src/data/calibrate.py with the --plot option. This line in .vscode/launch.json:

"args": ["--auv_name", "dorado", "--mission", "2023.324.00", "-v", "1", "--plot", "first10000"]

The differences are tiny.

[MBARIMike]

Examining these data more closely at a breakpoint in the hs2_calc_bb() function...

Typical (mean) values [channel 1 (blue) from 2023.324.00]:

print(f"{float(sigma.mean()) =:9.6f}")
float(sigma.mean()) = 1.004638

print(f"{float(b_b_corr.mean()) =:9.6f}")
float(b_b_corr.mean()) = 0.003301

print(f"{b_bw =:9.6f}")
b_bw = 0.003068

In comparing the calculations sigma * (b_b_uncorr - b_bw) vs. sigma * b_b_uncorr - b_bw we see that the difference is about 0.4% of b_bw, which seems to explain the tiny differences in the plot above. We had hoped that the fixed values would correct the negative values we often see in the backscatter data, but the values still go negative:

print(f"{float((b_b_corr - b_bw).min()) =:9.6f}")
float((b_b_corr - b_bw).min()) =-0.007809

print(f"{float(b_bp_corr.min()) =:9.6f}")
float(b_bp_corr.min()) =-0.007718

[MBARIMike]

Now, let's compare against coincident Ecopuck data...

First let's use STOQS to give us a linear regression of hs2_bb700 (m-1) vs. ecopuck_bbp700 (m^-1 sr^-1) from 2023.324.00 in stoqs_all_dorado: http://stoqs.mbari.org/p/DkUHQV4

There are many ecopuck_bbp700 (m^-1 sr^-1) outliers with values above 0.0003 that skew the regression. Removing these gives about a factor of about 6 between the two measures of "red" backscatter:

But, we're not sure if we are comparing apples to apples given that the units of these values don't match. It's possible that the ecopuck backscatter is volumetric backscatter vs. a backscatter coefficient.

[MBARIMike]

Rereading https://misclab.umeoce.maine.edu/OceanOpticsClass2021/wp-content/uploads/2021/07/Lecture_08_VSF_ScatModels_CM.pdf I see that there is no simple conversion factor for volumetric backscatter to backscatter coefficient. It's dependent on the medium and the nature of the scattering. This paper also discusses the WET Labs ECObb (ECOvsf) and HOBI Labs Hydroscat instruments: https://misclab.umeoce.maine.edu/ftp/classes/OO2015/Readings/Roesler_Boss_final_Ch5.pdf

Examining a ~20 minute period it seems that the backscatter from these two instruments co-vary reasonably:

[MBARIMike]

A simplistic look a the hs2_calc_bb() code shows this line, which applies a 6.7858 factor conversion from "beta" values to "b" values. Are the ecopuck values beta and the hs2 values b? (I feel out of my depth here with all of this optical oceanography.)

[MBARIMike]

The hs2_calc_bb() code appears to make some of the same calculations that _compute_backscatter() does, with some frustrating to understand differences. To simplify things lets try using _compute_backscatter() for computation of the particle backscatter coefficient from the hs2 beta_uncorr variable.

[MBARIMike]

With #10 the HS2 bbp calculation is simplified and more consistent with how ecopuck bbp is calculated. Here is a comparison of HS2 bbp700 and ecopuck beta700 (the plot label "ecopuck_bbp700" is incorrect, it's actually beta) from before this code change:

and after this code change:

It's somewhat comforting that the relationships appear similar indicating that mainly linear adjustments were made by the code change. The R-squared for the latter increases slightly, from 0.957 to 0.958. The units and data values between HS2 and ecopuck bbp are now in much better agreement with this regression:

OLS linear regression: hs2_bbp700 (m-1) = 0.8256 * ecopuck_bbp700 (m-1) + 0.0001897 (R-squared = 0.958)

[MBARIMike]

Going back to the original Sackmann code improves the slope a bit. It's now:

OLS linear regression: hs2_bbp700 (m-1) = 0.8539 * ecopuck_bbp700 (m-1) + 2.086e-05 (R-squared = 0.957)

[MBARIMike]

After going through the math to compute backscatter from HS2ManualRevI-2011-12.pdf @messiem found an error that's corrected in #15.