Repository: OSOceanAcoustics/echopype (language: python)

Issue:
Errors/issues processing sequentially pinged multifrequency broadband EK80 files
Hi all. I didn't see this in either of the current open or closed issues, but I am running into a data processing issue with Simrad *.raw files containing broadband multiple frequency data. I followed the threads and chain of issues/pulls in both [#308](https://github.com/OSOceanAcoustics/echopype/issues/308#issue-870570268) and [#657](https://github.com/OSOceanAcoustics/echopype/pull/657), but those didn't seem to directly help address the issues I was encountering. I apologize if I missed a thread that does, in fact, address this exact issue. Namely, some of the issues addressed between 0.5.6 and 0.6.0 didn't appear to resolve errors involving changing data values for `transmit_power`, `slope`, `transmit_duration_nominal`, `frequency_start`, and `frequency_end`.

**Background**

**OS:** Windows 10.0.19043
**Python version:** 3.7.9 _and_ 3.8.10
**Echopype version:** 0.5.6 _and_ 0.6.0
**Simrad EK80 software version:** 2.0.0

```
ed["Provenance"]

<xarray.Dataset>
Dimensions:           (filenames: 1)
Dimensions without coordinates: filenames
Data variables:
    source_filenames  (filenames) <U57 'C:/Users/Brandyn/Downloads/NYOS2105-D...
Attributes:
    conversion_software_name:     echopype
    conversion_software_version:  0.6.0
    conversion_time:              2022-06-28T18:22:15Z
```

An aside: I ran `Echopype 0.5.6` in `Python 3.7.9` and `Echopype 0.6.0` in `Python 3.8.10`. For `Python 3.7.x`, there is an error in trying to install `Echopype 0.6.0` due to certain dependencies requiring `Python >= 3.8` (namely `xarray`). So the recommended Python version in the [Installation Guide](https://echopype.readthedocs.io/en/stable/installation.html) may need to be amended for Windows. This was detected via: 

`pip install git+https://github.com/OSOceanAcoustics/echopype.git#egg=echopype`

**Issue**
I am trying to call in EK80 raw files that includes four frequencies (38, 70, 120, 200 kHz) where the 38 and 200 kHz come from a dual frequency transducer (38/200 Combi D). All four comprise broadband/complex measurements. The fileset I am currently processing is from May 2021. I haven't run into any significant issues when trying to open the raw files. However, I do run into an issue when trying to compute Sv or attempt really any other processing function. This also includes trying to visualize the echogram. I have been able to successfully run the examples from the provided [notebook](https://osoceanacoustics.github.io/echopype-examples/index.html) using EK60 data.

The issue appears to be related to how some of the arrays for certain variables are mapped to the `frequency` and `ping_time` coordinates. So the initial Error message I receive is as follows (**Note**: ultimately, the flagged error is identical in both `Echopype 0.5.6` and `Echopype 0.6.0`): 

```
import echopype as ep
from echopype import open_raw
import numpy as np

file = "C:/Users/Brandyn/Downloads/NYOS2105-D20210525-T213648.raw"
ed = open_raw(file, sonar_model='EK80')
ed_sv = ep.calibrate.compute_Sv(ed, waveform_mode="BB", encode_mode="complex")
```
_Resulting error_

```
  File "<stdin>", line 1, in <module>
  File "C:\Users\Brandyn\AppData\Roaming\Python\Python38\site-packages\echopype\calibrate\api.py", line 193, in compute_Sv
    return _compute_cal(cal_type="Sv", echodata=echodata, **kwargs)
  File "C:\Users\Brandyn\AppData\Roaming\Python\Python38\site-packages\echopype\calibrate\api.py", line 81, in _compute_cal
    cal_ds = cal_obj.compute_Sv(waveform_mode=waveform_mode, encode_mode=encode_mode)
  File "C:\Users\Brandyn\AppData\Roaming\Python\Python38\site-packages\echopype\calibrate\calibrate_ek.py", line 964, in compute_Sv
    return self._compute_cal(
  File "C:\Users\Brandyn\AppData\Roaming\Python\Python38\site-packages\echopype\calibrate\calibrate_ek.py", line 930, in _compute_cal
    ds_cal = self._cal_complex(cal_type=cal_type, waveform_mode=waveform_mode)
  File "C:\Users\Brandyn\AppData\Roaming\Python\Python38\site-packages\echopype\calibrate\calibrate_ek.py", line 684, in _cal_complex
    chirp, _, tau_effective = self.get_transmit_chirp(waveform_mode=waveform_mode)
  File "C:\Users\Brandyn\AppData\Roaming\Python\Python38\site-packages\echopype\calibrate\calibrate_ek.py", line 549, in get_transmit_chirp
    raise TypeError("File contains changing %s!" % p)
TypeError: File contains changing transmit_duration_nominal!
```

Looking under the hood in `beam` from my `EchoData` object (I used the equivalent indexing notation, e.g. `ed.beam`, for the `Echopype 0.5.6` release): 

```
print(ed['Sonar/Beam_group1'])

<xarray.Dataset>
Dimensions:                        (channel: 4, ping_time: 361, beam: 4,
                                    range_sample: 16729)
Coordinates:
  * channel                        (channel) <U36 'WBT 400039-15 ES120-7C_ES'...
  * ping_time                      (ping_time) datetime64[ns] 2021-05-25T21:3...
  * range_sample                   (range_sample) int64 0 1 2 ... 16727 16728
  * beam                           (beam) <U1 '1' '2' '3' '4'
Data variables: (12/21)
    frequency_nominal              (channel) float64 1.2e+05 7e+04 3.8e+04 2e+05
    beam_type                      (channel, ping_time) int32 1 1 1 1 ... 0 0 0
    beamwidth_twoway_alongship     (channel, ping_time, beam) float64 7.0 ......
    beamwidth_twoway_athwartship   (channel, ping_time, beam) float64 7.0 ......
    beam_direction_x               (channel, ping_time, beam) float64 0.0 ......
    beam_direction_y               (channel, ping_time, beam) float64 0.0 ......
    ...                             ...
    frequency_start                (channel, ping_time, beam) float64 nan ......
    frequency_end                  (channel, ping_time, beam) float64 nan ......
    sample_interval                (channel, ping_time) float64 nan ... 1.6e-05
    transmit_power                 (channel, ping_time) float64 nan ... 200.0
    transmit_duration_nominal      (channel, ping_time) float32 nan ... 0.000512
    slope                          (channel, ping_time) float64 nan ... 0.0514
Attributes:
    beam_mode:              vertical
    conversion_equation_t:  type_3
```

```
print(ed['Sonar/Beam_group1']['transmit_duration_nominal'])

<xarray.DataArray 'transmit_duration_nominal' (channel: 4, ping_time: 361)>
array([[     nan,      nan, 0.000512, ...,      nan, 0.000512,      nan],
       [     nan, 0.000512,      nan, ..., 0.000512,      nan,      nan],
       [0.000512,      nan,      nan, ...,      nan,      nan, 0.000512],
       [0.000512,      nan,      nan, ...,      nan,      nan, 0.000512]],
      dtype=float32)
Coordinates:
  * channel    (channel) <U36 'WBT 400039-15 ES120-7C_ES' ... 'WBT Mini 26118...
  * ping_time  (ping_time) datetime64[ns] 2021-05-25T21:36:48.192999936 ... 2...
Attributes:
    long_name:  Nominal bandwidth of transmitted pulse
    units:      s
    valid_min:  0.0
```
I wasn't sure if the padded `nan` values here were related to [#681](https://github.com/OSOceanAcoustics/echopype/issues/681#issuecomment-1133088359), but when checking whether there is a different transmit duration included in this array: 

```
ed['Sonar/Beam_group1']['transmit_duration_nominal'].shape
(4, 361)
np.unique(ed['Sonar/Beam_group1']['transmit_duration_nominal'].values)
array([0.000512,      nan], dtype=float32)
```

It makes sense since pings were sequential so the `ping_time` coordinate won't include actual data for each frequency. Looking at the first 6 pings: 

```
ed['Sonar/Beam_group1']["transmit_duration_nominal"].values[:, 0:6]

array([[0.000512,      nan,      nan, 0.000512,      nan,      nan],
       [     nan, 0.000512,      nan,      nan, 0.000512,      nan],
       [     nan,      nan, 0.000512,      nan,      nan, 0.000512],
       [0.000512,      nan,      nan, 0.000512,      nan,      nan]],
      dtype=float32)
```

So that matches what would be expected since the 38/200 pinged simultaneously followed by the 70 and 120 kHz transducers, respectively. So, presumably, this should be the same for other variables mapped with the same coordinates (i.e. `frequency` and `ping_time`), which it is (e.g. `slope`): 

```
ed['Sonar/Beam_group1']["slope"].values[:, 0:6]
array([[0.223214 ,       nan,       nan, 0.223214 ,       nan,       nan],
       [      nan, 0.0710227,       nan,       nan, 0.0710227,       nan],
       [      nan,       nan, 0.0434028,       nan,       nan, 0.0434028],
       [0.051398 ,       nan,       nan, 0.051398 ,       nan,       nan]])
```

So when I just replace the `nan` placeholders with the same `float32` value, the issue gets passed onto the next variable (so `slope` in this case): 

```
beam_tdn_copy = ed['Sonar/Beam_group1']['transmit_duration_nominal'].copy()
tdn_vals = np.unique(beam_tdn_copy.values.flatten(), return_index=True)
tdn_vals = tdn_vals[0][tdn_vals[1].argsort()]
tdn_vals = tdn_vals[~np.isnan(tdn_vals)]
ed['Sonar/Beam_group1']['transmit_duration_nominal'].values = np.full_like(beam_tdn_copy, tdn_vals[:, np.newaxis])
print(ed['Sonar/Beam_group1']["transmit_duration_nominal"])

<xarray.DataArray 'transmit_duration_nominal' (channel: 4, ping_time: 361)>
array([[0.000512, 0.000512, 0.000512, ..., 0.000512, 0.000512, 0.000512],
       [0.000512, 0.000512, 0.000512, ..., 0.000512, 0.000512, 0.000512],
       [0.000512, 0.000512, 0.000512, ..., 0.000512, 0.000512, 0.000512],
       [0.000512, 0.000512, 0.000512, ..., 0.000512, 0.000512, 0.000512]],
      dtype=float32)
Coordinates:
  * channel    (channel) <U36 'WBT 400039-15 ES120-7C_ES' ... 'WBT Mini 26118...
  * ping_time  (ping_time) datetime64[ns] 2021-05-25T21:36:48.192999936 ... 2...
Attributes:
    long_name:  Nominal bandwidth of transmitted pulse
    units:      s
    valid_min:  0.0
```
```
ed_sv = ep.calibrate.compute_Sv(ed, waveform_mode='BB', encode_mode='complex')
TypeError: File contains changing slope!
```
Changing slope: 

```
print(ed['Sonar/Beam_group1']["slope"])

<xarray.DataArray 'slope' (channel: 4, ping_time: 361)>
array([[      nan,       nan, 0.0434028, ...,       nan, 0.0434028,
              nan],
       [      nan, 0.0710227,       nan, ..., 0.0710227,       nan,
              nan],
       [0.223214 ,       nan,       nan, ...,       nan,       nan,
        0.223214 ],
       [0.051398 ,       nan,       nan, ...,       nan,       nan,
        0.051398 ]])
Coordinates:
  * channel    (channel) <U36 'WBT 400039-15 ES120-7C_ES' ... 'WBT Mini 26118...
  * ping_time  (ping_time) datetime64[ns] 2021-05-25T21:36:48.192999936 ... 2...
```

```
slope_vals = np.unique(ed['Sonar/Beam_group1']["slope"].values.flatten(), return_index=True)
slope_vals = slope_vals[0][slope_vals[1].argsort()]
slope_vals = slope_vals[~np.isnan(slope_vals)]
ed['Sonar/Beam_group1']["slope"].values = np.full_like(ed['Sonar/Beam_group1']["slope"], slope_vals[:, np.newaxis])
print(ed['Sonar/Beam_group1']["slope"])

<xarray.DataArray 'slope' (channel: 4, ping_time: 361)>
array([[0.0434028, 0.0434028, 0.0434028, ..., 0.0434028, 0.0434028,
        0.0434028],
       [0.0710227, 0.0710227, 0.0710227, ..., 0.0710227, 0.0710227,
        0.0710227],
       [0.223214 , 0.223214 , 0.223214 , ..., 0.223214 , 0.223214 ,
        0.223214 ],
       [0.051398 , 0.051398 , 0.051398 , ..., 0.051398 , 0.051398 ,
        0.051398 ]])
Coordinates:
  * channel    (channel) <U36 'WBT 400039-15 ES120-7C_ES' ... 'WBT Mini 26118...
  * ping_time  (ping_time) datetime64[ns] 2021-05-25T21:36:48.192999936 ... 2...
```

Same error, just for `slope` now:

```
ed_sv = ep.calibrate.compute_Sv(ed, waveform_mode='BB', encode_mode='complex')
TypeError: File contains changing transmit_power!
```
Repeating this by replacing the `nan` placeholders in `transmit_power`, `frequency_start`, and `frequency_end` (the latter two required a slightly different brute force approach since they apparently have a different dimension than the other data variables listed above) appears to resolve the issue since functions like `compute_Sv(...)` will now operate as expected:

```
ed_sv = ep.calibrate.compute_Sv(ed, waveform_mode='BB', encode_mode='complex')
print(ed_sv)
<xarray.Dataset>
Dimensions:                (ping_time: 361, channel: 4, range_sample: 16729,
                            filenames: 1, time3: 1)
Coordinates:
  * ping_time              (ping_time) datetime64[ns] 2021-05-25T21:36:48.192...
  * channel                (channel) <U36 'WBT 400039-15 ES120-7C_ES' ... 'WB...
  * range_sample           (range_sample) int64 0 1 2 3 ... 16726 16727 16728
  * time3                  (time3) datetime64[ns] 2021-05-25T21:36:48.192999936
Dimensions without coordinates: filenames
Data variables: (12/13)
    Sv                     (channel, ping_time, range_sample) float64 nan ......
    echo_range             (channel, ping_time, range_sample) float64 nan ......
    frequency_nominal      (channel) float64 1.2e+05 7e+04 3.8e+04 2e+05
    temperature            (ping_time) float64 10.0 nan nan nan ... nan nan nan
    salinity               (ping_time) float64 35.0 nan nan nan ... nan nan nan
    pressure               (ping_time) float64 100.0 nan nan nan ... nan nan nan
    ...                     ...
    sound_absorption       (ping_time, channel) float64 0.04094 0.02458 ... nan
    sa_correction          object None
    gain_correction        object None
    equivalent_beam_angle  (channel, ping_time) float64 -20.7 -20.7 ... -12.6
    source_filenames       (filenames) <U57 'C:/Users/Brandyn/Downloads/NYOS2...
    water_level            (time3) float64 0.0
Attributes:
    processing_software_name:     echopype
    processing_software_version:  0.6.0
    processing_time:              2022-06-28T18:33:22Z
    processing_function:          calibrate.compute_Sv
```

But this hack to brute force `compute_Sv(...)` unfortunately appears to potentially break other features (or this may be just a general 'feature' on how these *.raw files are formatted regardless of the above hack). For instance, I am able to generate the unbinned Sv echograms without much trouble (attached is the related figure for the 70 kHz): 

```
import matplotlib.pyplot as plt 
import echopype.visualize as epviz

gg = epviz.create_echogram(ed_sv, 
                      robust=True, 
                      get_range=True, 
                      range_kwargs=dict(waveform_mode='BB', encode_mode='complex'),
                      vmin=-70, 
                      vmax=-36)
plt.show(gg[1])
plt.close()
```
![Lucca_echogram_test_plot_unbinned_70kHz](https://user-images.githubusercontent.com/45861253/176272191-83fc0b62-1ed7-43c7-b814-36d1a4e7e7ad.jpeg)

I can manually toggle through each of them to pull the channel/frequency I would like. But as soon as I average by bin via `compute_MVBS(...)`, things seem to break for all frequencies except 38 kHz where the echogram is able to be generated. 

```
ed_mvbs = ep.preprocess.compute_MVBS(ed_sv, range_meter_bin=1, ping_time_bin='10S')
gg = epviz.create_echogram(ed_mvbs, 
                      robust=True, 
                      get_range=True, 
                      range_kwargs=dict(waveform_mode='BB', encode_mode='complex'),
                      vmin=-70, 
                      vmax=-36)
plt.show(gg[2])
plt.close()
```
![Lucca_echogram_test_plot_binnedmvbs_70kHz](https://user-images.githubusercontent.com/45861253/176273551-826c41ef-4f3c-46a0-9c7d-1acfc0f8393c.jpeg)

It's certainly quite possible that I broke something in the process and that is what is triggering the above issue, but I have replicated the exact issue and 'hack' in 10 unique EK *.raw files from other cruises where we also collected broadband measurements with sequential pings. Any guidance would be greatly appreciated! 
harmonize_env_param_time fails with index 0 error
Hi, 

When I attempt to `compute_Sv`:
'''
ep.calibrate.compute_Sv(echodata[0],waveform_mode='BB',encode_mode='complex')
'''

I am receiving an error: `IndexError: index 0 is out of bounds for axis 0 with size 0`, which is coming from the attempt to squeeze the environmental parameters object after dropping the `time1` column in `harmonize_env_param_time()` (line 56). 

I am new to `echopype` so I may be missing something obvious, or perhaps something is odd about our data files. Perhaps this view of the `environment` object will be instructive for someone with more knowledge:

![Notification_Center](https://github.com/OSOceanAcoustics/echopype/assets/2421975/d2d63477-840f-45bf-9697-4ab49740870e)



Discussion / hints:
@brandynlucca : Thanks for reporting the bugs and the very thorough investigation!! I believe you didn't break anything and it is the code that needs improvements. 

The NaN padding is a strategy we use to get data from all channels into a gridded form for slicing, and it is great (!) to know that the current calibration routine breaks with the interleaving NaNs for sequential pings. I quickly looked into the section where things broke and saw a TODO for exactly the problem you saw:
https://github.com/OSOceanAcoustics/echopype/blob/12693a72badcf54a18408e694b72e7b0e6d48023/echopype/calibrate/calibrate_ek.py#L528-L529

I will need probably until the weekend to get to this. If you wouldn't mind sharing a small example file that will really help with changing/debugging the implementation. This gives me a push to look into the EK80 cal code too, so thanks! (I think in addition to this issue the current implementation is also pretty slow due to the mechanism used to do pulse compression)

For the MVBS problem, it is a regression bug at v0.6.0 that should be fixed in #736 (but not merged yet). Could you give it a try to see if that resolved the problem?

Also, thanks for pointing out the >=3.8 problem. We'll get that fixed!
> @brandynlucca : Thanks for reporting the bugs and the very thorough investigation!! I believe you didn't break anything and it is the code that needs improvements.
> 
> The NaN padding is a strategy we use to get data from all channels into a gridded form for slicing, and it is great (!) to know that the current calibration routine breaks with the interleaving NaNs for sequential pings. I quickly looked into the section where things broke and saw a TODO for exactly the problem you saw:
> 
> https://github.com/OSOceanAcoustics/echopype/blob/12693a72badcf54a18408e694b72e7b0e6d48023/echopype/calibrate/calibrate_ek.py#L528-L529
> 
> I will need probably until the weekend to get to this. If you wouldn't mind sharing a small example file that will really help with changing/debugging the implementation. This gives me a push to look into the EK80 cal code too, so thanks! (I think in addition to this issue the current implementation is also pretty slow due to the mechanism used to do pulse compression)

No worries! I wonder if one plausible alternative would be to generate multidimensional coordinates that maps `ping_time` to a specific `channel`? Perhaps that would also allow frequency-specific data to be severable from the overall dataset (e.g. processing just 38/200 kHz backscatter via `Echopype` for something like dB differencing without requiring external manipulation or needing to process other unused frequency data that will just chew up memory and resources)? 

In terms of files, I unfortunately only have 100 MB RAW files, but I can upload on to my Google drive and share it if you'd like! 

> For the MVBS problem, it is a regression bug at v0.6.0 that should be fixed in #736 (but not merged yet). Could you give it a try to see if that resolved the problem?

I tested the dataflow by installing the updated version included in your repo but that didn't seem to fix my particular issue it seems. So I am able to generate the 38 kHz echograms without issue, but the other three frequencies still aren't bin-averaging correctly. I wonder if this is related to the 38 kHz data occupying both the first non-`nan`-padded column (i.e. ping 0) _and_ row? I know `xarray` can be particularly persnickety when it comes to reindexing things, so maybe that could explain why `compute_MVBS(...)` works for the first frequency channel and just happens to produce _something_ for the other three since they share the same geometries by coincidence? Running down the rabbit hole a bit seems to indicate that this may still be related to how `echo_range` is being calculated since the issue clears up when brute-forcing the code to resample vertically by `range_sample_num` and horizontally by `ping_time`. 

All that said, at least from the linked repo with the updated `compute_MVBS(...)` function I am actually able to manually construct echograms for the other frequencies by extracting their frequency-specific `ping_time` measurements with most of the same original code. The very unoptimized code below: 

_Original dataset_
```
sv = 10 ** (ed_sv["Sv"]/10)
print(sv)
```

```
<xarray.DataArray 'Sv' (channel: 4, ping_time: 361, range_sample: 16729)>
array([[[           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
         3.91095702e-07, 9.14494970e-07, 1.85589463e-07],
        ...,
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
         2.49507722e-06, 1.44740865e-06, 2.45585001e-06],
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan]],

       [[           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
...
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan]],

       [[           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        ...,
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan],
        [           nan,            nan,            nan, ...,
                    nan,            nan,            nan]]])
Coordinates:
  * ping_time     (ping_time) datetime64[ns] 2021-05-25T21:36:48.192999936 .....
  * channel       (channel) <U36 'WBT 400039-15 ES120-7C_ES' ... 'WBT Mini 26...
  * range_sample  (range_sample) int64 0 1 2 3 4 ... 16725 16726 16727 16728
```

I then glued together a bitmap function that subsets out the frequency-specific data's coordinates:

```
def _frequency_bitmap(sv_data, frequency): 
  #Add frequency as coordinate
  sv_data.coords["frequency"] = (
    ["channel"],
    ed_sv["frequency_nominal"].data
  )
  sv_data = sv_data.swap_dims({"channel": "frequency"})
  #Create bitmap of presence/absence of non-padded data
  bool_array = (~np.isnan(sv_data.sel(frequency=frequency).data)).sum(1)
  #Parse out ping times specific to frequency argument input
  ping_times_f = sv_data["ping_time"][bool_array > 0]
  #Select subset of original dataset
  return sv_data.loc[dict(frequency=frequency, ping_time=ping_times_f)]
```

So when applied to data to pull a specific frequency, the object ends up looking like: 

```
#Set frequency choice 
frequency_choice = 7.0e+4
test_obj = _frequency_bitmap(sv, frequency_choice)
print(test_obj)
```
```
<xarray.DataArray 'Sv' (ping_time: 120, range_meter: 16729)>
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]])
Coordinates:
  * ping_time     (ping_time) datetime64[ns] 2021-05-25T21:36:50.172000256 .....
    channel       <U36 'WBT 738184-15 ES70-7CD_ES'
    range_sample  (range_meter) int64 0 1 2 3 4 ... 16725 16726 16727 16728
  * range_meter   (range_meter) float64 0.0 0.0 0.0 0.0 0.0 ... nan nan nan nan
    frequency     float64 7e+04
```

So with that operationalized, I was able to generate the corresponding MVBS echogram for the other frequencies via:

```
#####
#Set frequency choice 
frequency_choice = 7.0e+4
#Set range bin (m)
range_meter_bin = 1.0 
#Set ping time bin (S)
ping_bin = "10S"
#####
sv = 10 ** (ed_sv["Sv"] / 10)
freq_ind_mask = np.unique(ds_Sv.frequency_nominal.values.flatten(), return_index=True)[0] == frequency_choice
freq_ind = np.unique(ds_Sv.frequency_nominal.values.flatten(), return_index=True)[1][freq_ind_mask]
sorted_freq_ind = np.argsort(ds_Sv.frequency_nominal)
sv.coords["range_meter"] = (
    ["range_sample"],
    ed_sv["echo_range"].isel(channel=sorted_freq_ind[freq_ind[0]], ping_time=freq_ind[0]).data,
)
sv = sv.swap_dims({"range_sample": "range_meter"})
range_interval = np.arange(0, ed_sv["echo_range"].max() + range_meter_bin, range_meter_bin)
rint = range_interval
pbin = ping_bin
sv_groupby_bins = (
  _frequency_bitmap(sv, frequency_choice)
  .groupby_bins("range_meter", bins=rint, right=False, include_lowest=True)
  .mean()
  .resample(ping_time=pbin, skipna=True)
  .mean()
)
sv_groupby_bins.coords["echo_range"] = (["range_meter_bins"], rint[:-1])
sv_groupby_bins = sv_groupby_bins.swap_dims({"range_meter_bins": "echo_range"})
sv_groupby_bins = sv_groupby_bins.drop_vars("range_meter_bins")
sv_groupby_bins = 10 * np.log10(sv_groupby_bins)
gg = sv_groupby_bins.plot.pcolormesh(x='ping_time', y='echo_range', cmap='jet', vmin=-80, vmax=-44)
plt.ylim(45, 2)
plt.show(gg)
plt.close()
```
![Lucca Echopype_test_echogram_MVBS_single_frequency_70kHz](https://user-images.githubusercontent.com/45861253/176537408-abf9cb82-82cc-40eb-9620-0e5447dc807b.jpeg)

Hopefully there is a tidbit in there that may at least provide some extra inkling into the source of the issue with this dataset? Again, I can provide a link to the example file if needed.
> No worries! I wonder if one plausible alternative would be to generate multidimensional coordinates that maps `ping_time` to a specific `channel`? Perhaps that would also allow frequency-specific data to be severable from the overall dataset (e.g. processing just 38/200 kHz backscatter via `Echopype` for something like dB differencing without requiring external manipulation or needing to process other unused frequency data that will just chew up memory and resources)?

Yes! this is perhaps the nice use case for using dask.delayed via accessing already existing zarr store (so just slice out the frequencies needed) for distributing the computation. In a way once we are processing large files an intermediate zarr store may always be needed. In this particular case we could slice out just 38/200 and then drop all NaNs for the calculation. Great! I'll give this a try when revising the calibration code, perhaps it'll be some combination with groupby to make it channel-independent and delay-friendly. 

> 
> In terms of files, I unfortunately only have 100 MB RAW files, but I can upload on to my Google drive and share it if you'd like!

Please do! 100 MB is not a problem for developing. I may be able to reduce the file size via replaying it through EK80 (if you have access to the EK80 software perhaps you could try that too), though I don't know if that may change some detailed info of the file (like software version etc) -- but I guess for using it as a test file in the CI for making sure the functionalities work it would be fine to have those different.


And thanks for providing the code that made things work! I am sure somewhere in there the issue is illuminated. I will look more closely in a couple of days.
> Please do! 100 MB is not a problem for developing. I may be able to reduce the file size via replaying it through EK80 (if you have access to the EK80 software perhaps you could try that too), though I don't know if that may change some detailed info of the file (like software version etc) -- but I guess for using it as a test file in the CI for making sure the functionalities work it would be fine to have those different.

Here is a Google Drive link to the *.raw file I was using:

https://drive.google.com/file/d/1c1d0rUgBa9N-2UfZOM7_rnCSvuhb0g0H/view?usp=sharing

Let me know if there are any issues accessing it! 
This issue might be related to https://github.com/OSOceanAcoustics/echopype/issues/720#issuecomment-1165568926
@brandynlucca : we fixed a few regression bugs from v0.6.0 in the most recent release (v0.6.1 this morning), but I haven't had a chance to fix the part related to the NaN issues during BB calibration. Stay tuned -- I'll ping when there's PR for that.
I have a similar problem with WBAT-mini data (EK80 format), in CW mode,

```
>>> file='data/MultiFrequency-4min-38N200W-continuous38WN-200W-38N-SOFS-8-Phase0-D20190401-T000014-0.raw'
>>> ed = open_raw(file, sonar_model='EK80')
/Users/pete/DWM/echo/venv/lib/python3.8/site-packages/xarray/coding/times.py:228: RuntimeWarning: invalid value encountered in cast
  flat_num_dates_ns_int = (flat_num_dates * _NS_PER_TIME_DELTA[delta]).astype(
/Users/pete/DWM/echo/venv/lib/python3.8/site-packages/xarray/coding/times.py:228: RuntimeWarning: invalid value encountered in cast
  flat_num_dates_ns_int = (flat_num_dates * _NS_PER_TIME_DELTA[delta]).astype(
/Users/pete/DWM/echo/venv/lib/python3.8/site-packages/xarray/coding/times.py:228: RuntimeWarning: invalid value encountered in cast
  flat_num_dates_ns_int = (flat_num_dates * _NS_PER_TIME_DELTA[delta]).astype(
>>> ed_sv = ep.calibrate.compute_Sv(ed, waveform_mode="CW", encode_mode="complex")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/api.py", line 205, in compute_Sv
    return _compute_cal(cal_type="Sv", echodata=echodata, **kwargs)
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/api.py", line 89, in _compute_cal
    cal_ds = cal_obj.compute_Sv(waveform_mode=waveform_mode, encode_mode=encode_mode)
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/calibrate_ek.py", line 997, in compute_Sv
    return self._compute_cal(
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/calibrate_ek.py", line 963, in _compute_cal
    ds_cal = self._cal_complex(cal_type=cal_type, waveform_mode=waveform_mode)
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/calibrate_ek.py", line 693, in _cal_complex
    chirp, _, tau_effective = self.get_transmit_chirp(waveform_mode=waveform_mode)
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/calibrate_ek.py", line 554, in get_transmit_chirp
    raise TypeError("File contains changing %s!" % p)
TypeError: File contains changing transmit_duration_nominal!
>>> print(ed['Sonar/Beam_group1'])
<xarray.Dataset>
Dimensions:                        (channel: 3, ping_time: 4, beam: 4,
                                    range_sample: 5750)
Coordinates:
  * channel                        (channel) <U29 'EKA 266958-0F ES38D' ... '...
  * ping_time                      (ping_time) datetime64[ns] 2019-04-01T00:0...
  * range_sample                   (range_sample) int64 0 1 2 ... 5747 5748 5749
  * beam                           (beam) <U1 '1' '2' '3' '4'
Data variables: (12/19)
    frequency_nominal              (channel) float64 3.8e+04 3.8e+04 2e+05
    beam_type                      (channel, ping_time) int64 1 1 1 1 ... 0 0 0
    beamwidth_twoway_alongship     (channel, ping_time, beam) float64 7.0 ......
    beamwidth_twoway_athwartship   (channel, ping_time, beam) float64 7.0 ......
    beam_direction_x               (channel, ping_time, beam) float64 nan ......
    beam_direction_y               (channel, ping_time, beam) float64 nan ......
    ...                             ...
    backscatter_r                  (channel, ping_time, range_sample, beam) float32 ...
    backscatter_i                  (channel, ping_time, range_sample, beam) float32 ...
    sample_interval                (channel, ping_time) float64 nan ... nan
    transmit_power                 (channel, ping_time) float64 nan ... nan
    transmit_duration_nominal      (channel, ping_time) float32 nan ... nan
    slope                          (channel, ping_time) float64 nan ... nan
Attributes:
    beam_mode:              vertical
    conversion_equation_t:  type_3
>>> 
```

Would an option be to have a channel option to calibrate.compute_Sv to process a single channel to Sv? to Combine them after they are put onto the same grid with compute_MVBS.

Attached 2 example data files 
[MultiFrequency-WBAT-mini.zip](https://github.com/OSOceanAcoustics/echopype/files/10952259/MultiFrequency-WBAT-mini.zip)

Get a different error with the current dev branch,

```
ds_Sv = ep.calibrate.compute_Sv(ed, waveform_mode='CW', encode_mode='complex')
Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/api.py", line 204, in compute_Sv
    return _compute_cal(cal_type="Sv", echodata=echodata, **kwargs)
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/api.py", line 58, in _compute_cal
    cal_ds = cal_obj.compute_Sv()
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/calibrate_ek.py", line 490, in compute_Sv
    return self._compute_cal(cal_type="Sv")
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/calibrate_ek.py", line 474, in _compute_cal
    ds_cal = self._cal_complex_samples(cal_type=cal_type, complex_ed_group=self.ed_group)
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/calibrate_ek.py", line 355, in _cal_complex_samples
    tx_coeff = get_filter_coeff(vend)
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/ek80_complex.py", line 158, in get_filter_coeff
    coeff[ch_id]["wbt_fil"] = get_vend_filter_EK80(vend, ch_id, "WBT", "coeff")
  File "/Users/pete/DWM/echo/echopype/echopype/calibrate/ek80_complex.py", line 130, in get_vend_filter_EK80
    v = vend.attrs[f"{channel_id} {filter_name} filter_r"] + 1j * np.array(
KeyError: 'EKA 266958-0F ES38D WBT filter_r'
print(ep.__version__)
0.6.4.dev86+g4ee4844
```
@valschmidt : I wonder if a bug got accidentally introduced in a recent PR attempting to address some performance issue with `harmonize_env_param_time`. Could you send us a .raw file or a saved EchoData zarr/netcdf file?
Hey @valschmidt : Pinging again to see if you can send us an example file. We cannot reproduce this error using our own datasets.
Sure! 

Here’s a link to on:
https://universitysystemnh-my.sharepoint.com/:i:/g/personal/vyk2_usnh_edu/EZZHUF_AKZdJsOvwtvKbi-gBrBx_wv1dZzUm5hEmru_tSw?e=tsOfls￼

> On Apr 1, 2024, at 13:33, Wu-Jung Lee ***@***.***> wrote:
> 
> 
> Hey @valschmidt <https://github.com/valschmidt> : Pinging again to see if you can send us an example file. We cannot reproduce this error using our own datasets.
> 
> —
> Reply to this email directly, view it on GitHub <https://github.com/OSOceanAcoustics/echopype/issues/1287#issuecomment-2030215509>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AASPJV7HMTM4NLPJV4W2653Y3GK6VAVCNFSM6AAAAABFBOGDOWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAMZQGIYTKNJQHE>.
> You are receiving this because you were mentioned.
> 

------------------------------------------------------
Val Schmidt (he/him)
Center for Coastal and Ocean Mapping / Joint Hydrographic Center
University of New Hampshire
Chase Ocean Engineering Lab
24 Colovos Road
Durham, NH 03824
e: vschmidt [AT] ccom.unh.edu
m: 614.286.3726





@valschmidt Thanks for the example file! This was pretty fun to dig into.

The error stems from the `transmit_frequency_start` and `transmit_frequency_end` in `ed["Sonar/Beam_group1"]` having NaNs, which in turn messes up the calculation of the `sound_absorption` data variable, and that causes the error in `harmonize_env_param_time`. These NaNs originate from the 2 channels that make up `transmit_frequency_start` and `transmit_frequency_end` having inconsistent `ping_time` values, and this inconsistency originates from the EK80 parser/datagram itself. I created a Gist that expands on this a bit more with the file you sent us: https://gist.github.com/ctuguinay/202141a9c263d6e3d2e7a0e8ff748508

@leewujung Let me know what you both think of this. I think one way to fix this in order to run `compute_Sv` is by manually (post-`.open_raw`) filling in the NaN variables in `transmit_frequency_start` and `transmit_frequency_end` with the appropriate frequency values. As for why the EK80 parser/datagram has this `ping_time` inconsistency, I am not sure.
Oh! This is because we are using the multiplexed ping option of the wbt mini. That’s why ping times for the two channels are different. Before I forget, thank you so much for looking out to this. I really appreciate it.
Ah I see, I think this is the same as or similar to #743. It's so long ago but I haven't been able to carve out time to do. :/

@ctuguinay : could you take a look at that issue as well? This is time for us to address this "bug" since sequential pinging sequences are more common now. I think we probably should to take a 2-stage approach to resolve this:
- first resolve the bug so that at least the Sv values can be computed correctly for this type of data
- at a later time we need to address the performance side of things -- the padded NaNs in the raw-converted data may create a lot of potentially wasteful computations (we can see if that is the case and how bad that is), and I wonder if dropna is properly delayed if the computations would only be on the entries that have actual values. This is likely more involved though, as we know from trying to scale other functions, hence my recommendation to separate the scaling into a separate issue.
@leewujung I agree with this approach and I'll look at Brandyn's issue.

Fix patch:
diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py
index c162e5752..2991ae474 100644
--- a/echopype/calibrate/calibrate_ek.py
+++ b/echopype/calibrate/calibrate_ek.py
@@ -473,6 +473,10 @@ def _get_B_theta_phi_m(self):
         ) ** 2
         B_theta_phi_m = 0.5 * 6.0206 * (fac_along + fac_athwart - 0.18 * fac_along * fac_athwart)
 
+        # Zero out NaNs that appear due to 1) multiplex ping patterns and 2) single beam transducer
+        # systems that don't contain angle offset information.
+        B_theta_phi_m = B_theta_phi_m.fillna(0)
+
         return B_theta_phi_m
 
     def _cal_complex_samples(self, cal_type: str) -> xr.Dataset:
diff --git a/echopype/calibrate/ek80_complex.py b/echopype/calibrate/ek80_complex.py
index 5c93118b1..82bf99d89 100644
--- a/echopype/calibrate/ek80_complex.py
+++ b/echopype/calibrate/ek80_complex.py
@@ -226,11 +226,10 @@ def get_transmit_signal(
     # but keeping this here for use as standalone function
     if waveform_mode == "BB" and np.all(beam["transmit_type"] == "CW"):
         raise TypeError("File does not contain BB mode complex samples!")
-
     # Generate all transmit replica
     y_all = {}
     y_time_all = {}
-    # TODO: expand to deal with the case with varying tx param across ping_time
+    # TODO: expand to deal with the case with varying non-NaN tx param across ping_time
     tx_param_names = [
         "transmit_duration_nominal",
         "slope",
@@ -240,16 +239,18 @@ def get_transmit_signal(
     for ch in beam["channel"].values:
         tx_params = {}
         for p in tx_param_names:
-            tx_params[p] = np.unique(beam[p].sel(channel=ch))
+            # Extract beam values and filter out NaNs
+            beam_values = np.unique(beam[p].sel(channel=ch))
+            # Filter out NaN values
+            beam_values_without_nan = beam_values[~np.isnan(beam_values)]
+            tx_params[p] = beam_values_without_nan
             if tx_params[p].size != 1:
                 raise TypeError("File contains changing %s!" % p)
         fs_chan = fs.sel(channel=ch).data if isinstance(fs, xr.DataArray) else fs
         tx_params["fs"] = fs_chan
         y_ch, _ = tapered_chirp(**tx_params)
-
         # Filter and decimate chirp template
         y_ch, y_tmp_time = filter_decimate_chirp(coeff_ch=coeff[ch], y_ch=y_ch, fs=fs_chan)
-
         # Fill into output dict
         y_all[ch] = y_ch
         y_time_all[ch] = y_tmp_time
@@ -278,9 +279,21 @@ def compress_pulse(backscatter: xr.DataArray, chirp: Dict) -> xr.DataArray:
         # Select channel `chan` and drop the specific beam dimension if all of the values are nan.
         backscatter_chan = backscatter.sel(channel=chan).dropna(dim="beam", how="all")
 
+        # Create NaN mask
+        # If `backscatter_chan` is lazy loaded, then `nan_mask` too will be lazy loaded.
+        nan_mask = np.isnan(backscatter_chan)
+
+        # Zero out backscatter NaN values
+        # If `nan_mask` is lazy loaded, then resulting `backscatter_chan` will be lazy loaded.
+        backscatter_chan = xr.where(nan_mask, 0.0 + 0j, backscatter_chan)
+
+        # Extract transmit values
         tx = chirp[str(chan.values)]
+
+        # Compute complex conjugate of transmit values and reverse order of elements
         replica = np.flipud(np.conj(tx))
 
+        # Apply convolve on backscatter (along range sample dimension) and replica
         pc = xr.apply_ufunc(
             lambda m: (signal.convolve(m, replica, mode="full")[replica.size - 1 :]),
             backscatter_chan,
@@ -291,6 +304,11 @@ def compress_pulse(backscatter: xr.DataArray, chirp: Dict) -> xr.DataArray:
             output_dtypes=[np.complex64],
         ).compute()
 
+        # Restore NaN values in the resulting array.
+        # Computing of `nan_mask` here is necessary in the case when `nan_mask` is lazy loaded
+        # or else the resulting `pc` will also be lazy loaded.
+        pc = xr.where(nan_mask.compute(), np.nan, pc)
+
         pc_all.append(pc)
 
     pc_all = xr.concat(pc_all, dim="channel")
diff --git a/echopype/calibrate/env_params.py b/echopype/calibrate/env_params.py
index 49be4d12c..d19bb18fb 100644
--- a/echopype/calibrate/env_params.py
+++ b/echopype/calibrate/env_params.py
@@ -50,10 +50,13 @@ def harmonize_env_param_time(
         if "time1" not in p.coords:
             return p
 
-        # If there's only 1 time1 value,
-        # or if after dropping NaN there's only 1 time1 value
-        if p["time1"].size == 1 or p.dropna(dim="time1").size == 1:
-            return p.dropna(dim="time1").squeeze(dim="time1").drop("time1")
+        # If there's only 1 time1 value:
+        if p["time1"].size == 1:
+            return p.squeeze(dim="time1").drop_vars("time1")
+
+        # If after dropping NaN along time1 dimension there's only 1 time1 value:
+        if p.dropna(dim="time1").size == 1:
+            return p.dropna(dim="time1").squeeze(dim="time1").drop_vars("time1")
 
         if ping_time is None:
             raise ValueError(


Test patch:
diff --git a/echopype/tests/calibrate/test_calibrate_ek80.py b/echopype/tests/calibrate/test_calibrate_ek80.py
index 53722a2b2..a4b4011d7 100644
--- a/echopype/tests/calibrate/test_calibrate_ek80.py
+++ b/echopype/tests/calibrate/test_calibrate_ek80.py
@@ -177,7 +177,6 @@ def test_ek80_BB_range(ek80_cal_path, ek80_ext_path):
     pyel_vals = pyel_BB_p_data["range"]
     assert np.allclose(pyel_vals, ep_vals)
 
-
 @pytest.mark.parametrize(
     ("raw_data_path,raw_file_name,pyecholab_data_path,pyecholab_file_path, dask_array"),
     [
@@ -253,7 +252,6 @@ def test_ek80_BB_power_from_complex(
     )
     assert np.allclose(pyel_vals[idx_to_cmp], ep_vals[idx_to_cmp])
 
-
 @pytest.mark.parametrize(
     ("raw_data_path,raw_file_name,pyecholab_data_path,pyecholab_file_path, dask_array"),
     [
@@ -371,10 +369,9 @@ def test_ek80_BB_power_echoview(ek80_path):
     # Below only compare the first ping
     ev_vals = df_real.values[:, :]
     ep_vals = pc_mean.values.real[:, :]
-    assert np.allclose(ev_vals[:, 69:8284], ep_vals[:, 69:], atol=1e-4)
-    assert np.allclose(ev_vals[:, 90:8284], ep_vals[:, 90:], atol=1e-5)
+    assert np.allclose(ev_vals[:, 69:], ep_vals[:, 69:], atol=1e-4)
+    assert np.allclose(ev_vals[:, 90:], ep_vals[:, 90:], atol=1e-5)
 
-    
 def test_ek80_CW_complex_Sv_receiver_sampling_freq(ek80_path):
     ek80_raw_path = str(ek80_path.joinpath("D20230804-T083032.raw"))
     ed = ep.open_raw(ek80_raw_path, sonar_model="EK80")
@@ -390,3 +387,67 @@ def test_ek80_CW_complex_Sv_receiver_sampling_freq(ek80_path):
     # receiver_sampling_frequency substituted with default value in compute_Sv
     assert ds_Sv["receiver_sampling_frequency"] is not None
     assert np.allclose(ds_Sv["receiver_sampling_frequency"].data, 1500000)
+
+
+@pytest.mark.parametrize(
+    ("raw_data_path,target_channel_ping_pattern"),
+    [
+        (
+            "NYOS2105-D20210525-T213648.raw",
+            # (3,4), 2, 1, (3,4), 2, 1, etc...
+            np.array(
+                [
+                    [np.nan, np.nan,  1., np.nan, np.nan,  1., np.nan, np.nan,  1., np.nan],
+                    [np.nan,  1., np.nan, np.nan,  1., np.nan, np.nan,  1., np.nan, np.nan],
+                    [ 1., np.nan, np.nan,  1., np.nan, np.nan,  1., np.nan, np.nan,  1.],
+                    [ 1., np.nan, np.nan,  1., np.nan, np.nan,  1., np.nan, np.nan,  1.]
+                ]
+            ),
+        ),
+        (
+            "DRIX08-D20231003-T120051.raw",
+            # 2, 1, 2, 1, 2, 1, etc...
+            np.array(
+                [
+                    [np.nan,  1., np.nan, 1,  np.nan, 1, np.nan,  1., np.nan, 1],
+                    [1, np.nan,  1., np.nan, 1,  np.nan, 1, np.nan,  1., np.nan],
+                ]
+            ),
+        ),
+    ],
+)
+@pytest.mark.integration
+def test_ek80_BB_complex_multiplex_NaNs_and_non_NaNs(raw_data_path, target_channel_ping_pattern):
+    # Extract bb complex multiplex EK80 data
+    ed = ep.open_raw(f"echopype/test_data/ek80_bb_complex_multiplex/{raw_data_path}", sonar_model="EK80")
+
+    # Compute Sv
+    ds_Sv = ep.calibrate.compute_Sv(ed,waveform_mode='BB',encode_mode='complex')
+    
+    # Extract mask for both real backscatter and calibrated Sv:
+    
+    # EchoData Real (component) Backscatter mask should be 1 if all values along 
+    # a specific range sample and beam dimension are non-NaN; else NaN.
+    ed_backscatter_r_mask = xr.where(
+        ~np.isnan(ed["Sonar/Beam_group1"]["backscatter_r"]
+    ).all(dim=['range_sample', "beam"]), 1, np.nan)
+    # Calibrated Sv mask should be 1 if all values along 
+    # a specific range sample are non-NaN; else NaN.
+    calibrated_Sv_mask = xr.where(
+        ~np.isnan(ds_Sv["Sv"]
+    ).all(dim=['range_sample']), 1, np.nan)
+
+    # Check that they are equal
+    assert np.array_equal(
+        ed_backscatter_r_mask.data,
+        calibrated_Sv_mask.data,
+        equal_nan=True
+    )
+
+    # Check that a slice of calibrated Sv mask is equal to the following array
+    # with the target channel ping pattern:
+    assert np.array_equal(
+        calibrated_Sv_mask.isel(ping_time=slice(0,10)).data,
+        target_channel_ping_pattern,
+        equal_nan=True
+    )
diff --git a/echopype/tests/calibrate/test_ek80_complex.py b/echopype/tests/calibrate/test_ek80_complex.py
index 66abe8539..16a3a95a5 100644
--- a/echopype/tests/calibrate/test_ek80_complex.py
+++ b/echopype/tests/calibrate/test_ek80_complex.py
@@ -70,4 +70,4 @@ def test_get_vend_filter_EK80(ch_num, filter_len, has_nan):
 
             assert sel_vend[var_df].values == get_vend_filter_EK80(
                 vend, channel_id=ch, filter_name=filter_name, param_type="decimation"
-            )
\ No newline at end of file
+            )
diff --git a/echopype/tests/calibrate/test_env_params.py b/echopype/tests/calibrate/test_env_params.py
index 205376032..98e31d74b 100644
--- a/echopype/tests/calibrate/test_env_params.py
+++ b/echopype/tests/calibrate/test_env_params.py
@@ -3,6 +3,7 @@
 import dask.array
 import numpy as np
 import xarray as xr
+import pandas as pd
 
 import echopype as ep
 from echopype.calibrate.env_params import (
@@ -132,6 +133,17 @@ def test_harmonize_env_param_time():
     # .all computes dask array under the hood
     assert (p_new.data == [0.5, 2880.5]).all()
 
+@pytest.mark.unit
+def test_harmonize_env_param_time_only_one_non_NaN_along_time1():
+    # Create data array with time1 dimension with only 1 non-NaN value
+    data = np.array([1, np.nan, np.nan])
+    time = pd.date_range(start='2024-01-01', periods=3, freq='D')
+    da = xr.DataArray(data, dims='time1', coords={'time1': time})
+
+    # Check that output da has only 1 value and no dimension
+    output_da = harmonize_env_param_time(da, None)
+    assert output_da == 1, "Output data array should just be 1"
+    assert 'time1' not in output_da.dims, "```harmonize_env_param_time``` should have dropped 'time1' dimension"
 
 @pytest.mark.parametrize(
     ("user_dict", "channel", "out_dict"),