diff --git a/echopype/convert/set_groups_base.py b/echopype/convert/set_groups_base.py index f0985cd58..7c028c817 100644 --- a/echopype/convert/set_groups_base.py +++ b/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 @@ -10,7 +11,11 @@ from ..utils.coding import COMPRESSION_SETTINGS, DEFAULT_TIME_ENCODING, set_time_encodings from ..utils.prov import echopype_prov_attrs, source_files_vars -NMEA_SENTENCE_DEFAULT = ["GGA", "GLL", "RMC"] +NMEA_SENTENCE_LOCATION = ["GGA", "GLL", "RMC"] +NMEA_SENTENCE_SPEED = ["RMC", "VTG"] +NMEA_SENTENCE_HEADING = ["HDT"] + +KNOTS_TO_M_PER_S = 0.51444444444 class SetGroupsBase(abc.ABC): @@ -182,10 +187,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 = [] @@ -199,6 +206,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_LOCATION) if nmea_msg: lat, lon = [], [] for x in nmea_msg: @@ -220,28 +264,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 @@ -467,9 +542,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, diff --git a/echopype/convert/set_groups_ek60.py b/echopype/convert/set_groups_ek60.py index bcb75047d..0dba497b2 100644 --- a/echopype/convert/set_groups_ek60.py +++ b/echopype/convert/set_groups_ek60.py @@ -187,6 +187,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": ""} @@ -195,8 +197,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( { @@ -250,6 +254,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": ( @@ -271,6 +285,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.", + }, + ), }, ) diff --git a/echopype/convert/set_groups_ek80.py b/echopype/convert/set_groups_ek80.py index d1a271389..34e4e5ec3 100644 --- a/echopype/convert/set_groups_ek80.py +++ b/echopype/convert/set_groups_ek80.py @@ -334,11 +334,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() @@ -356,6 +360,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( @@ -461,17 +485,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])), @@ -482,6 +496,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": ( @@ -520,6 +539,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) diff --git a/echopype/echodata/convention/1.0.yml b/echopype/echodata/convention/1.0.yml index 66068300b..ea7785105 100644 --- a/echopype/echodata/convention/1.0.yml +++ b/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 diff --git a/echopype/echodata/echodata.py b/echopype/echodata/echodata.py index 24ac2473f..57b6c3208 100644 --- a/echopype/echodata/echodata.py +++ b/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)} diff --git a/echopype/tests/convert/test_convert_ek60.py b/echopype/tests/convert/test_convert_ek60.py index 8bcd74798..79b710fe3 100644 --- a/echopype/tests/convert/test_convert_ek60.py +++ b/echopype/tests/convert/test_convert_ek60.py @@ -367,3 +367,32 @@ def test_open_raw_channels_invalid_ek60(ek60_path): sonar_model="EK60", channels=["nonexistent-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']))) diff --git a/echopype/tests/convert/test_convert_ek80.py b/echopype/tests/convert/test_convert_ek80.py index b508a68d7..421179f21 100644 --- a/echopype/tests/convert/test_convert_ek80.py +++ b/echopype/tests/convert/test_convert_ek80.py @@ -42,6 +42,11 @@ def ek80_new_path(test_path): def ek80_multiplex_path(test_path): return test_path["EK80_MULTIPLEX"] +@pytest.fixture +def ek80_heading_path(test_path): + return test_path["EK80_HEADING"] + + def pytest_generate_tests(metafunc): """Dynamically parameterize tests for EK80 .raw files.""" from echopype.tests import conftest as ct @@ -504,7 +509,7 @@ def test_parse_mru0_mru1(ek80_path): # Check dimensions assert ( echodata["Platform"].sizes - == {'channel': 1, 'time1': 1, 'time2': 43, 'time3': 43} + == {'channel': 1, 'time1': 1, 'time2': 43, 'time3': 43, 'time10': 1, 'time11': 1} ) # Check no NaN values in MRU data @@ -520,6 +525,35 @@ def test_parse_mru0_mru1(ek80_path): assert not np.any(np.isnan(echodata["Platform"][mru_var_name])) +@pytest.mark.unit +def test_parse_speed_over_ground(ek80_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=ek80_path/'vessel_speed'/'khr2405-D20241001-T024415.raw', + sonar_model='EK80' + ) + + # Check that there are data that are not NaN + assert (echodata["Platform"]['speed_over_ground'].sizes == {'time11': 220}) + assert (not np.any(np.isnan(echodata["Platform"]['speed_over_ground']))) + + +@pytest.mark.unit +def test_parse_NMEA_heading(ek80_heading_path): + """Make sure we parse NMEA heading from a RAW file when MRU heading is not present.""" + + echodata = open_raw( + raw_file=ek80_heading_path/'D20260613-T230914.raw', + sonar_model='ES80' + ) + + # Check that there are non-NaN data + assert (echodata["Platform"]['heading'].sizes == {'time10': 911}) + assert (not np.any(np.isnan(echodata["Platform"]['heading']))) + + @pytest.mark.unit def test_parse_missing_sound_velocity_profile(ek80_missing_sound_path): """ diff --git a/echopype/tests/echodata/test_echodata.py b/echopype/tests/echodata/test_echodata.py index 0bc2a9e87..32a8002fc 100644 --- a/echopype/tests/echodata/test_echodata.py +++ b/echopype/tests/echodata/test_echodata.py @@ -541,21 +541,21 @@ def test_update_platform( # times have max interval of 2s # check times are > min(ed["Sonar/Beam_group1"]["ping_time"]) - 2s assert ( - ed["Platform"]["time3"] + ed["Platform"]["time12"] > ed["Sonar/Beam_group1"]["ping_time"].min() - np.timedelta64(2, "s") ).all() # check there is only 1 time < min(ed["Sonar/Beam_group1"]["ping_time"]) assert ( - np.count_nonzero(ed["Platform"]["time3"] < ed["Sonar/Beam_group1"]["ping_time"].min()) <= 1 + np.count_nonzero(ed["Platform"]["time12"] < ed["Sonar/Beam_group1"]["ping_time"].min()) <= 1 ) # check times are < max(ed["Sonar/Beam_group1"]["ping_time"]) + 2s assert ( - ed["Platform"]["time3"] + ed["Platform"]["time12"] < ed["Sonar/Beam_group1"]["ping_time"].max() + np.timedelta64(2, "s") ).all() # check there is only 1 time > max(ed["Sonar/Beam_group1"]["ping_time"]) assert ( - np.count_nonzero(ed["Platform"]["time3"] > ed["Sonar/Beam_group1"]["ping_time"].max()) <= 1 + np.count_nonzero(ed["Platform"]["time12"] > ed["Sonar/Beam_group1"]["ping_time"].max()) <= 1 ) @@ -598,8 +598,8 @@ def test_update_platform_multidim(test_path): # Number of dimensions in Platform group and addition of time3 and time4 assert len(ed["Platform"].dims) == len(platform_preexisting_dims) + 2 - assert "time3" in ed["Platform"].dims - assert "time4" in ed["Platform"].dims + assert "time12" in ed["Platform"].dims + assert "time13" in ed["Platform"].dims # Dimension assignment assert ed["Platform"]["longitude"].dims[0] == ed["Platform"]["latitude"].dims[0] diff --git a/echopype/tests/echodata/test_echodata_combine.py b/echopype/tests/echodata/test_echodata_combine.py index c7e829412..c5af723e3 100644 --- a/echopype/tests/echodata/test_echodata_combine.py +++ b/echopype/tests/echodata/test_echodata_combine.py @@ -176,7 +176,17 @@ def test_combine_echodata(raw_datasets): eds = [echopype.open_raw(file, sonar_model, xml_file) for file in files] - append_dims = {"filenames", "time1", "time2", "time3", "nmea_time", "ping_time", "filter_time"} + append_dims = { + "filenames", + "time1", + "time2", + "time3", + "time10", + "time11", + "nmea_time", + "ping_time", + "filter_time", + } combined = echopype.combine_echodata(eds) @@ -227,7 +237,6 @@ def attr_time_to_dt(time_str): for dim in concat_dims: drop_dims = [c_dim for c_dim in concat_dims if c_dim != dim] - diff_concats.append(xr.concat([ed_subset.drop_dims(drop_dims) for ed_subset in eds_groups], dim=dim, # noqa: E501 coords="minimal", data_vars="minimal", join="outer"))