Implementation of ASN regions of interest (ROIs)

xdas/io/asn.py

@@ -1,4 +1,5 @@
 import json
+from bisect import bisect_left
 
 import h5py
 import numpy as np
@@ -10,16 +11,45 @@
 from ..virtual import VirtualSource
 
 
-def read(fname):
+def read(fname: str) -> DataArray:
     with h5py.File(fname, "r") as file:
         header = file["header"]
+        demod = file["demodSpec"]
+
         t0 = np.datetime64(round(header["time"][()] * 1e9), "ns")
         dt = np.timedelta64(round(1e9 * header["dt"][()]), "ns")
-        dx = header["dx"][()] * np.median(np.diff(header["channels"]))
+        dx = float(header["dx"][()])  # Note: dx before (internal) downsampling!
         data = VirtualSource(file["data"])
-    nt, nx = data.shape
+
+        # Get the optical distance for all the recorded channels (after downsampling)
+        # Note that this vector is not continuous for more than one ROI
+        all_dists = file["cableSpec"]["sensorDistances"][...]
+
+        # Buffer for the data index at which each ROI starts/stops
+        dist_tie_inds = []
+        # Buffer for the optical distance at which each ROI starts/stops
+        dist_tie_vals = []
+
+        # Loop over ROIs, get the start/stop index before downsampling
+        for n_start, n_end in zip(demod["roiStart"], demod["roiEnd"]):
+            # Get the index where the ROI starts based on the position in the
+            # distance vector. This solves the issue of rounding during decimation
+            i = bisect_left(all_dists, n_start * dx)
+            # Append the data index and optical distance to the buffers
+            dist_tie_inds.append(i)
+            dist_tie_vals.append(float(all_dists[i]))
+
+            # Repeat the procedure for the index/distance at which the ROI ends.
+            # A "discontinuity" in the interpolation scheme is created in the
+            # following way: n_roi = [start, stop-1, stop, start, stop-1, stop, ...]
+            i = bisect_left(all_dists, n_end * dx)
+            for j in reversed(range(2)):
+                dist_tie_inds.append(i-j)
+                dist_tie_vals.append(float(all_dists[i-j]))
+
+    nt = data.shape[0]
     time = {"tie_indices": [0, nt - 1], "tie_values": [t0, t0 + (nt - 1) * dt]}
-    distance = {"tie_indices": [0, nx - 1], "tie_values": [0.0, (nx - 1) * dx]}
+    distance = {"tie_indices": dist_tie_inds, "tie_values": dist_tie_vals}
     return DataArray(data, {"time": time, "distance": distance})