Problem reading RDI Pinnacle 45 burst data with read_rdi

[stinawahlgren]

Describe the bug:

I get the following error message when I try to read a pd0-file from a RDI Pinnacle 45 instrument using read_rdi with mhkit version 0.9.0:

ZeroDivisionError                         Traceback (most recent call last)
Cell In[3], line 1
----> 1 data = read_rdi(file)

File [...]/MHKiT-python/mhkit/dolfyn/io/rdi.py:75, in read_rdi(filename, userdata, nens, debug, vmdas_search, winriver, search_num, **kwargs)
     66 # Reads into a dictionary of dictionaries using netcdf naming conventions
     67 # Should be easier to debug
     68 rdr = _RDIReader(
     69     filename,
     70     debug=debug,
   (...)
     73     search_num=search_num,
     74 )
---> 75 datNB, datBB = rdr.load_data(nens=nens)
     77 dats = [dat for dat in [datNB, datBB] if dat is not None]
     79 # Read in userdata

File [...]/MHKiT-python/mhkit/dolfyn/io/rdi.py:355, in _RDIReader.load_data(self, nens)
    353 for dat, cfg in zip(datl, cfgl):
    354     dat, cfg = self.cleanup(dat, cfg)
--> 355     dat = self.finalize(dat, cfg)
    356     if "vel_bt" in dat["data_vars"]:
    357         cfg["rotate_vars"].append("vel_bt")

File [...]/MHKiT-python/mhkit/dolfyn/io/rdi.py:1047, in _RDIReader.finalize(self, dat, cfg)
   1045     cfg["fs"] = round(1 / np.median(np.diff(dat["coords"]["time"])), 2)
   1046 else:
-> 1047     cfg["fs"] = 1 / (cfg["sec_between_ping_groups"] * cfg["pings_per_ensemble"])
   1049 # Save configuration data as attributes
   1050 dat["attrs"] = cfg

ZeroDivisionError: float division by zero

The data is collected in burst mode and with both water column data and bottom track. The reason for the division-by-zero error is that cfg["sec_between_ping_groups"] is 0.

If I comment out line 1047 in mhkit/dolfyn/io/rdi.py (or use line 1045 instead) the function works and the file is read correctly.

Desktop:

• MHKiT Version: 0.9.0

[simmsa]

@stinawahlgren, thank you for reporting! I've identified the likely cause of this division by zero error.

The issue appears to be that either sec_between_ping_groups or pings_per_ensemble is reading as 0 from your RDI file, possibly due to burst mode configuration differences.

To help us fix this could you provide the RDI file that causes this error (or a subset file that is smaller)? Since GitHub has file size limits, feel free to share via Google Drive or similar cloud storage.

I need to ask @jmcvey3, the primary developer of the ADCP/dolfyn code, for help in determining the root cause and a developing fix for this bug.

[simmsa]

@jmcvey3, here's some context for this bug:

The issue is in this sampling rate calculation:

MHKiT-Python/mhkit/dolfyn/io/rdi.py

Lines 1040 to 1047 in 6a474e9

	# Need to figure out how to differentiate burst mode from averaging mode
	if (
	("source_program" in cfg)
	and (cfg["source_program"].lower() in ["vmdas", "winriver", "winriver2"])
	) or ("sentinelv" in cfg["inst_model"].lower()):
	cfg["fs"] = round(1 / np.median(np.diff(dat["coords"]["time"])), 2)
	else:
	cfg["fs"] = 1 / (cfg["sec_between_ping_groups"] * cfg["pings_per_ensemble"])

Two potential causes:

1. sec_between_ping_groups reading as 0:

MHKiT-Python/mhkit/dolfyn/io/rdi_defs.py

Lines 463 to 465 in 6a474e9

	cfg["sec_between_ping_groups"] = round(
	float(np.sum(np.array(fd.read_ui8(3)) * [60.0, 1.0, 0.01])), 3
	)

Could RDI burst mode output zero for all time values?

2. pings_per_ensemble reading as 0:

MHKiT-Python/mhkit/dolfyn/io/rdi_defs.py

Line 449 in 6a474e9

cfg["pings_per_ensemble"] = fd.read_ui16(1)

Is there some reason the instrument would output 0 here?

Questions:

• Do you have any ideas/clues on the 0 outputs above. It seems possible that there is an RDI specific edge case with burst mode, but these seem hard to troubleshoot without having actual files to work with?
• Would falling back to time-based fs calculation be technically correct for burst mode if the ping-based method fails?

[stinawahlgren]

@simmsa I realized that I made a typo in my bug report, it's corrected now. It is indeed cfg["sec_between_ping_groups"] that is 0, while cfg["pings_per_ensemble"] is 1. We were running the instrument in continuous burst mode (essentially pinging as often as possible for the whole deployment).

[stinawahlgren]

Some thoughts as a Pinnacle user (feel free to ignore if it is not relevant):

If I understand correctly, setting sec_between_ping_groups to 0 while configurating the instrument will make it ping as often as possible. In our case, where we used both bottom track and water column data, this actually causes a variable ping rate, so using the median to compute cfg["fs"] (as is done on line 1045) would perhaps be a bit misleading? On the other hand, cfg["fs"] is not well defined in that case anyway.

To summarize:

• There is a use case for setting the config variable behind sec_between_ping_groups to zero when configurating the instrument.
• There are cases where the actual time between pings is irregular

[simmsa]

@stinawahlgren, no worries on the typo. And the practical perspective as a Pinnacle user is valuable for us to hear and understand.

Based on the specified instrument settings the data is being read by dolfyn correctly, but there is an edge case in burst mode where pings_per_ensemble outputs a zero value, which likely indicates that burst mode is active, and also likely indicates that the underlying sample time delta may vary (or is unlikely to be consistent) which causes any calculation of an average sample rate to be incorrect. Does that seem like an accurate summary here?

To fix, MHKiT/dolfyn should understand and handle this edge case gracefully. How exactly we handle this edge case is the next decision we need to make.

@jmcvey3, do you have any guidance here?

I think a solution is to have a fallback value for fs, either np.nan or 0 and warn the user instead of the undesired raising of the ZeroDivisionError. This is likely to cause some issues during data processing, and it possible that there is a better solution.

[stinawahlgren]

@simmsa Yep, I think that's an accurate summary.

[jmcvey3]

@stinawahlgren Thanks for bringing this up! RDI doesn't explicitly encode a sampling rate because "reasons", so when you have an instrument in burst mode they don't save it clearly, if at all.

@simmsa For better or worse, the best way I've found of doing this is to just calculate the difference between timesteps and average them.

[simmsa]

@jmcvey3 sounds good. I updated the corresponding PR #409 accordingly. Thanks for pointing us in the right direction.

[jmcvey3]

Closing now that #409 is merged