Support speed and NMEA heading from Simrad .raw files

echopype/convert/set_groups_base.py

@@ -1,5 +1,6 @@
 import abc
 import warnings
+from collections.abc import Iterable
 from typing import List, Set
 
 import numpy as np
@@ -11,6 +12,10 @@
 from ..utils.prov import echopype_prov_attrs, source_files_vars
 
 NMEA_SENTENCE_DEFAULT = ["GGA", "GLL", "RMC"]
+NMEA_SENTENCE_SPEED = ["RMC", "VTG"]
+NMEA_SENTENCE_HEADING = ["HDT"]
+
+KNOTS_TO_M_PER_S = 0.51444444444
 
 
 class SetGroupsBase(abc.ABC):
@@ -177,10 +182,12 @@ def set_vendor(self) -> xr.Dataset:
         raise NotImplementedError
 
     # TODO: move this to be part of parser as it is not a "set" operation
-    def _extract_NMEA_latlon(self):
-        """Get the lat and lon values from the raw nmea data"""
+    def _extract_selected_NMEA(
+        self, nmea_sentence_types: list
+    ) -> tuple[Iterable, Iterable, Iterable]:
+        """Parse out the selected NMEA messages."""
         messages = [string[3:6] for string in self.parser_obj.nmea["nmea_string"]]
-        idx_loc = np.argwhere(np.isin(messages, NMEA_SENTENCE_DEFAULT)).squeeze()
+        idx_loc = np.argwhere(np.isin(messages, nmea_sentence_types)).squeeze()
         if idx_loc.size == 1:  # in case of only 1 matching message
             idx_loc = np.expand_dims(idx_loc, axis=0)
         nmea_msg = []
@@ -194,6 +201,43 @@ def _extract_NMEA_latlon(self):
                 pynmea2.ParseError,
             ):
                 nmea_msg.append(None)
+
+        msg_type = (
+            [x.sentence_type if hasattr(x, "sentence_type") else np.nan for x in nmea_msg]
+            if nmea_msg
+            else [np.nan]
+        )
+
+        if nmea_msg:
+            time, _, _ = xr.coding.times.encode_cf_datetime(
+                np.array(self.parser_obj.nmea["timestamp"])[idx_loc],
+                **{
+                    "units": DEFAULT_TIME_ENCODING["units"],
+                    "calendar": DEFAULT_TIME_ENCODING["calendar"],
+                },
+            )
+            time = xr.coding.times.decode_cf_datetime(
+                time,
+                units=DEFAULT_TIME_ENCODING["units"],
+                calendar=DEFAULT_TIME_ENCODING["calendar"],
+            )
+        else:
+            time = [np.nan]
+
+        # There can be duplicate timestamps both due to a problem in earlier Simrad raw
+        # files and if multiple NMEA sentences are used with the same
+        # timestamp. Remove them here.
+        if nmea_msg:
+            time, indices = np.unique(time, return_index=True, sorted=False)
+            nmea_msg = list(np.array(nmea_msg)[indices])
+            msg_type = list(np.array(msg_type)[indices])
+
+        return nmea_msg, time, msg_type
+
+    # TODO: move this to be part of parser as it is not a "set" operation
+    def _extract_NMEA_latlon(self):
+        """Get the lat and lon values from the raw nmea data"""
+        nmea_msg, time, msg_type = self._extract_selected_NMEA(NMEA_SENTENCE_DEFAULT)
         if nmea_msg:
             lat, lon = [], []
             for x in nmea_msg:
@@ -215,28 +259,59 @@ def _extract_NMEA_latlon(self):
                     )
         else:
             lat, lon = [np.nan], [np.nan]
-        msg_type = (
-            [x.sentence_type if hasattr(x, "sentence_type") else np.nan for x in nmea_msg]
-            if nmea_msg
-            else [np.nan]
-        )
+
+        return time, msg_type, lat, lon
+
+    # TODO: move this to be part of parser as it is not a "set" operation
+    def _extract_NMEA_speed(self):
+        """Get the speed over ground values from the raw nmea data"""
+        nmea_msg, time, msg_type = self._extract_selected_NMEA(NMEA_SENTENCE_SPEED)
         if nmea_msg:
-            time1, _, _ = xr.coding.times.encode_cf_datetime(
-                np.array(self.parser_obj.nmea["timestamp"])[idx_loc],
-                **{
-                    "units": DEFAULT_TIME_ENCODING["units"],
-                    "calendar": DEFAULT_TIME_ENCODING["calendar"],
-                },
-            )
-            time1 = xr.coding.times.decode_cf_datetime(
-                time1,
-                units=DEFAULT_TIME_ENCODING["units"],
-                calendar=DEFAULT_TIME_ENCODING["calendar"],
-            )
+            sog = []
+            for x in nmea_msg:
+                try:
+                    # pynmea2 has different names for speed over ground, depending on the NMEA
+                    # message that it comes from
+                    if x.sentence_type == "VTG":
+                        # VTG speed is returned as a Decimal, so fix that
+                        sog.append(
+                            float(x.spd_over_grnd_kts) * KNOTS_TO_M_PER_S
+                            if hasattr(x, "spd_over_grnd_kts") and x.spd_over_grnd_kts is not None
+                            else np.nan
+                        )
+                    else:  # only RMC so far
+                        sog.append(
+                            x.spd_over_grnd * KNOTS_TO_M_PER_S
+                            if hasattr(x, "spd_over_grnd") and x.spd_over_grnd is not None
+                            else np.nan
+                        )
+                except ValueError:
+                    sog.append(np.nan)
+        else:
+            sog = [np.nan]
+
+        return time, msg_type, sog
+
+    # TODO: move this to be part of parser as it is not a "set" operation
+    def _extract_NMEA_heading(self):
+        """Get heading values from the raw nmea data"""
+        nmea_msg, time, msg_type = self._extract_selected_NMEA(NMEA_SENTENCE_HEADING)
+        if nmea_msg:
+            heading = []
+            for x in nmea_msg:
+                try:
+                    # HDG speed is returned as a Decimal, so fix that
+                    heading.append(
+                        float(x.heading)
+                        if hasattr(x, "heading") and x.heading is not None
+                        else np.nan
+                    )
+                except ValueError:
+                    heading.append(np.nan)
         else:
-            time1 = [np.nan]
+            heading = [np.nan]
 
-        return time1, msg_type, lat, lon
+        return time, msg_type, heading
 
     def _beam_groups_vars(self):
         """Stage beam_group coordinate and beam_group_descr variables sharing
@@ -462,9 +537,7 @@ def _add_index_data_to_platform_ds(
             }
         )
 
-        return platform_ds.transpose(
-            "channel", "time1", "time2", "time3", "time4", missing_dims="ignore"
-        )
+        return platform_ds.transpose(missing_dims="ignore")
 
     def _add_seafloor_detection_data_to_vendor_ds(
         self,

echopype/convert/set_groups_ek60.py

@@ -180,6 +180,8 @@ def set_platform(self) -> xr.Dataset:
         # Collect variables
         # Read lat/long from NMEA datagram
         time1, msg_type, lat, lon = self._extract_NMEA_latlon()
+        time10, msg_type_heading, heading = self._extract_NMEA_heading()
+        time11, msg_type_sog, sog = self._extract_NMEA_speed()
 
         # NMEA dataset: variables filled with np.nan if they do not exist
         platform_dict = {"platform_name": "", "platform_type": "", "platform_code_ICES": ""}
@@ -188,8 +190,10 @@ def set_platform(self) -> xr.Dataset:
         # are identical across channels
         ch = list(self.sorted_channel.keys())[0]
 
-        # Handle potential nan timestamp for time1 and time2
+        # Handle potential nan timestamp
         time1 = self._nan_timestamp_handler(time1)
+        time10 = self._nan_timestamp_handler(time10)
+        time11 = self._nan_timestamp_handler(time11)
 
         ds = xr.Dataset(
             {
@@ -243,6 +247,16 @@ def set_platform(self) -> xr.Dataset:
                         "position_offset_z",
                     ]
                 },
+                "heading": (
+                    ["time10"],
+                    np.array(heading),
+                    self._varattrs["platform_var_default"]["heading"],
+                ),
+                "speed_over_ground": (
+                    ["time11"],
+                    np.array(sog),
+                    self._varattrs["platform_var_default"]["speed_over_ground"],
+                ),
             },
             coords={
                 "time1": (
@@ -264,6 +278,22 @@ def set_platform(self) -> xr.Dataset:
                         "orientation data.",
                     },
                 ),
+                "time10": (
+                    ["time10"],
+                    time10,
+                    {
+                        **self._varattrs["platform_coord_default"]["time1"],
+                        "comment": "Time coordinate corresponding to NMEA heading data.",
+                    },
+                ),
+                "time11": (
+                    ["time11"],
+                    time11,
+                    {
+                        **self._varattrs["platform_coord_default"]["time1"],
+                        "comment": "Time coordinate corresponding to NMEA speed data.",
+                    },
+                ),
             },
         )
 

echopype/convert/set_groups_ek80.py

@@ -329,11 +329,15 @@ def set_platform(self) -> xr.Dataset:
         time2 = np.array(time2) if time2 is not None else [np.nan]
         time3 = self.parser_obj.mru1.get("timestamp", None)
         time3 = np.array(time3) if time3 is not None else [np.nan]
+        time10, msg_type_heading, heading_nmea = self._extract_NMEA_heading()
+        time11, msg_type_sog, sog_nmea = self._extract_NMEA_speed()
 
-        # Handle potential nan timestamp for time1, time2, and time3
+        # Handle potential nan timestamps
         time1 = self._nan_timestamp_handler(time1)
         time2 = self._nan_timestamp_handler(time2)
         time3 = self._nan_timestamp_handler(time3)
+        time10 = self._nan_timestamp_handler(time10)
+        time11 = self._nan_timestamp_handler(time11)
 
         # Set MRU1 lat lon attributes
         latitude_mru1_attrs = self._varattrs["platform_var_default"]["latitude"].copy()
@@ -351,6 +355,26 @@ def set_platform(self) -> xr.Dataset:
             }
         ),
 
+        # If there is no heading data from an MRU but there is from the NMEA data, use that instead
+        if "heading" in self.parser_obj.mru0:
+            hdg_data = (
+                ["time2"],
+                np.array(self.parser_obj.mru0.get("heading", [np.nan])),
+                self._varattrs["platform_var_default"]["heading"],
+            )
+        elif len(heading_nmea) > 0:
+            hdg_data = (
+                ["time10"],
+                heading_nmea,
+                self._varattrs["platform_var_default"]["heading"],
+            )
+        else:
+            hdg_data = (
+                ["time2"],
+                [np.nan],
+                self._varattrs["platform_var_default"]["heading"],
+            )
+
         # Assemble variables into a dataset: variables filled with nan if do not exist
         platform_dict = {"platform_name": "", "platform_type": "", "platform_code_ICES": ""}
         ds = xr.Dataset(
@@ -456,17 +480,7 @@ def set_platform(self) -> xr.Dataset:
                         "standard_name": "sound_frequency",
                     },
                 ),
-                "heading": (
-                    ["time2"],
-                    np.array(self.parser_obj.mru0.get("heading", [np.nan])),
-                    {
-                        "long_name": "Platform heading (true)",
-                        "standard_name": "platform_orientation",
-                        "units": "degrees_north",
-                        "valid_min": 0.0,
-                        "valid_max": 360.0,
-                    },
-                ),
+                "heading": hdg_data,
                 "latitude_mru1": (
                     ["time3"],
                     np.array(self.parser_obj.mru1.get("latitude", [np.nan])),
@@ -477,6 +491,11 @@ def set_platform(self) -> xr.Dataset:
                     np.array(self.parser_obj.mru1.get("longitude", [np.nan])),
                     longitude_mru1_attrs,
                 ),
+                "speed_over_ground": (
+                    ["time11"],
+                    np.array(sog_nmea),
+                    self._varattrs["platform_var_default"]["speed_over_ground"],
+                ),
             },
             coords={
                 "channel": (
@@ -515,6 +534,22 @@ def set_platform(self) -> xr.Dataset:
                         "orientation data from the Kongsberg Maritime Binary Datagram.",
                     },
                 ),
+                "time10": (
+                    ["time10"],
+                    time10,
+                    {
+                        **self._varattrs["platform_coord_default"]["time1"],
+                        "comment": "Time coordinate corresponding to NMEA heading data.",
+                    },
+                ),
+                "time11": (
+                    ["time11"],
+                    time11,
+                    {
+                        **self._varattrs["platform_coord_default"]["time1"],
+                        "comment": "Time coordinate corresponding to NMEA speed data.",
+                    },
+                ),
             },
         )
         ds = ds.assign_attrs(platform_dict)

echopype/echodata/convention/1.0.yml

@@ -156,5 +156,16 @@ variable_and_varattributes:
     water_level:
       long_name: Distance from the platform coordinate system origin to the nominal water level along the z-axis
       units: m
+    speed_over_ground:
+      long_name: Speed over ground of the platform
+      standard_name: platform_speed_wrt_ground
+      units: m/s
+      valid_min: 0.0
+    heading:
+      long_name: Platform heading (true)
+      standard_name: platform_orientation
+      units: degrees_north
+      valid_min: 0.0
+      valid_max: 360.0
     sentence_type:
       long_name: NMEA sentence type

echopype/echodata/echodata.py

@@ -439,7 +439,7 @@ def update_platform(
                 if v["ext_time_dim_name"] != "scalar"
             }
         )
-        time_dims_max = max([int(dim[-1]) for dim in platform.dims if dim.startswith("time")])
+        time_dims_max = max([int(dim[4:]) for dim in platform.dims if dim.startswith("time")])
         new_time_dims = [f"time{time_dims_max + i + 1}" for i in range(len(ext_time_dims))]
         # Map each new time dim name to the external time dim name:
         new_time_dims_mappings = {new: ext for new, ext in zip(new_time_dims, ext_time_dims)}

echopype/tests/convert/test_convert_ek60.py

@@ -264,3 +264,33 @@ def test_converting_ek60_raw_with_missing_channel_power(ek60_missing_channel_pow
     # Check that all empty power channels do not exist in the EchoData Beam group
     for _, empty_power_channel_name in empty_power_chs.items():
         assert empty_power_channel_name not in ed["Sonar/Beam_group1"]["channel"]
+
+
+@pytest.mark.unit
+def test_parse_speed_over_ground(ek60_path):
+    """Make sure we parse speed over ground from a RAW file."""
+
+    # This raw file has speed in NMEA VTG and RMC messages
+    echodata = open_raw(
+        raw_file=ek60_path/'NBP_B050N-D20180118-T090228.raw',
+        sonar_model='EK60'
+    )
+
+    # Check that there are data that are not NaN
+    assert (echodata["Platform"]['speed_over_ground'].sizes == {'time11': 584})
+    # this .raw file has nan's in the speed over ground data 
+    # assert (not np.any(np.isnan(echodata["Platform"]['speed_over_ground'])))
+
+
+@pytest.mark.unit
+def test_parse_NMEA_heading(ek60_path):
+    """Make sure we parse NMEA heading from a RAW file when MRU heading is not present."""
+
+    echodata = open_raw(
+        raw_file=ek60_path/'NBP_B050N-D20180118-T090228.raw',
+        sonar_model='EK60'
+    )
+
+    # Check that there are non-NaN data
+    assert (echodata["Platform"]['heading'].sizes == {'time10': 584})
+    assert (not np.any(np.isnan(echodata["Platform"]['heading'])))