fix(commongrid): Fixed to include lat lon

Fixed 'compute_MVBS' function to now
include latitude and longitude if the variables
exists in the Sv dataset.

Additionally, the flox method and keyword arguments
are now exposed within the 'compute_MVBS' function.

Ref: Issue OSOceanAcoustics#1002

echopype/commongrid/api.py

@@ -7,6 +7,7 @@
 import pandas as pd
 import xarray as xr
 
+from ..consolidate.api import POSITION_VARIABLES
 from ..utils.prov import add_processing_level, echopype_prov_attrs, insert_input_processing_level
 from .mvbs import get_MVBS_along_channels
 from .nasc import (
@@ -71,7 +72,9 @@ def _set_MVBS_attrs(ds):
 
 
 @add_processing_level("L3*")
-def compute_MVBS(ds_Sv, range_meter_bin=20, ping_time_bin="20S"):
+def compute_MVBS(
+    ds_Sv, range_meter_bin=20, ping_time_bin="20S", method="map-reduce", **flox_kwargs
+):
     """
     Compute Mean Volume Backscattering Strength (MVBS)
     based on intervals of range (``echo_range``) and ``ping_time`` specified in physical units.
@@ -88,6 +91,13 @@ def compute_MVBS(ds_Sv, range_meter_bin=20, ping_time_bin="20S"):
         bin size along ``echo_range`` in meters, default to ``20``
     ping_time_bin : str
         bin size along ``ping_time``, default to ``20S``
+    method: str
+        The flox strategy for reduction of dask arrays only.
+        See flox `documentation <https://flox.readthedocs.io/en/latest/implementation.html>`_
+        for more details.
+    **kwargs
+        Additional keyword arguments to be passed
+        to flox reduction function.
 
     Returns
     -------
@@ -112,19 +122,26 @@ def compute_MVBS(ds_Sv, range_meter_bin=20, ping_time_bin="20S"):
     ping_interval = d_index.union([d_index[-1] + pd.Timedelta(ping_time_bin)])
 
     # calculate the MVBS along each channel
-    raw_MVBS = get_MVBS_along_channels(ds_Sv, range_interval, ping_interval)
+    raw_MVBS = get_MVBS_along_channels(
+        ds_Sv, range_interval, ping_interval, method=method, **flox_kwargs
+    )
 
     # create MVBS dataset
     # by transforming the binned dimensions to regular coords
     ds_MVBS = xr.Dataset(
-        data_vars={"Sv": (["channel", "ping_time", "echo_range"], raw_MVBS.data)},
+        data_vars={"Sv": (["channel", "ping_time", "echo_range"], raw_MVBS["Sv"].data)},
         coords={
             "ping_time": np.array([v.left for v in raw_MVBS.ping_time_bins.values]),
             "channel": raw_MVBS.channel.values,
             "echo_range": np.array([v.left for v in raw_MVBS.echo_range_bins.values]),
         },
     )
 
+    # Add the position variables
+    if raw_MVBS.attrs.get("has_positions", False):
+        for var in POSITION_VARIABLES:
+            ds_MVBS[var] = (["ping_time"], raw_MVBS[var].data)
+
     # ping_time_bin parsing and conversions
     # Need to convert between pd.Timedelta and np.timedelta64 offsets/frequency strings
     # https://xarray.pydata.org/en/stable/generated/xarray.Dataset.resample.html

echopype/commongrid/mvbs.py

@@ -10,6 +10,8 @@
 import xarray as xr
 from flox.xarray import xarray_reduce
 
+from ..consolidate.api import POSITION_VARIABLES
+
 
 def get_bin_indices(
     echo_range: np.ndarray, bins_er: np.ndarray, times: np.ndarray, bins_time: np.ndarray
@@ -467,6 +469,21 @@ def get_MVBS_along_channels(
     # average should be done in linear domain
     sv = ds_Sv["Sv"].pipe(_linear_transform)
 
+    # Get positions if exists
+    # otherwise just use an empty dataset
+    ds_Pos = xr.Dataset(attrs={"has_positions": False})
+    if all(v in ds_Sv for v in POSITION_VARIABLES):
+        ds_Pos = xarray_reduce(
+            ds_Sv[POSITION_VARIABLES],
+            ds_Sv["ping_time"],
+            func="nanmean",
+            expected_groups=(ping_interval),
+            isbin=True,
+            engine="flox",
+            method=method,
+        )
+        ds_Pos.attrs["has_positions"] = True
+
     # reduce along ping_time and echo_range
     # by binning and averaging
     mvbs = xarray_reduce(
@@ -483,4 +500,5 @@ def get_MVBS_along_channels(
     )
 
     # apply inverse mapping to get back to the original domain and store values
-    return mvbs.pipe(_linear_transform, inverse=True)
+    da_MVBS = mvbs.pipe(_linear_transform, inverse=True)
+    return xr.merge([ds_Pos, da_MVBS])

echopype/consolidate/api.py

@@ -12,6 +12,8 @@
 from ..utils.prov import add_processing_level
 from .split_beam_angle import add_angle_to_ds, get_angle_complex_samples, get_angle_power_samples
 
+POSITION_VARIABLES = ["latitude", "longitude"]
+
 
 def swap_dims_channel_frequency(ds: xr.Dataset) -> xr.Dataset:
     """
@@ -185,7 +187,7 @@ def sel_interp(var, time_dim_name):
         f"{datetime.datetime.utcnow()} +00:00. "
         "Interpolated or propagated from Platform latitude/longitude."  # noqa
     )
-    for da_name in ["latitude", "longitude"]:
+    for da_name in POSITION_VARIABLES:
         interp_ds[da_name] = interp_ds[da_name].assign_attrs({"history": history_attr})
 
     if time_dim_name in interp_ds: