From 5cf4daf8fc20e53bc547706ce50f7407f42ae955 Mon Sep 17 00:00:00 2001 From: ctuguinay Date: Thu, 25 Jul 2024 19:39:15 +0000 Subject: [PATCH] add align_to_ping_time function, and replace the other interp functions --- echopype/calibrate/env_params.py | 5 +- echopype/consolidate/ek_depth_utils.py | 30 ++---------- echopype/consolidate/loc_utils.py | 3 +- echopype/tests/calibrate/test_env_params.py | 11 +---- .../tests/consolidate/test_add_location.py | 6 ++- echopype/tests/utils/test_align.py | 48 +++++++++++++++++++ echopype/utils/align.py | 42 ++++++++++++++++ 7 files changed, 103 insertions(+), 42 deletions(-) create mode 100644 echopype/tests/utils/test_align.py create mode 100644 echopype/utils/align.py diff --git a/echopype/calibrate/env_params.py b/echopype/calibrate/env_params.py index 811b3476f..ab2b8499f 100644 --- a/echopype/calibrate/env_params.py +++ b/echopype/calibrate/env_params.py @@ -6,6 +6,7 @@ from ..echodata import EchoData from ..utils import uwa +from ..utils.align import align_to_ping_time from .cal_params import param2da ENV_PARAMS = ( @@ -68,9 +69,7 @@ def harmonize_env_param_time( return p.rename({"time1": "ping_time"}) # Interpolate `p` to `ping_time` - return ( - p.dropna(dim="time1").interp(time1=ping_time).ffill(dim="ping_time").drop_vars("time1") - ) + return align_to_ping_time(p.dropna(dim="time1"), "time1", ping_time, method="linear") return p diff --git a/echopype/consolidate/ek_depth_utils.py b/echopype/consolidate/ek_depth_utils.py index 6a83f4d09..4b6e9084e 100644 --- a/echopype/consolidate/ek_depth_utils.py +++ b/echopype/consolidate/ek_depth_utils.py @@ -2,6 +2,7 @@ import xarray as xr from scipy.spatial.transform import Rotation as R +from ..utils.align import align_to_ping_time from ..utils.log import _init_logger logger = _init_logger(__name__) @@ -26,31 +27,6 @@ def _check_and_log_nans( ) -def _var_time2_to_ping_time(var_with_time2, ping_time_da): - """ - If `time2` does not differ from `var`, we rename `time2` to 'ping_time', - else interpolate `transducer_depth`'s `time2` dimension to `ping_time`. - - Useful for handling EK60 and EK80 platform data: - - EK80 will have platform variables with time2 dimension that does not match - Beam group ping time values, while EK60 will have time2 dimension that - matches Beam group ping time values. - """ - if not ping_time_da.equals(var_with_time2["time2"].rename({"time2": "ping_time"})): - var_with_ping_time = var_with_time2.interp( - {"time2": ping_time_da}, - method="nearest", - # More details for `fill_value` and `extrapolate` can be found here: - # https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html # noqa - kwargs={"fill_value": "extrapolate"}, - ).drop_vars("time2") - else: - var_with_ping_time = var_with_time2.rename({"time2": "ping_time"}) - - return var_with_ping_time - - def ek_use_platform_vertical_offsets( platform_ds: xr.Dataset, ping_time_da: xr.DataArray, @@ -72,7 +48,7 @@ def ek_use_platform_vertical_offsets( # Compute z translation for transducer position vs water level transducer_depth = transducer_offset_z - (water_level + vertical_offset) - return _var_time2_to_ping_time(transducer_depth, ping_time_da) + return align_to_ping_time(transducer_depth, "time2", ping_time_da) def ek_use_platform_angles(platform_ds: xr.Dataset, ping_time_da: xr.DataArray) -> xr.DataArray: @@ -95,7 +71,7 @@ def ek_use_platform_angles(platform_ds: xr.Dataset, ping_time_da: xr.DataArray) echo_range_scaling, dims="time2", coords={"time2": platform_ds["time2"]} ) - return _var_time2_to_ping_time(echo_range_scaling, ping_time_da) + return align_to_ping_time(echo_range_scaling, "time2", ping_time_da) def ek_use_beam_angles( diff --git a/echopype/consolidate/loc_utils.py b/echopype/consolidate/loc_utils.py index 20b77ec41..daeab0f09 100644 --- a/echopype/consolidate/loc_utils.py +++ b/echopype/consolidate/loc_utils.py @@ -4,6 +4,7 @@ import xarray as xr from ..echodata import EchoData +from ..utils.align import align_to_ping_time from ..utils.log import _init_logger logger = _init_logger(__name__) @@ -162,4 +163,4 @@ def sel_interp( ) else: # Values may be nan if there are ping_time values outside the time_dim_name range - return position_var.interp(**{time_dim_name: ds["ping_time"]}) + return align_to_ping_time(position_var, time_dim_name, ds["ping_time"], "linear") diff --git a/echopype/tests/calibrate/test_env_params.py b/echopype/tests/calibrate/test_env_params.py index 98e31d74b..a359a68f7 100644 --- a/echopype/tests/calibrate/test_env_params.py +++ b/echopype/tests/calibrate/test_env_params.py @@ -30,6 +30,7 @@ def ek80_cal_path(test_path): return test_path["EK80_CAL"] +@pytest.mark.unit def test_harmonize_env_param_time(): # Scalar p = 10.05 @@ -67,16 +68,6 @@ def test_harmonize_env_param_time(): assert (p_new.data == p.data).all() # ping_time target requires actual interpolation - ping_time_target = xr.DataArray( - data=[1, 2], - coords={ - "ping_time": np.array( - ["2017-06-20T01:00:15", "2017-06-20T01:00:30"], dtype="datetime64[ns]" - ) - }, - dims=["ping_time"], - ) - ping_time_target = xr.DataArray( data=[1], coords={ diff --git a/echopype/tests/consolidate/test_add_location.py b/echopype/tests/consolidate/test_add_location.py index edc9a418a..81495fb7a 100644 --- a/echopype/tests/consolidate/test_add_location.py +++ b/echopype/tests/consolidate/test_add_location.py @@ -162,7 +162,11 @@ def _tests(ds_test, location_type, nmea_sentence=None): position_var = ed["Platform"][ed_position] if nmea_sentence: position_var = position_var[ed["Platform"]["sentence_type"] == nmea_sentence] - position_interp = position_var.interp(time1=ds_test["ping_time"]) + position_interp = position_var.interp( + {"time1": ds_test["ping_time"]}, + method="nearest", + kwargs={"fill_value": "extrapolate"}, + ) # interpolated values are identical assert np.allclose(ds_test[ds_position].values, position_interp.values, equal_nan=True) # noqa elif location_type == "fixed-location": diff --git a/echopype/tests/utils/test_align.py b/echopype/tests/utils/test_align.py new file mode 100644 index 000000000..ef625fd26 --- /dev/null +++ b/echopype/tests/utils/test_align.py @@ -0,0 +1,48 @@ +import xarray as xr +import numpy as np +import pytest + +import echopype as ep +from echopype.utils.align import align_to_ping_time + + +@pytest.mark.unit +def test_align_to_ping_time(): + """ + Test aligning external pitch data to Echosounder ping time + """ + # Open RAW and extract ping time + ping_time_da = ep.open_raw( + raw_file="echopype/test_data/azfp/rutgers_glider_external_nc/18011107.01A", + xml_path="echopype/test_data/azfp/rutgers_glider_external_nc/18011107.XML", + sonar_model="azfp" + )["Sonar/Beam_group1"]["ping_time"] + + # Open external glider dataset + glider_ds = xr.open_dataset( + "echopype/test_data/azfp/rutgers_glider_external_nc/ru32-20180109T0531-profile-sci-delayed-subset.nc", + engine="netcdf4" + ) + + # Drop NaNs from pitch and align data array + aligned_da = align_to_ping_time(glider_ds["m_pitch"].dropna("time"), "time", ping_time_da) + + # Extract earliest non-NaN time of the aligned dataset + earliest_non_nan_time_in_aligned_da = aligned_da.dropna("ping_time")["ping_time"].min() + + # Grab all values past the earliest non-NaN time in the aligned DataArray + subset_aligned_da = aligned_da.where(aligned_da["ping_time"] >= earliest_non_nan_time_in_aligned_da, drop=True) + + # Check that all values past the earliest non-NaN time in the aligned DataArray are non-NaN + assert np.all(~np.isnan(subset_aligned_da)) + + # Test if aligned_da matches interpolation + assert np.allclose( + aligned_da, + glider_ds["m_pitch"].dropna("time").interp( + {"time": ping_time_da}, + method="nearest", + kwargs={"fill_value": "extrapolate"}, + ), + equal_nan=True, + ) diff --git a/echopype/utils/align.py b/echopype/utils/align.py new file mode 100644 index 000000000..9d0c8c4aa --- /dev/null +++ b/echopype/utils/align.py @@ -0,0 +1,42 @@ +import xarray as xr + + +def align_to_ping_time( + external_da: xr.DataArray, + external_time_name: str, + ping_time_da: xr.DataArray, + method: str = "nearest", +) -> xr.DataArray: + """ + Aligns an external DataArray to align time-wise with the Echosounder ping time DataArray. + + A wrapper function for https://docs.xarray.dev/en/stable/generated/xarray.DataArray.interp.html. + + Parameters + ---------- + external_da : xr.DataArray + External Non-Echosounder data. + external_time_name : str, default 'nearest' + Time variable name of the External Non-Echosounder data. + ping_time_da : xr.DataArray + Echosounder ping time. + + Returns + ------- + aligned_da : xr.DataArray + External Non-Echosounder data that is now aligned with the Echosounder ping time. + """ + # Interpolate only if ping time and external time are not equal + if not ping_time_da.equals( + external_da[external_time_name].rename({external_time_name: "ping_time"}) + ): + aligned_da = external_da.interp( + {external_time_name: ping_time_da}, + method=method, + # More details for `fill_value` and `extrapolate` can be found here: + # https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html # noqa + kwargs={"fill_value": "extrapolate"}, + ).drop_vars(external_time_name) + else: + aligned_da = external_da.rename({external_time_name: "ping_time"}) + return aligned_da