diff --git a/examples/data/dolfyn/dual_profile.ad2cp b/examples/data/dolfyn/dual_profile.ad2cp new file mode 100644 index 000000000..7839b3285 Binary files /dev/null and b/examples/data/dolfyn/dual_profile.ad2cp differ diff --git a/examples/data/dolfyn/test_data/BenchFile01.nc b/examples/data/dolfyn/test_data/BenchFile01.nc index 57a73e44d..3b2af8fc4 100644 Binary files a/examples/data/dolfyn/test_data/BenchFile01.nc and b/examples/data/dolfyn/test_data/BenchFile01.nc differ diff --git a/examples/data/dolfyn/test_data/BenchFile01_avg.nc b/examples/data/dolfyn/test_data/BenchFile01_avg.nc index 0fa4e236b..24d488138 100644 Binary files a/examples/data/dolfyn/test_data/BenchFile01_avg.nc and b/examples/data/dolfyn/test_data/BenchFile01_avg.nc differ diff --git a/examples/data/dolfyn/test_data/BenchFile01_rotate_beam2inst.nc b/examples/data/dolfyn/test_data/BenchFile01_rotate_beam2inst.nc index d36b432d1..2004de5f4 100644 Binary files a/examples/data/dolfyn/test_data/BenchFile01_rotate_beam2inst.nc and b/examples/data/dolfyn/test_data/BenchFile01_rotate_beam2inst.nc differ diff --git a/examples/data/dolfyn/test_data/BenchFile01_rotate_earth2principal.nc b/examples/data/dolfyn/test_data/BenchFile01_rotate_earth2principal.nc index 2a1ff37df..a71cbbbdd 100644 Binary files a/examples/data/dolfyn/test_data/BenchFile01_rotate_earth2principal.nc and b/examples/data/dolfyn/test_data/BenchFile01_rotate_earth2principal.nc differ diff --git a/examples/data/dolfyn/test_data/BenchFile01_rotate_inst2earth.nc b/examples/data/dolfyn/test_data/BenchFile01_rotate_inst2earth.nc index 194dee8a2..bbfeaf37e 100644 Binary files a/examples/data/dolfyn/test_data/BenchFile01_rotate_inst2earth.nc and b/examples/data/dolfyn/test_data/BenchFile01_rotate_inst2earth.nc differ diff --git a/examples/data/dolfyn/test_data/dual_profile.nc b/examples/data/dolfyn/test_data/dual_profile.nc new file mode 100644 index 000000000..8e3a98dbe Binary files /dev/null and b/examples/data/dolfyn/test_data/dual_profile.nc differ diff --git a/mhkit/dolfyn/io/base.py b/mhkit/dolfyn/io/base.py index a4414cbe7..545035cdb 100644 --- a/mhkit/dolfyn/io/base.py +++ b/mhkit/dolfyn/io/base.py @@ -113,69 +113,65 @@ def _handle_nan(data): def _create_dataset(data): - """Creates an xarray dataset from dictionary created from binary + """ + Creates an xarray dataset from dictionary created from binary readers. Direction 'dir' coordinates are set in `set_coords` """ - ds = xr.Dataset() - tag = ["_avg", "_b5", "_echo", "_bt", "_gps", "_altraw", "_sl"] - FoR = {} - try: - beams = data["attrs"]["n_beams"] - except: + tag = ["_avg", "_b5", "_echo", "_bt", "_gps", "_altraw", "_altraw_avg", "_sl"] + + ds_dict = {} + for key in data["coords"]: + ds_dict[key] = {"dims": (key), "data": data["coords"][key]} + + # Set various coordinate frames + if "n_beams_avg" in data["attrs"]: beams = data["attrs"]["n_beams_avg"] + else: + beams = data["attrs"]["n_beams"] n_beams = max(min(beams, 4), 3) beams = np.arange(1, n_beams + 1, dtype=np.int32) - FoR["beam"] = xr.DataArray( - beams, - dims=["beam"], - name="beam", - attrs={"units": "1", "long_name": "Beam Reference Frame"}, - ) - FoR["dir"] = xr.DataArray( - beams, - dims=["dir"], - name="dir", - attrs={"units": "1", "long_name": "Reference Frame"}, - ) + ds_dict["beam"] = {"dims": ("beam"), "data": beams} + ds_dict["dir"] = {"dims": ("dir"), "data": beams} + data["units"].update({"beam": "1", "dir": "1"}) + data["long_name"].update({"beam": "Beam Reference Frame", "dir": "Reference Frame"}) + + # Iterate through data variables and add them to new dictionary for key in data["data_vars"]: # orientation matrices if "mat" in key: if "inst" in key: # beam2inst & inst2head orientation matrices - ds[key] = xr.DataArray( - data["data_vars"][key], - coords={"x1": beams, "x2": beams}, - dims=["x1", "x2"], - attrs={"units": "1", "long_name": "Rotation Matrix"}, - ) + if "x1" not in ds_dict: + ds_dict["x1"] = {"dims": ("x1"), "data": beams} + ds_dict["x2"] = {"dims": ("x2"), "data": beams} + + ds_dict[key] = {"dims": ("x1", "x2"), "data": data["data_vars"][key]} + data["units"].update({key: "1"}) + data["long_name"].update({key: "Rotation Matrix"}) + elif "orientmat" in key: # earth2inst orientation matrix if any(val in key for val in tag): tg = "_" + key.rsplit("_")[-1] else: tg = "" - earth = xr.DataArray( - ["E", "N", "U"], - dims=["earth"], - name="earth", - attrs={"units": "1", "long_name": "Earth Reference Frame"}, - ) - inst = xr.DataArray( - ["X", "Y", "Z"], - dims=["inst"], - name="inst", - attrs={"units": "1", "long_name": "Instrument Reference Frame"}, + + ds_dict["earth"] = {"dims": ("earth"), "data": ["E", "N", "U"]} + ds_dict["inst"] = {"dims": ("inst"), "data": ["X", "Y", "Z"]} + ds_dict[key] = { + "dims": ("earth", "inst", "time" + tg), + "data": data["data_vars"][key], + } + data["units"].update( + {"earth": "1", "inst": "1", key: data["units"]["orientmat"]} ) - time = data["coords"]["time" + tg] - ds[key] = xr.DataArray( - data["data_vars"][key], - coords={"earth": earth, "inst": inst, "time" + tg: time}, - dims=["earth", "inst", "time" + tg], - attrs={ - "units": data["units"]["orientmat"], - "long_name": data["long_name"]["orientmat"], - }, + data["long_name"].update( + { + "earth": "Earth Reference Frame", + "inst": "Instrument Reference Frame", + key: data["long_name"]["orientmat"], + } ) # quaternion units never change @@ -184,67 +180,43 @@ def _create_dataset(data): tg = "_" + key.rsplit("_")[-1] else: tg = "" - q = xr.DataArray( - ["w", "x", "y", "z"], - dims=["q"], - name="q", - attrs={"units": "1", "long_name": "Quaternion Vector Components"}, - ) - time = data["coords"]["time" + tg] - ds[key] = xr.DataArray( - data["data_vars"][key], - coords={"q": q, "time" + tg: time}, - dims=["q", "time" + tg], - attrs={ - "units": data["units"]["quaternions"], - "long_name": data["long_name"]["quaternions"], - }, - ) - else: - # Assign each variable to a dataArray - ds[key] = xr.DataArray(data["data_vars"][key]) - # Assign metadata to each dataArray - for md in ["units", "long_name", "standard_name"]: - if key in data[md]: - ds[key].attrs[md] = data[md][key] - try: # make sure ones with tags get units - tg = "_" + key.rsplit("_")[-1] - if any(val in key for val in tag): - ds[key].attrs[md] = data[md][key[: -len(tg)]] - except: - pass - # Fill in dimensions and coordinates for each dataArray + if "q" not in ds_dict: + ds_dict["q"] = {"dims": ("q"), "data": ["w", "x", "y", "z"]} + data["units"].update({"q": "1"}) + data["long_name"].update({"q": "Quaternion Vector Components"}) + + ds_dict[key] = {"dims": ("q", "time" + tg), "data": data["data_vars"][key]} + data["units"].update({key: data["units"]["quaternions"]}) + data["long_name"].update({key: data["long_name"]["quaternions"]}) + + else: shp = data["data_vars"][key].shape - l = len(shp) - if l == 1: # 1D variables - if any(val in key for val in tag): + if len(shp) == 1: # 1D variables + if "_altraw_avg" in key: + tg = "_altraw_avg" + elif any(val in key for val in tag): tg = "_" + key.rsplit("_")[-1] else: tg = "" - ds[key] = ds[key].rename({"dim_0": "time" + tg}) - ds[key] = ds[key].assign_coords( - {"time" + tg: data["coords"]["time" + tg]} - ) + ds_dict[key] = {"dims": ("time" + tg), "data": data["data_vars"][key]} - elif l == 2: # 2D variables + elif len(shp) == 2: # 2D variables if key == "echo": - ds[key] = ds[key].rename( - {"dim_0": "range_echo", "dim_1": "time_echo"} - ) - ds[key] = ds[key].assign_coords( - { - "range_echo": data["coords"]["range_echo"], - "time_echo": data["coords"]["time_echo"], - } - ) - elif key == "samp_altraw": # raw altimeter samples - ds[key] = ds[key].rename( - {"dim_0": "n_altraw", "dim_1": "time_altraw"} - ) - ds[key] = ds[key].assign_coords( - {"time_altraw": data["coords"]["time_altraw"]} - ) + ds_dict[key] = { + "dims": ("range_echo", "time_echo"), + "data": data["data_vars"][key], + } + elif key == "samp_altraw": + ds_dict[key] = { + "dims": ("n_altraw", "time_altraw"), + "data": data["data_vars"][key], + } + elif key == "samp_altraw_avg": + ds_dict[key] = { + "dims": ("n_altraw_avg", "time_altraw_avg"), + "data": data["data_vars"][key], + } # ADV/ADCP instrument vector data, bottom tracking elif shp[0] == n_beams and not any(val in key for val in tag[:3]): @@ -259,10 +231,11 @@ def _create_dataset(data): dim0 = "beam" else: dim0 = "dir" - ds[key] = ds[key].rename({"dim_0": dim0, "dim_1": "time" + tg}) - ds[key] = ds[key].assign_coords( - {dim0: FoR[dim0], "time" + tg: data["coords"]["time" + tg]} - ) + ds_dict[key] = { + "dims": (dim0, "time" + tg), + "data": data["data_vars"][key], + } + # ADCP IMU data elif shp[0] == 3: if not any(val in key for val in tag): @@ -270,20 +243,24 @@ def _create_dataset(data): else: tg = [val for val in tag if val in key] tg = tg[0] - dirIMU = xr.DataArray( - [1, 2, 3], - dims=["dirIMU"], - name="dirIMU", - attrs={"units": "1", "long_name": "Reference Frame"}, - ) - ds[key] = ds[key].rename({"dim_0": "dirIMU", "dim_1": "time" + tg}) - ds[key] = ds[key].assign_coords( - {"dirIMU": dirIMU, "time" + tg: data["coords"]["time" + tg]} - ) - ds[key].attrs["coverage_content_type"] = "physicalMeasurement" + if "dirIMU" not in ds_dict: + ds_dict["dirIMU"] = {"dims": ("dirIMU"), "data": [1, 2, 3]} + data["units"].update({"dirIMU": "1"}) + data["long_name"].update({"dirIMU": "Reference Frame"}) + + ds_dict[key] = { + "dims": ("dirIMU", "time" + tg), + "data": data["data_vars"][key], + } - elif l == 3: # 3D variables + elif "b5" in tg: + ds_dict[key] = { + "dims": ("range_b5", "time_b5"), + "data": data["data_vars"][key], + } + + elif len(shp) == 3: # 3D variables if "vel" in key: dim0 = "dir" else: # amp, corr, prcnt_gd, status @@ -294,43 +271,43 @@ def _create_dataset(data): tg = "_avg" else: tg = "" - ds[key] = ds[key].rename( - {"dim_0": dim0, "dim_1": "range" + tg, "dim_2": "time" + tg} - ) - ds[key] = ds[key].assign_coords( - { - dim0: FoR[dim0], - "range" + tg: data["coords"]["range" + tg], - "time" + tg: data["coords"]["time" + tg], - } - ) + ds_dict[key] = { + "dims": (dim0, "range" + tg, "time" + tg), + "data": data["data_vars"][key], + } + elif "b5" in key: - # xarray can't handle coords of length 1 - ds[key] = ds[key][0] - ds[key] = ds[key].rename({"dim_1": "range_b5", "dim_2": "time_b5"}) - ds[key] = ds[key].assign_coords( - { - "range_b5": data["coords"]["range_b5"], - "time_b5": data["coords"]["time_b5"], - } - ) + # "vel_b5" sometimes stored as (1, range_b5, time_b5) + ds_dict[key] = { + "dims": ("range_b5", "time_b5"), + "data": data["data_vars"][key][0], + } elif "sl" in key: - ds[key] = ds[key].rename( - {"dim_0": dim0, "dim_1": "range_sl", "dim_2": "time"} - ) - ds[key] = ds[key].assign_coords( - { - "range_sl": data["coords"]["range_sl"], - "time": data["coords"]["time"], - } - ) + ds_dict[key] = { + "dims": (dim0, "range_sl", "time"), + "data": data["data_vars"][key], + } else: - ds = ds.drop_vars(key) warnings.warn(f"Variable not included in dataset: {key}") + # Create dataset + ds = xr.Dataset.from_dict(ds_dict) + + # Assign data array attributes + for key in ds.variables: + for md in ["units", "long_name", "standard_name"]: + if key in data[md]: + ds[key].attrs[md] = data[md][key] + if len(ds[key].shape) > 1: ds[key].attrs["coverage_content_type"] = "physicalMeasurement" + try: # make sure ones with tags get units + tg = "_" + key.rsplit("_")[-1] + if any(val in key for val in tag): + ds[key].attrs[md] = data[md][key[: -len(tg)]] + except: + pass - # coordinate attributes + # Assign coordinate attributes for ky in ds.dims: ds[ky].attrs["coverage_content_type"] = "coordinate" r_list = [r for r in ds.coords if "range" in r] @@ -344,7 +321,7 @@ def _create_dataset(data): ds[ky].attrs["long_name"] = "Time" ds[ky].attrs["standard_name"] = "time" - # dataset metadata + # Set dataset metadata ds.attrs = data["attrs"] return ds diff --git a/mhkit/dolfyn/io/nortek2.py b/mhkit/dolfyn/io/nortek2.py index c29035a92..fd984a817 100644 --- a/mhkit/dolfyn/io/nortek2.py +++ b/mhkit/dolfyn/io/nortek2.py @@ -15,7 +15,13 @@ def read_signature( - filename, userdata=True, nens=None, rebuild_index=False, debug=False, **kwargs + filename, + userdata=True, + nens=None, + rebuild_index=False, + debug=False, + dual_profile=False, + **kwargs ): """ Read a Nortek Signature (.ad2cp) datafile @@ -34,6 +40,8 @@ def read_signature( Default = False debug : bool Logs debugger ouput if true. Default = False + dual_profile : bool + Set to true if instrument is running multiple profiles. Default = False Returns ------- @@ -68,9 +76,13 @@ def read_signature( userdata = _find_userdata(filename, userdata) - rdr = _Ad2cpReader(filename, rebuild_index=rebuild_index, debug=debug) + rdr = _Ad2cpReader( + filename, rebuild_index=rebuild_index, debug=debug, dual_profile=dual_profile + ) d = rdr.readfile(nens[0], nens[1]) rdr.sci_data(d) + if rdr._dp: + _clean_dp_skips(d) out = _reorg(d) _reduce(out) @@ -120,12 +132,22 @@ def read_signature( logging.root.removeHandler(handler) handler.close() - return ds + # Return two datasets if dual profile + if rdr._dp: + return split_dp_datasets(ds) + else: + return ds class _Ad2cpReader: def __init__( - self, fname, endian=None, bufsize=None, rebuild_index=False, debug=False + self, + fname, + endian=None, + bufsize=None, + rebuild_index=False, + debug=False, + dual_profile=False, ): self.fname = fname self.debug = debug @@ -133,8 +155,13 @@ def __init__( self.f.seek(0, 2) # Seek to end self._eof = self.f.tell() self.start_pos = self._check_header() - self._index = lib.get_index( - fname, pos=self.start_pos, eof=self._eof, reload=rebuild_index, debug=debug + self._index, self._dp = lib.get_index( + fname, + pos=self.start_pos, + eof=self._eof, + rebuild=rebuild_index, + debug=debug, + dp=dual_profile, ) self._reopen(bufsize) self.filehead_config = self._read_filehead_config_string() @@ -257,18 +284,23 @@ def init_data(self, ens_start, ens_stop): outdat = {} nens = int(ens_stop - ens_start) - # ID 26 usually only recorded in first ensemble - n26 = ( - (self._index["ID"] == 26) - & (self._index["ens"] >= ens_start) - & (self._index["ens"] < ens_stop) - ).sum() - if not n26 and 26 in self._burst_readers: + # ID 26 and 31 recorded infrequently + def n_id(id): + return ( + (self._index["ID"] == id) + & (self._index["ens"] >= ens_start) + & (self._index["ens"] < ens_stop) + ).sum() + + n_altraw = {26: n_id(26), 31: n_id(31)} + if not n_altraw[26] and 26 in self._burst_readers: self._burst_readers.pop(26) + if not n_altraw[31] and 31 in self._burst_readers: + self._burst_readers.pop(31) for ky in self._burst_readers: - if ky == 26: - n = n26 + if (ky == 26) or (ky == 31): + n = n_altraw[ky] ens = np.zeros(n, dtype="uint32") else: ens = np.arange(ens_start, ens_stop).astype("uint32") @@ -310,7 +342,7 @@ def readfile(self, ens_start=0, ens_stop=None): outdat["filehead_config"] = self.filehead_config print("Reading file %s ..." % self.fname) c = 0 - c26 = 0 + c_altraw = {26: 0, 31: 0} self.f.seek(self._ens_pos[ens_start], 0) while True: try: @@ -322,9 +354,9 @@ def readfile(self, ens_start=0, ens_stop=None): # "avg data record" (vel_avg + ast_avg), "bottom track data record" (bt), # "interleaved burst data record" (vel_b5), "echosounder record" (echo) self._read_burst(id, outdat[id], c) - elif id in [26]: - # "burst altimeter raw record" (alt_raw) - recorded on nens==0 - rdr = self._burst_readers[26] + elif id in [26, 31]: + # "burst altimeter raw record" (_altraw), "avg altimeter raw record" (_altraw_avg) + rdr = self._burst_readers[id] if not hasattr(rdr, "_nsamp_index"): first_pass = True tmp_idx = rdr._nsamp_index = rdr._names.index("nsamp_alt") @@ -350,8 +382,8 @@ def readfile(self, ens_start=0, ens_stop=None): "<" + "{}H".format(int(rdr.nbyte // 2)) ) # Initialize the array - outdat[26]["samp_alt"] = defs._nans( - [rdr._N[tmp_idx], len(outdat[26]["samp_alt"])], dtype=np.uint16 + outdat[id]["samp_alt"] = defs._nans( + [rdr._N[tmp_idx], len(outdat[id]["samp_alt"])], dtype=np.uint16 ) else: if sz != rdr._N[tmp_idx]: @@ -359,12 +391,12 @@ def readfile(self, ens_start=0, ens_stop=None): "The number of samples in this 'Altimeter Raw' " "burst is different from prior bursts." ) - self._read_burst(id, outdat[id], c26) - outdat[id]["ensemble"][c26] = c - c26 += 1 + self._read_burst(id, outdat[id], c_altraw[id]) + outdat[id]["ensemble"][c_altraw[id]] = c + c_altraw[id] += 1 - elif id in [27, 29, 30, 31, 35, 36]: # unknown how to handle - # "bottom track record", DVL, "altimeter record", "avg altimeter raw record", + elif id in [27, 29, 30, 35, 36]: # unknown how to handle + # "bottom track record", DVL, "altimeter record", # "raw echosounder data record", "raw echosounder transmit data record" if self.debug: logging.debug("Skipped ID: 0x{:02X} ({:02d})\n".format(id, id)) @@ -430,18 +462,37 @@ def sci_data(self, dat): "float32" ) - def __exit__( - self, - type, - value, - trace, - ): - self.f.close() - def __enter__( - self, - ): - return self +def _altraw_reorg(outdat, tag=""): + """Submethod for `_reorg` particular to raw altimeter pings (ID 26 and 31)""" + for ky in list(outdat["data_vars"]): + if ky.endswith("raw" + tag) and not ky.endswith("_altraw" + tag): + outdat["data_vars"].pop(ky) + outdat["coords"]["time_altraw" + tag] = outdat["coords"].pop("timeraw" + tag) + # convert "signed fractional" to float + outdat["data_vars"]["samp_altraw" + tag] = ( + outdat["data_vars"]["samp_altraw" + tag].astype("float32") / 2**8 + ) + + # Read altimeter status + outdat["data_vars"].pop("status_altraw" + tag) + status_alt = lib._alt_status2data(outdat["data_vars"]["status_alt" + tag]) + for ky in status_alt: + outdat["attrs"][ky + tag] = lib._collapse( + status_alt[ky].astype("uint8"), name=ky + ) + outdat["data_vars"].pop("status_alt" + tag) + + # Power level index + power = {0: "high", 1: "med-high", 2: "med-low", 3: "low"} + outdat["attrs"]["power_level_alt" + tag] = power[ + outdat["attrs"].pop("power_level_idx_alt" + tag) + ] + + # Other attrs + for ky in list(outdat["attrs"]): + if ky.endswith("raw" + tag): + outdat["attrs"][ky.split("raw")[0] + "_alt" + tag] = outdat["attrs"].pop(ky) def _reorg(dat): @@ -473,6 +524,7 @@ def _reorg(dat): (24, "_b5"), (26, "raw"), (28, "_echo"), + (31, "raw_avg"), ]: if id in [24, 26]: collapse_exclude = [0] @@ -503,7 +555,7 @@ def _reorg(dat): lib._collapse( dnow["beam_config"], exclude=collapse_exclude, name="beam_config" ), - 21, + 21, # always 21 here ) cfg["n_cells" + tag] = tmp["n_cells"] cfg["coord_sys_axes" + tag] = tmp["cy"] @@ -575,26 +627,9 @@ def _reorg(dat): # Move 'altimeter raw' data to its own down-sampled structure if 26 in dat: - for ky in list(outdat["data_vars"]): - if ky.endswith("raw") and not ky.endswith("_altraw"): - outdat["data_vars"].pop(ky) - outdat["coords"]["time_altraw"] = outdat["coords"].pop("timeraw") - outdat["data_vars"]["samp_altraw"] = ( - outdat["data_vars"]["samp_altraw"].astype("float32") / 2**8 - ) # convert "signed fractional" to float - - # Read altimeter status - outdat["data_vars"].pop("status_altraw") - status_alt = lib._alt_status2data(outdat["data_vars"]["status_alt"]) - for ky in status_alt: - outdat["attrs"][ky] = lib._collapse(status_alt[ky].astype("uint8"), name=ky) - outdat["data_vars"].pop("status_alt") - - # Power level index - power = {0: "high", 1: "med-high", 2: "med-low", 3: "low"} - outdat["attrs"]["power_level_alt"] = power[ - outdat["attrs"].pop("power_level_idx_alt") - ] + _altraw_reorg(outdat) + if 31 in dat: + _altraw_reorg(outdat, tag="_avg") # Read status data status0_vars = [x for x in outdat["data_vars"] if "status0" in x] @@ -641,7 +676,7 @@ def _reorg(dat): for ky in status0_data: outdat["attrs"][ky] = lib._collapse(status0_data[ky].astype("uint8"), name=ky) - # Remove status0 variables - keep status variables as they useful for finding missing pings + # Remove status0 variables - keep status variables as they are useful for finding missing pings [outdat["data_vars"].pop(var) for var in status0_vars] # Set coordinate system @@ -668,8 +703,25 @@ def _reorg(dat): return outdat +def _clean_dp_skips(data): + """ + Removes zeros from interwoven measurements taken in a dual profile + configuration. + """ + + for id in data: + if id == "filehead_config": + continue + # Check where 'ver' is zero (should be 1 (for bt) or 3 (everything else)) + skips = np.where(data[id]["ver"] != 0) + for var in data[id]: + if var not in ["units", "long_name", "standard_name"]: + data[id][var] = np.squeeze(data[id][var][..., skips], axis=-2) + + def _reduce(data): - """This function takes the output from `reorg`, and further simplifies the + """ + This function takes the output from `reorg`, and further simplifies the data. Mostly this is combining system, environmental, and orientation data --- from different data structures within the same ensemble --- by averaging. @@ -697,13 +749,15 @@ def _reduce(data): dc["range_avg"] = (np.arange(dv["vel_avg"].shape[1]) + 1) * da[ "cell_size_avg" ] + da["blank_dist_avg"] - dv["orientmat"] = dv.pop("orientmat_avg") + if "orientmat" not in dv: + dv["orientmat"] = dv.pop("orientmat_avg") tmat = da["filehead_config"]["XFAVG"] da["fs"] = da["filehead_config"]["PLAN"]["MIAVG"] da["avg_interval_sec"] = da["filehead_config"]["AVG"]["AI"] da["bandwidth"] = da["filehead_config"]["AVG"]["BW"] if "vel_b5" in dv: - dc["range_b5"] = (np.arange(dv["vel_b5"].shape[1]) + 1) * da[ + # vel_b5 is sometimes shape 2 and sometimes shape 3 + dc["range_b5"] = (np.arange(dv["vel_b5"].shape[-2]) + 1) * da[ "cell_size_b5" ] + da["blank_dist_b5"] if "echo_echo" in dv: @@ -736,3 +790,65 @@ def _reduce(data): if "time" in val: time = val dc["time"] = dc[time] + + +def split_dp_datasets(ds): + """ + Splits a dataset containing dual profiles into individual profiles + """ + + # Figure out which variables belong to which profile based on length of time variables + t_dict = {} + for t in ds.coords: + if "time" in t: + t_dict[t] = ds[t].size + + other_coords = [] + for key, val in t_dict.items(): + if val != t_dict["time"]: + if key.endswith("altraw"): + # altraw goes with burst, altraw_avg goes with avg + continue + other_coords.append(key) + # Fetch variables, coordinates, and attrs for second profiling configuration + other_vars = [ + v for v in ds.data_vars if any(x in ds[v].coords for x in other_coords) + ] + other_tags = [s.split("_")[-1] for s in other_coords] + other_coords += [v for v in ds.coords if any(x in v for x in other_tags)] + other_attrs = [s for s in ds.attrs if any(x in s for x in other_tags)] + critical_attrs = [ + "inst_model", + "inst_make", + "inst_type", + "fs", + "orientation", + "orient_status", + "has_imu", + "beam_angle", + ] + + # Create second dataset + ds2 = type(ds)() + for a in other_attrs + critical_attrs: + ds2.attrs[a] = ds.attrs[a] + for v in other_vars: + ds2[v] = ds[v] + # Set rotate_vars + rotate_vars2 = [v for v in ds.attrs["rotate_vars"] if v in other_vars] + ds2.attrs["rotate_vars"] = rotate_vars2 + # Set orientation matricies + ds2["beam2inst_orientmat"] = ds["beam2inst_orientmat"] + ds2 = ds2.rename({"orientmat_" + other_tags[0]: "orientmat"}) + # Set original coordinate system + cy = ds2.attrs["coord_sys_axes_" + other_tags[0]] + ds2.attrs["coord_sys"] = {"XYZ": "inst", "ENU": "earth", "beam": "beam"}[cy] + ds2 = _set_coords(ds2, ref_frame=ds2.coord_sys) + + # Clean up first dataset + [ds.attrs.pop(ky) for ky in other_attrs] + ds = ds.drop_vars(other_vars + other_coords) + for itm in rotate_vars2: + ds.attrs["rotate_vars"].remove(itm) + + return ds, ds2 diff --git a/mhkit/dolfyn/io/nortek2_lib.py b/mhkit/dolfyn/io/nortek2_lib.py index 047194fac..3336651f5 100644 --- a/mhkit/dolfyn/io/nortek2_lib.py +++ b/mhkit/dolfyn/io/nortek2_lib.py @@ -157,9 +157,11 @@ def _create_index(infile, outfile, init_pos, eof, debug): fin.seek(seek_2ens[dat[2]], 1) ens[idk] = struct.unpack(" 0: - # Covers all id keys saved in "burst mode" - ens[idk] = last_ens[idk] + 1 + if last_ens[idk] > 0: + if (ens[idk] == 1) or (ens[idk] < last_ens[idk]): + # Covers all id keys saved in "burst mode" + # Covers ID keys not saved in sequential order + ens[idk] = last_ens[idk] + 1 if last_ens[idk] > 0 and last_ens[idk] != ens[idk]: N[idk] += 1 @@ -188,6 +190,7 @@ def _create_index(infile, outfile, init_pos, eof, debug): last_ens[idk] = ens[idk] if debug: + # File Position: Valid ID keys (1A, 10), Hex ID, Length in bytes, Ensemble #, Last Ensemble Found' # hex: [18, 15, 1C, 17] = [vel_b5, vel, echo, bt] logging.info( "%10d: %02X, %d, %02X, %d, %d, %d, %d\n" @@ -213,7 +216,7 @@ def _create_index(infile, outfile, init_pos, eof, debug): print(" Done.") -def _check_index(idx, infile, fix_hw_ens=False): +def _check_index(idx, infile, fix_hw_ens=False, dp=False): uid = np.unique(idx["ID"]) if fix_hw_ens: hwe = idx["hw_ens"] @@ -223,13 +226,29 @@ def _check_index(idx, infile, fix_hw_ens=False): ens = idx["ens"] N_id = len(uid) FLAG = False + + # Are there better ways to detect dual profile? + if (21 in uid) and (22 in uid): + warnings.warn("Dual Profile detected... Two datasets will be returned.") + dp = True + # This loop fixes 'skips' inside the file for id in uid: # These are the indices for this ID inds = np.nonzero(idx["ID"] == id)[0] - # These are bad steps in the indices for this ID ibad = np.nonzero(np.diff(inds) > N_id)[0] + # Check if spacing is equal for dual profiling ADCPs + if dp: + skip_size = np.diff(ibad) + n_skip, count = np.unique(skip_size, return_counts=True) + # If multiple skips are of the same size, assume okay + for n, c in zip(n_skip, count): + if c > 1: + skip_size[skip_size == n] = 0 + # assume last "ibad" element is always good for dp's + mask = np.append(skip_size, 0).astype(bool) if any(skip_size) else [] + ibad = ibad[mask] for ib in ibad: FLAG = True # The ping number reported here may not be quite right if @@ -242,16 +261,7 @@ def _check_index(idx, infile, fix_hw_ens=False): hwe[inds[(ib + 1) :]] += 1 ens[inds[(ib + 1) :]] += 1 - # This block fixes skips that originate from before this file. - delta = max(hwe[:N_id]) - hwe[:N_id] - for d, id in zip(delta, idx["ID"][:N_id]): - if d != 0: - FLAG = True - hwe[id == idx["ID"]] += d - ens[id == idx["ID"]] += d - - if np.any(np.diff(ens) > 1) and FLAG: - idx["ens"] = np.unwrap(hwe.astype(np.int64), period=period) - hwe[0] + return dp def _boolarray_firstensemble_ping(index): @@ -264,7 +274,7 @@ def _boolarray_firstensemble_ping(index): return dens -def get_index(infile, pos=0, eof=2**32, reload=False, debug=False): +def get_index(infile, pos=0, eof=2**32, rebuild=False, debug=False, dp=False): """ This function reads ad2cp.index files @@ -284,7 +294,7 @@ def get_index(infile, pos=0, eof=2**32, reload=False, debug=False): """ index_file = infile + ".index" - if not path.isfile(index_file) or reload: + if not path.isfile(index_file) or rebuild or debug: _create_index(infile, index_file, pos, eof, debug) f = open(_abspath(index_file), "rb") file_head = f.read(12) @@ -296,8 +306,8 @@ def get_index(infile, pos=0, eof=2**32, reload=False, debug=False): f.seek(0, 0) out = np.fromfile(f, dtype=_index_dtype[index_ver]) f.close() - _check_index(out, infile) - return out + dp = _check_index(out, infile, dp=dp) + return out, dp def crop_ensembles(infile, outfile, range): @@ -319,7 +329,7 @@ def crop_ensembles(infile, outfile, range): 2 element list of start and end ensemble (or time index) """ - idx = get_index(infile) + idx, dp = get_index(infile) with open(_abspath(infile), "rb") as fin: with open(_abspath(outfile), "wb") as fout: fout.write(fin.read(idx["pos"][0])) @@ -485,7 +495,8 @@ def _beams_cy_int2dict(val, id): """Convert the beams/coordinate-system bytes to a dict of values.""" if id == 28: # 0x1C (echosounder) return dict(n_cells=val) - + elif id in [26, 31]: + return dict(n_cells=val & (2**10 - 1), cy="beam", n_beams=1) return dict( n_cells=val & (2**10 - 1), cy=["ENU", "XYZ", "beam", None][val >> 10 & 3], @@ -544,22 +555,32 @@ def _calc_config(index): ids = np.unique(index["ID"]) config = {} for id in ids: - if id not in [21, 22, 23, 24, 26, 28]: + if id not in [21, 22, 23, 24, 26, 28, 31]: continue if id == 23: type = "bt" - elif id == 22: + elif (id == 22) or (id == 31): type = "avg" else: type = "burst" inds = index["ID"] == id _config = index["config"][inds] _beams_cy = index["beams_cy"][inds] + # Check that these variables are consistent if not _isuniform(_config): raise Exception("config are not identical for id: 0x{:X}.".format(id)) if not _isuniform(_beams_cy): - raise Exception("beams_cy are not identical for id: 0x{:X}.".format(id)) + err = True + if id == 23: + # change in "n_cells" doesn't matter + lob = np.unique(_beams_cy) + beams = list(map(_beams_cy_int2dict, lob, 23 * np.ones(lob.size))) + if all([d["cy"] for d in beams]) and all([d["n_beams"] for d in beams]): + err = False + if err: + raise Exception("beams_cy are not identical for id: 0x{:X}.".format(id)) + # Now that we've confirmed they are the same: config[id] = _headconfig_int2dict(_config[0], mode=type) config[id].update(_beams_cy_int2dict(_beams_cy[0], id)) diff --git a/mhkit/tests/dolfyn/test_read_adp.py b/mhkit/tests/dolfyn/test_read_adp.py index 30ad34d98..b4eceae01 100644 --- a/mhkit/tests/dolfyn/test_read_adp.py +++ b/mhkit/tests/dolfyn/test_read_adp.py @@ -34,6 +34,7 @@ dat_sig_skip = load("Sig_SkippedPings01.nc") dat_sig_badt = load("Sig1000_BadTime01.nc") dat_sig5_leiw = load("Sig500_last_ensemble_is_whole.nc") +dat_sig_dp2 = load("dual_profile.nc") class io_adp_testcase(unittest.TestCase): @@ -91,12 +92,16 @@ def test_io_nortek(self): def test_io_nortek2(self): nens = 100 - td_sig = read("BenchFile01.ad2cp", nens=nens) - td_sig_i = read("Sig1000_IMU.ad2cp", userdata=False, nens=nens) - td_sig_i_ud = read("Sig1000_IMU.ad2cp", nens=nens) - td_sig_ieb = read("VelEchoBT01.ad2cp", nens=nens) - td_sig_ie = read("Sig500_Echo.ad2cp", nens=nens) - td_sig_tide = read("Sig1000_tidal.ad2cp", nens=nens) + td_sig = read("BenchFile01.ad2cp", nens=nens, rebuild_index=True) + td_sig_i = read( + "Sig1000_IMU.ad2cp", userdata=False, nens=nens, rebuild_index=True + ) + td_sig_i_ud = read("Sig1000_IMU.ad2cp", nens=nens, rebuild_index=True) + td_sig_ieb = read("VelEchoBT01.ad2cp", nens=nens, rebuild_index=True) + td_sig_ie = read("Sig500_Echo.ad2cp", nens=nens, rebuild_index=True) + td_sig_tide = read("Sig1000_tidal.ad2cp", nens=nens, rebuild_index=True) + # Only need to test 2nd dataset + td_sig_dp1, td_sig_dp2 = read("dual_profile.ad2cp") with pytest.warns(UserWarning): # This issues a warning... @@ -117,6 +122,7 @@ def test_io_nortek2(self): os.remove(tb.exdt("Sig_SkippedPings01.ad2cp.index")) os.remove(tb.exdt("Sig500_last_ensemble_is_whole.ad2cp.index")) os.remove(tb.rfnm("Sig1000_BadTime01.ad2cp.index")) + os.remove(tb.exdt("dual_profile.ad2cp.index")) if make_data: save(td_sig, "BenchFile01.nc") @@ -128,6 +134,7 @@ def test_io_nortek2(self): save(td_sig_skip, "Sig_SkippedPings01.nc") save(td_sig_badt, "Sig1000_BadTime01.nc") save(td_sig5_leiw, "Sig500_last_ensemble_is_whole.nc") + save(td_sig_dp2, "dual_profile.nc") return assert_allclose(td_sig, dat_sig, atol=1e-6) @@ -139,6 +146,7 @@ def test_io_nortek2(self): assert_allclose(td_sig5_leiw, dat_sig5_leiw, atol=1e-6) assert_allclose(td_sig_skip, dat_sig_skip, atol=1e-6) assert_allclose(td_sig_badt, dat_sig_badt, atol=1e-6) + assert_allclose(td_sig_dp2, dat_sig_dp2, atol=1e-6) def test_nortek2_crop(self): # Test file cropping function