canvod.readers API Reference¶

class GNSSDataReader(BaseModel, ABC):
    """Abstract base class for all GNSS data format readers.

    All readers must:
    1. Inherit from this class
    2. Implement all abstract methods
    3. Return xarray.Dataset that passes :func:`validate_dataset`
    4. Provide file hash for deduplication

    This ensures compatibility with:
    - canvod-vod: VOD calculation
    - canvod-store: MyIcechunkStore storage
    - canvod-grids: Grid projection operations

    Subclasses may override ``model_config`` to set ``frozen``, ``extra``,
    etc.  The base class provides ``arbitrary_types_allowed=True`` which is
    needed by readers that use ``pint.Quantity`` or similar third-party types.

    Examples
    --------
    >>> class Rnxv3Obs(GNSSDataReader):
    ...     def to_ds(self, **kwargs) -> xr.Dataset:
    ...         # Implementation
    ...         return dataset
    ...
    >>> reader = Rnxv3Obs(fpath="station.24o")
    >>> ds = reader.to_ds()
    >>> validate_dataset(ds)
    """

    model_config = ConfigDict(arbitrary_types_allowed=True)

    fpath: Path

    @field_validator("fpath")
    @classmethod
    def _validate_fpath(cls, v: Path) -> Path:
        """Validate that the file path points to an existing file."""
        v = Path(v)
        if not v.is_file():
            raise FileNotFoundError(f"File not found: {v}")
        return v

    @property
    def source_format(self) -> str:
        """Return the format identifier for this reader (e.g. ``"rinex3"``, ``"sbf"``)."""
        return "rinex3"

    @property
    @abstractmethod
    def file_hash(self) -> str:
        """Return SHA256 hash of file for deduplication.

        Used by MyIcechunkStore to avoid duplicate ingestion.
        Must be deterministic and reproducible.

        Returns
        -------
        str
            Short hash (16 chars) or full hash of file content
        """

    @abstractmethod
    def to_ds(
        self,
        keep_data_vars: list[str] | None = None,
        **kwargs: object,
    ) -> xr.Dataset:
        """Convert data to xarray.Dataset.

        Must return Dataset with structure:
        - Dims: (epoch, sid)
        - Coords: epoch, sid, sv, system, band, code, freq_*
        - Data vars: At minimum SNR
        - Attrs: Must include "File Hash"

        Parameters
        ----------
        keep_data_vars : list of str, optional
            Data variables to include. If None, includes all available.
        **kwargs
            Implementation-specific parameters

        Returns
        -------
        xr.Dataset
            Dataset that passes :func:`validate_dataset`.
        """

    @abstractmethod
    def iter_epochs(self) -> Iterator[object]:
        """Iterate over epochs in the file.

        Yields
        ------
        Epoch
            Parsed epoch with satellites and observations.
        """

    def to_ds_and_auxiliary(
        self,
        keep_data_vars: list[str] | None = None,
        **kwargs: object,
    ) -> tuple[xr.Dataset, dict[str, xr.Dataset]]:
        """Produce the obs dataset and any auxiliary datasets in a single call.

        Default: calls ``to_ds(**kwargs)`` and returns an empty auxiliary dict.
        Readers that produce metadata (e.g. SBF) override this to collect both
        in a single file scan.

        Returns
        -------
        tuple[xr.Dataset, dict[str, xr.Dataset]]
            ``(obs_ds, {"name": aux_ds, ...})``.  Auxiliary dict is empty for
            readers with no extra data (RINEX v2/v3).
        """
        return self.to_ds(keep_data_vars=keep_data_vars, **kwargs), {}

    def _build_attrs(self) -> dict[str, str]:
        """Build standard global attributes for the output Dataset.

        Reads institution/author from config, adds timestamp, version,
        and the file hash.

        Returns
        -------
        dict[str, str]
            Ready-to-use attrs dict.
        """
        from canvod.readers.gnss_specs.metadata import get_global_attrs
        from canvod.readers.gnss_specs.utils import get_version_from_pyproject

        attrs = get_global_attrs()
        attrs["Created"] = datetime.now(UTC).isoformat()
        attrs["Software"] = (
            f"{attrs['Software']}, Version: {get_version_from_pyproject()}"
        )
        attrs["File Hash"] = self.file_hash
        return attrs

    @property
    @abstractmethod
    def start_time(self) -> datetime:
        """Return start time of observations.

        Returns
        -------
        datetime
            First observation timestamp in the file.
        """

    @property
    @abstractmethod
    def end_time(self) -> datetime:
        """Return end time of observations.

        Returns
        -------
        datetime
            Last observation timestamp in the file.
        """

    @property
    @abstractmethod
    def systems(self) -> list[str]:
        """Return list of GNSS systems in file.

        Returns
        -------
        list of str
            System identifiers: 'G', 'R', 'E', 'C', 'J', 'S', 'I'
        """

    @property
    def num_epochs(self) -> int:
        """Return number of epochs in file.

        Default implementation iterates epochs.  Subclasses may override
        with a faster approach.

        Returns
        -------
        int
            Total number of observation epochs.
        """
        return sum(1 for _ in self.iter_epochs())

    @property
    @abstractmethod
    def num_satellites(self) -> int:
        """Return total number of unique satellites observed.

        Returns
        -------
        int
            Count of unique satellite vehicles across all systems.
        """

    def __repr__(self) -> str:
        """Return the string representation."""
        return f"{self.__class__.__name__}(file='{self.fpath.name}')"

class SignalID(BaseModel):
    """Validated signal identifier (SV + band + code).

    >>> sid = SignalID(sv="G01", band="L1", code="C")
    >>> str(sid)
    'G01|L1|C'
    >>> sid.system
    'G'
    """

    model_config = ConfigDict(frozen=True)

    sv: str
    band: str
    code: str

    @field_validator("sv")
    @classmethod
    def _validate_sv(cls, v: str) -> str:
        if not SV_PATTERN.match(v):
            raise ValueError(
                f"Invalid SV: {v!r} — expected system letter + 2-digit PRN "
                f"(e.g. 'G01'). Valid systems: G, R, E, C, J, S, I"
            )
        return v

    @property
    def system(self) -> str:
        """GNSS system letter (e.g. 'G' for GPS)."""
        return self.sv[0]

    @property
    def sid(self) -> str:
        """Full signal ID string ('SV|band|code')."""
        return f"{self.sv}|{self.band}|{self.code}"

    def __str__(self) -> str:
        return self.sid

    def __hash__(self) -> int:
        return hash(self.sid)

    def __eq__(self, other: object) -> bool:
        if isinstance(other, SignalID):
            return self.sid == other.sid
        return NotImplemented

    @classmethod
    def from_string(cls, sid_str: str) -> SignalID:
        """Parse a signal ID string ('SV|band|code') into a SignalID.

        Parameters
        ----------
        sid_str : str
            Signal ID in 'SV|band|code' format (e.g. 'G01|L1|C').

        Returns
        -------
        SignalID
            Validated signal identifier.

        Raises
        ------
        ValueError
            If the string does not have exactly three pipe-separated parts.
        """
        parts = sid_str.split("|")
        if len(parts) != 3:
            raise ValueError(
                f"Invalid SID format: {sid_str!r} — expected 'SV|band|code'"
            )
        return cls(sv=parts[0], band=parts[1], code=parts[2])

class DatasetBuilder:
    """Guided builder for constructing valid GNSSDataReader output Datasets.

    Handles coordinate arrays, dtype enforcement, frequency resolution,
    and contract validation automatically.

    Parameters
    ----------
    reader : GNSSDataReader
        The reader instance (used for ``_build_attrs()`` and file hash).
    aggregate_glonass_fdma : bool, optional
        Whether to aggregate GLONASS FDMA channels (default True).

    Examples
    --------
    >>> builder = DatasetBuilder(reader)
    >>> for epoch in reader.iter_epochs():
    ...     ei = builder.add_epoch(epoch.timestamp)
    ...     for obs in epoch.observations:
    ...         sig = builder.add_signal(sv="G01", band="L1", code="C")
    ...         builder.set_value(ei, sig, "SNR", 42.0)
    >>> ds = builder.build()   # validated Dataset
    """

    def __init__(
        self,
        reader: GNSSDataReader,
        *,
        aggregate_glonass_fdma: bool = True,
    ) -> None:
        self._reader = reader
        self._mapper = SignalIDMapper(aggregate_glonass_fdma=aggregate_glonass_fdma)
        self._signals: dict[str, SignalID] = {}
        self._epochs: list[datetime] = []
        self._values: dict[str, dict[tuple[int, str], float]] = {}

    def add_epoch(self, timestamp: datetime) -> int:
        """Register an epoch timestamp. Returns epoch index."""
        self._epochs.append(timestamp)
        return len(self._epochs) - 1

    def add_signal(self, sv: str, band: str, code: str) -> SignalID:
        """Register a signal (idempotent). Returns validated SignalID."""
        sig = SignalID(sv=sv, band=band, code=code)
        self._signals[sig.sid] = sig
        return sig

    def set_value(
        self,
        epoch_idx: int,
        signal: SignalID | str,
        var: str,
        value: float,
    ) -> None:
        """Set a data value for a given epoch, signal, and variable.

        Parameters
        ----------
        epoch_idx : int
            Index returned by :meth:`add_epoch`.
        signal : SignalID or str
            Signal identifier (SignalID or 'SV|band|code' string).
        var : str
            Variable name (e.g. 'SNR', 'Pseudorange', 'Phase').
        value : float
            The observation value.
        """
        sid = str(signal)
        if var not in self._values:
            self._values[var] = {}
        self._values[var][(epoch_idx, sid)] = value

    def build(
        self,
        keep_data_vars: list[str] | None = None,
        extra_attrs: dict[str, str] | None = None,
    ) -> xr.Dataset:
        """Build, validate, and return the Dataset.

        1. Sorts signals alphabetically
        2. Resolves frequencies from band names via SignalIDMapper
        3. Constructs coordinate arrays with correct dtypes (float32 for freq)
        4. Attaches CF-compliant metadata from COORDS_METADATA
        5. Calls validate_dataset() before returning

        Parameters
        ----------
        keep_data_vars : list of str, optional
            If provided, only include these data variables.  If ``None``,
            includes all variables that had values set.
        extra_attrs : dict, optional
            Additional global attributes to merge into the Dataset.

        Returns
        -------
        xr.Dataset
            Validated Dataset with dimensions ``(epoch, sid)``.
        """
        sorted_sids = sorted(self._signals)
        sid_to_idx = {sid: i for i, sid in enumerate(sorted_sids)}
        n_epochs = len(self._epochs)
        n_sids = len(sorted_sids)

        # --- Coordinate arrays ---
        epoch_arr = [
            np.datetime64(ts.replace(tzinfo=None) if ts.tzinfo else ts, "ns")
            for ts in self._epochs
        ]
        sv_arr = np.array([self._signals[s].sv for s in sorted_sids], dtype=object)
        system_arr = np.array(
            [self._signals[s].system for s in sorted_sids], dtype=object
        )
        band_arr = np.array([self._signals[s].band for s in sorted_sids], dtype=object)
        code_arr = np.array([self._signals[s].code for s in sorted_sids], dtype=object)

        # Frequency resolution via SignalIDMapper
        freq_center = np.array(
            [
                self._mapper.get_band_frequency(self._signals[s].band) or np.nan
                for s in sorted_sids
            ],
            dtype=np.float32,
        )
        bandwidths = np.array(
            [
                self._mapper.get_band_bandwidth(self._signals[s].band) or 0.0
                for s in sorted_sids
            ],
            dtype=np.float32,
        )
        freq_min = (freq_center - bandwidths / 2).astype(np.float32)
        freq_max = (freq_center + bandwidths / 2).astype(np.float32)

        # --- Determine which variables to include ---
        all_vars = set(self._values.keys())
        if keep_data_vars is not None:
            vars_to_build = [v for v in keep_data_vars if v in all_vars]
        else:
            vars_to_build = sorted(all_vars)

        # --- Data variable arrays ---
        data_vars: dict[str, tuple] = {}
        for var in vars_to_build:
            dtype = DTYPES.get(var, np.dtype("float32"))
            fill = np.nan if np.issubdtype(dtype, np.floating) else -1
            arr = np.full((n_epochs, n_sids), fill, dtype=dtype)

            for (ei, sid_str), val in self._values[var].items():
                if sid_str in sid_to_idx:
                    arr[ei, sid_to_idx[sid_str]] = val

            meta = _VAR_METADATA.get(var, {})
            data_vars[var] = (("epoch", "sid"), arr, meta)

        # --- Coordinates ---
        coords = {
            "epoch": ("epoch", epoch_arr, COORDS_METADATA["epoch"]),
            "sid": xr.DataArray(
                np.array(sorted_sids, dtype=object),
                dims=["sid"],
                attrs=COORDS_METADATA["sid"],
            ),
            "sv": ("sid", sv_arr, COORDS_METADATA["sv"]),
            "system": ("sid", system_arr, COORDS_METADATA["system"]),
            "band": ("sid", band_arr, COORDS_METADATA["band"]),
            "code": ("sid", code_arr, COORDS_METADATA["code"]),
            "freq_center": ("sid", freq_center, COORDS_METADATA["freq_center"]),
            "freq_min": ("sid", freq_min, COORDS_METADATA["freq_min"]),
            "freq_max": ("sid", freq_max, COORDS_METADATA["freq_max"]),
        }

        # --- Global attributes ---
        attrs = self._reader._build_attrs()
        if extra_attrs:
            attrs.update(extra_attrs)

        ds = xr.Dataset(data_vars=data_vars, coords=coords, attrs=attrs)

        # Validate before returning
        validate_dataset(ds, required_vars=keep_data_vars)

        return ds

class DatasetStructureValidator(BaseModel):
    """Validates that an xarray.Dataset meets the GNSSDataReader contract.

    Wraps a Dataset and checks it against the contract constants above.
    Use this in tests and reader implementations to catch structural errors
    early with clear messages.

    Examples
    --------
    >>> validator = DatasetStructureValidator(dataset=ds)
    >>> validator.validate_all()          # raises ValueError on any violation
    >>> validator.validate_dimensions()   # check just one aspect
    """

    model_config = ConfigDict(arbitrary_types_allowed=True)

    dataset: xr.Dataset

    def validate_all(self, required_vars: list[str] | None = None) -> None:
        """Run all validations, collecting **all** errors.

        Delegates to :func:`validate_dataset` so the logic lives in one place.
        """
        validate_dataset(self.dataset, required_vars=required_vars)

    def validate_dimensions(self) -> None:
        """Check that required dimensions (epoch, sid) exist."""
        missing = set(REQUIRED_DIMS) - set(self.dataset.dims)
        if missing:
            raise ValueError(f"Missing required dimensions: {missing}")

    def validate_coordinates(self) -> None:
        """Check that required coordinates exist with correct dtypes."""
        for coord, expected_dtype in REQUIRED_COORDS.items():
            if coord not in self.dataset.coords:
                raise ValueError(f"Missing required coordinate: {coord}")
            actual = str(self.dataset[coord].dtype)
            if expected_dtype == "object":
                is_valid_string = actual == "object" or actual.startswith("StringDType")
                if not is_valid_string:
                    raise ValueError(
                        f"Coordinate {coord}: expected string (object/StringDType), got {actual}"
                    )
            elif expected_dtype not in actual:
                raise ValueError(
                    f"Coordinate {coord}: expected {expected_dtype}, got {actual}"
                )

    def validate_data_variables(self, required_vars: list[str] | None = None) -> None:
        """Check that required data variables exist with correct dims."""
        if required_vars is None:
            required_vars = list(DEFAULT_REQUIRED_VARS)
        missing = set(required_vars) - set(self.dataset.data_vars)
        if missing:
            raise ValueError(f"Missing required data variables: {missing}")
        for var in self.dataset.data_vars:
            if self.dataset[var].dims != REQUIRED_DIMS:
                raise ValueError(
                    f"Variable {var}: expected dims {REQUIRED_DIMS}, "
                    f"got {self.dataset[var].dims}"
                )

    def validate_attributes(self) -> None:
        """Check that required global attributes are present."""
        missing = REQUIRED_ATTRS - set(self.dataset.attrs.keys())
        if missing:
            raise ValueError(f"Missing required attributes: {missing}")

class Rnxv3Obs(GNSSDataReader):
    """RINEX v3.04 observation reader.

    Attributes
    ----------
    fpath : Path
        Path to the RINEX observation file.
    polarization : str, default "RHCP"
        Polarization label for observables.
    completeness_mode : {"strict", "warn", "off"}, default "strict"
        Behavior when epoch completeness checks fail.
    expected_dump_interval : str or pint.Quantity, optional
        Expected file dump interval for completeness validation.
    expected_sampling_interval : str or pint.Quantity, optional
        Expected sampling interval for completeness validation.
    apply_overlap_filter : bool, default False
        Whether to filter overlapping signal groups.
    overlap_preferences : dict[str, str], optional
        Preferred signals for overlap resolution.
    aggregate_glonass_fdma : bool, optional
        Whether to aggregate GLONASS FDMA channels.

    Notes
    -----
    Inherits ``fpath``, its validator, and ``arbitrary_types_allowed``
    from :class:`GNSSDataReader`.

    """

    model_config = ConfigDict(frozen=True)

    polarization: str = "RHCP"

    completeness_mode: Literal["strict", "warn", "off"] = "strict"
    expected_dump_interval: str | pint.Quantity | None = None
    expected_sampling_interval: str | pint.Quantity | None = None

    apply_overlap_filter: bool = False
    overlap_preferences: dict[str, str] | None = None

    aggregate_glonass_fdma: bool = True

    _header: Rnxv3Header = PrivateAttr()
    _signal_mapper: SignalIDMapper = PrivateAttr()

    _lines: list[str] = PrivateAttr()
    _file_hash: str = PrivateAttr()
    _cached_epoch_batches: list[tuple[int, int]] | None = PrivateAttr(default=None)

    @model_validator(mode="after")
    def _post_init(self) -> Self:
        """Initialize derived state after validation."""
        # Load header once
        self._header = Rnxv3Header.from_file(self.fpath)

        # Initialize signal mapper
        self._signal_mapper = SignalIDMapper(
            aggregate_glonass_fdma=self.aggregate_glonass_fdma
        )

        # Optionally auto-check completeness
        if self.completeness_mode != "off":
            try:
                self.validate_epoch_completeness(
                    dump_interval=self.expected_dump_interval,
                    sampling_interval=self.expected_sampling_interval,
                )
            except MissingEpochError as e:
                if self.completeness_mode == "strict":
                    raise
                warnings.warn(str(e), RuntimeWarning, stacklevel=2)

        # Cache file lines
        self._lines = self._load_file()

        return self

    @property
    def header(self) -> Rnxv3Header:
        """Expose validated header (read-only).

        Returns
        -------
        Rnxv3Header
            Parsed and validated RINEX header.

        """
        return self._header

    def __str__(self) -> str:
        """Return a human-readable summary."""
        return (
            f"{self.__class__.__name__}:\n"
            f"  File Path: {self.fpath}\n"
            f"  Header: {self.header}\n"
            f"  Polarization: {self.polarization}\n"
        )

    def __repr__(self) -> str:
        """Return a concise representation for debugging."""
        return f"{self.__class__.__name__}(fpath={self.fpath})"

    def _load_file(self) -> list[str]:
        """Read file once, cache lines, and compute hash.

        Returns
        -------
        list[str]
            File contents split into lines.

        """
        if not hasattr(self, "_lines"):
            h = hashlib.sha256()
            with self.fpath.open("rb") as f:  # binary mode for consistent hash
                data = f.read()
                h.update(data)
                self._lines = data.decode("utf-8", errors="replace").splitlines()
            self._file_hash = h.hexdigest()[:16]  # short hash for storage
        return self._lines

    @property
    def file_hash(self) -> str:
        """Return cached SHA256 short hash of the file content.

        Returns
        -------
        str
            16-character short hash for deduplication.

        """
        return self._file_hash

    @property
    def start_time(self) -> datetime:
        """Return start time of observations from header.

        Returns
        -------
        datetime
            First observation timestamp.

        """
        return min(self.header.t0.values())

    @property
    def end_time(self) -> datetime:
        """Return end time of observations from last epoch.

        Returns
        -------
        datetime
            Last observation timestamp.

        """
        last_epoch = None
        for epoch in self.iter_epochs():
            last_epoch = epoch
        if last_epoch:
            return self.get_datetime_from_epoch_record_info(last_epoch.info)
        return self.start_time

    @property
    def systems(self) -> list[str]:
        """Return list of GNSS systems in file.

        Returns
        -------
        list of str
            System identifiers (G, R, E, C, J, S, I).

        """
        if self.header.systems == "M":
            return list(self.header.obs_codes_per_system.keys())
        return [self.header.systems]

    @property
    def num_epochs(self) -> int:
        """Return number of epochs in file.

        Returns
        -------
        int
            Total epoch count.

        """
        return len(list(self.get_epoch_record_batches()))

    @property
    def num_satellites(self) -> int:
        """Return total number of unique satellites observed.

        Returns
        -------
        int
            Count of unique satellite vehicles across all systems.

        """
        satellites = set()
        for epoch in self.iter_epochs():
            for sat in epoch.data:
                satellites.add(sat.sv)
        return len(satellites)

    def get_epoch_record_batches(
        self, epoch_record_indicator: str = EPOCH_RECORD_INDICATOR
    ) -> list[tuple[int, int]]:
        """Get the start and end line numbers for each epoch in the file.

        Parameters
        ----------
        epoch_record_indicator : str, default '>'
            Character marking epoch record lines.

        Returns
        -------
        list of tuple of int
            List of (start_line, end_line) pairs for each epoch.

        """
        if self._cached_epoch_batches is not None:
            return self._cached_epoch_batches

        lines = self._load_file()
        starts = [
            i for i, line in enumerate(lines) if line.startswith(epoch_record_indicator)
        ]
        starts.append(len(lines))  # Add EOF
        self._cached_epoch_batches = [
            (start, starts[i + 1])
            for i, start in enumerate(starts)
            if i + 1 < len(starts)
        ]
        return self._cached_epoch_batches

    def parse_observation_slice(
        self,
        slice_text: str,
    ) -> tuple[float | None, int | None, int | None]:
        """Parse a RINEX observation slice into value, LLI, and SSI.

        Enhanced to handle both standard 16-character format and
        variable-length records.

        Parameters
        ----------
        slice_text : str
            Observation slice to parse.

        Returns
        -------
        tuple[float | None, int | None, int | None]
            Parsed (value, LLI, SSI) tuple.

        """
        if not slice_text or not slice_text.strip():
            return None, None, None

        try:
            # Method 1: Standard RINEX format with decimal at position -6
            if (
                len(slice_text) >= OBS_SLICE_MIN_LEN
                and len(slice_text) <= OBS_SLICE_MAX_LEN
                and slice_text[OBS_SLICE_DECIMAL_POS] == "."
            ):
                slice_chars = list(slice_text)
                ssi = slice_chars.pop(-1) if len(slice_chars) > 0 else ""
                lli = slice_chars.pop(-1) if len(slice_chars) > 0 else ""

                # Convert LLI and SSI
                lli = int(lli) if lli.strip() and lli.isdigit() else None
                ssi = int(ssi) if ssi.strip() and ssi.isdigit() else None

                # Convert value
                value_str = "".join(slice_chars).strip()
                if value_str:
                    value = float(value_str)
                    return value, lli, ssi

        except ValueError, IndexError:
            pass

        try:
            # Method 2: Flexible parsing for variable-length records
            slice_trimmed = slice_text.strip()
            if not slice_trimmed:
                return None, None, None

            # Look for a decimal point to identify the numeric value
            if "." in slice_trimmed:
                # Find the main numeric value (supports negative numbers)
                number_match = re.search(r"(-?\d+\.\d+)", slice_trimmed)

                if number_match:
                    value = float(number_match.group(1))

                    # Check for LLI/SSI indicators after the number
                    remaining_part = slice_trimmed[number_match.end() :].strip()
                    lli = None
                    ssi = None

                    # Parse remaining characters as potential LLI/SSI
                    if remaining_part:
                        # Could be just SSI, or LLI followed by SSI
                        if len(remaining_part) == 1:
                            # Just one indicator - assume it's SSI
                            if remaining_part.isdigit():
                                ssi = int(remaining_part)
                        elif len(remaining_part) >= LLI_SSI_PAIR_LEN:
                            # Two or more characters - take last two as LLI, SSI
                            lli_char = remaining_part[-2]
                            ssi_char = remaining_part[-1]

                            if lli_char.isdigit():
                                lli = int(lli_char)
                            if ssi_char.isdigit():
                                ssi = int(ssi_char)

                    return value, lli, ssi

        except ValueError, IndexError:
            pass

        # Method 3: Last resort - try simple float parsing
        try:
            simple_value = float(slice_text.strip())
            return simple_value, None, None
        except ValueError:
            pass

        return None, None, None

    def process_satellite_data(self, s: str) -> Satellite:
        """Process satellite data line into a Satellite object with observations.

        Handles variable-length observation records correctly by adaptively parsing
        based on the actual line length and content.
        """
        sv = s[:3].strip()
        satellite = Satellite(sv=sv)
        bands_tbe = [f"{sv}|{b}" for b in self.header.obs_codes_per_system[sv[0]]]

        # Get the data part (after sv identifier)
        data_part = s[3:]

        # Process each observation adaptively
        for i, band in enumerate(bands_tbe):
            start_idx = i * 16
            end_idx = start_idx + 16

            # Check if we have enough data for this observation
            if start_idx >= len(data_part):
                # No more data available - create empty observation
                observation = Observation(
                    obs_type=band.split("|")[1][0],
                    value=None,
                    lli=None,
                    ssi=None,
                )
                satellite.add_observation(observation)
                continue

            # Extract the slice, but handle variable length
            if end_idx <= len(data_part):
                # Full 16-character slice available
                slice_data = data_part[start_idx:end_idx]
            else:
                # Partial slice - pad with spaces to maintain consistency
                available_slice = data_part[start_idx:]
                slice_data = available_slice.ljust(16)  # Pad with spaces if needed

            value, lli, ssi = self.parse_observation_slice(slice_data)

            observation = Observation(
                obs_type=band.split("|")[1][0],
                value=value,
                lli=lli,
                ssi=ssi,
            )
            satellite.add_observation(observation)

        return satellite

    @property
    def epochs(self) -> list[Rnxv3ObsEpochRecord]:
        """Materialize all epochs (legacy compatibility).

        Returns
        -------
        list of Rnxv3ObsEpochRecord
            All epochs in memory (use iter_epochs for efficiency)

        """
        return list(self.iter_epochs())

    def iter_epochs(self) -> Iterator[Rnxv3ObsEpochRecord]:
        """Yield epochs one by one instead of materializing the whole list.

        Returns
        -------
        Generator
            Generator yielding Rnxv3ObsEpochRecord objects

        Yields
        ------
        Rnxv3ObsEpochRecord
            Each epoch with timestamp and satellite observations

        """
        for start, end in self.get_epoch_record_batches():
            try:
                info = Rnxv3ObsEpochRecordLineModel.model_validate(
                    {"epoch": self._lines[start]}
                )

                # Skip event epochs (flag 2-6: special records, not observations)
                if info.epoch_flag > 1:
                    continue

                # Filter out blank/whitespace-only lines from data slice
                data = [line for line in self._lines[start + 1 : end] if line.strip()]
                epoch = Rnxv3ObsEpochRecord(
                    info=info,
                    data=[self.process_satellite_data(line) for line in data],
                )
                yield epoch
            except InvalidEpochError, IncompleteEpochError, ValueError:
                # Skip epochs with validation errors (invalid SV, malformed data,
                # pydantic ValidationError inherits from ValueError)
                pass

    def iter_epochs_in_range(
        self,
        start: datetime,
        end: datetime,
    ) -> Iterable[Rnxv3ObsEpochRecord]:
        """Yield epochs lazily that fall into the given datetime range.

        Parameters
        ----------
        start : datetime
            Start of time range (inclusive)
        end : datetime
            End of time range (inclusive)

        Returns
        -------
        Generator
            Generator yielding epochs in the specified range

        Yields
        ------
        Rnxv3ObsEpochRecord
            Epochs within the time range

        """
        for epoch in self.iter_epochs():
            dt = self.get_datetime_from_epoch_record_info(epoch.info)
            if start <= dt <= end:
                yield epoch

    def get_datetime_from_epoch_record_info(
        self,
        epoch_record_info: Rnxv3ObsEpochRecordLineModel,
    ) -> datetime:
        """Convert epoch record info to datetime object.

        Parameters
        ----------
        epoch_record_info : Rnxv3ObsEpochRecordLineModel
            Parsed epoch record line

        Returns
        -------
        datetime
            Timestamp from epoch record

        """
        return datetime(
            year=int(epoch_record_info.year),
            month=int(epoch_record_info.month),
            day=int(epoch_record_info.day),
            hour=int(epoch_record_info.hour),
            minute=int(epoch_record_info.minute),
            second=int(epoch_record_info.seconds),
            tzinfo=UTC,
        )

    @staticmethod
    def epochrecordinfo_dt_to_numpy_dt(
        epch: Rnxv3ObsEpochRecord,
    ) -> np.datetime64:
        """Convert Python datetime to numpy datetime64[ns].

        Parameters
        ----------
        epch : Rnxv3ObsEpochRecord
            Epoch record containing timestamp info

        Returns
        -------
        np.datetime64
            Numpy datetime64 with nanosecond precision

        """
        dt = datetime(
            year=int(epch.info.year),
            month=int(epch.info.month),
            day=int(epch.info.day),
            hour=int(epch.info.hour),
            minute=int(epch.info.minute),
            second=int(epch.info.seconds),
            tzinfo=UTC,
        )
        # np.datetime64 doesn't support timezone info, but datetime is already UTC
        # Convert to naive datetime (UTC) to avoid warning
        return np.datetime64(dt.replace(tzinfo=None), "ns")

    def _epoch_datetimes(self) -> list[datetime]:
        """Extract epoch datetimes from the file.

        Uses the same epoch parsing logic already implemented.
        """
        dts: list[datetime] = []

        for start, _end in self.get_epoch_record_batches():
            info = Rnxv3ObsEpochRecordLineModel.model_validate(
                {"epoch": self._lines[start]}
            )
            dts.append(
                datetime(
                    year=int(info.year),
                    month=int(info.month),
                    day=int(info.day),
                    hour=int(info.hour),
                    minute=int(info.minute),
                    second=int(info.seconds),
                    tzinfo=UTC,
                )
            )
        return dts

    def infer_sampling_interval(self) -> pint.Quantity | None:
        """Infer sampling interval from consecutive epoch deltas.

        Returns
        -------
        pint.Quantity or None
            Sampling interval in seconds, or None if cannot be inferred

        """
        dts = self._epoch_datetimes()
        if len(dts) < MIN_EPOCHS_FOR_INTERVAL:
            return None
        # Compute deltas
        deltas: list[timedelta] = [b - a for a, b in pairwise(dts) if b >= a]
        if not deltas:
            return None
        # Pick the most common delta (robust to an occasional missing epoch)
        seconds = Counter(
            int(dt.total_seconds()) for dt in deltas if dt.total_seconds() > 0
        )
        if not seconds:
            return None
        mode_seconds, _ = seconds.most_common(1)[0]
        return (mode_seconds * UREG.second).to(UREG.seconds)

    def infer_dump_interval(
        self, sampling_interval: pint.Quantity | None = None
    ) -> pint.Quantity | None:
        """Infer the intended dump interval for the RINEX file.

        Parameters
        ----------
        sampling_interval : pint.Quantity, optional
            Known sampling interval. If provided, returns (#epochs * sampling_interval)

        Returns
        -------
        pint.Quantity or None
            Dump interval in seconds, or None if cannot be inferred

        """
        idx = self.get_epoch_record_batches()
        n_epochs = len(idx)
        if n_epochs == 0:
            return None

        if sampling_interval is not None:
            return (n_epochs * sampling_interval).to(UREG.seconds)

        # Fallback: time coverage inclusive (last - first) + typical step
        dts = self._epoch_datetimes()
        if len(dts) == 0:
            return None
        if len(dts) == 1:
            # single epoch: treat as 1 * unknown step (cannot infer)
            return None

        # Estimate step from data
        est_step = self.infer_sampling_interval()
        if est_step is None:
            return None

        # Inclusive coverage often equals (n_epochs - 1) * step; intended
        # dump interval is n_epochs * step.
        return (n_epochs * est_step.to(UREG.seconds)).to(UREG.seconds)

    def validate_epoch_completeness(
        self,
        dump_interval: str | pint.Quantity | None = None,
        sampling_interval: str | pint.Quantity | None = None,
    ) -> None:
        """Validate that the number of epochs matches the expected dump interval.

        Parameters
        ----------
        dump_interval : str or pint.Quantity, optional
            Expected file dump interval. If None, inferred from epochs.
        sampling_interval : str or pint.Quantity, optional
            Expected sampling interval. If None, inferred from epochs.

        Returns
        -------
        None

        Raises
        ------
        MissingEpochError
            If total sampling time doesn't match dump interval
        ValueError
            If intervals cannot be inferred

        """
        # Normalize/Infer sampling interval
        if sampling_interval is None:
            inferred = self.infer_sampling_interval()
            if inferred is None:
                msg = "Could not infer sampling interval from epochs"
                raise ValueError(msg)
            sampling_interval = inferred
        # normalize to pint
        elif not isinstance(sampling_interval, pint.Quantity):
            sampling_interval = UREG.Quantity(sampling_interval).to(UREG.seconds)

        # Normalize/Infer dump interval
        if dump_interval is None:
            inferred_dump = self.infer_dump_interval(
                sampling_interval=sampling_interval
            )
            if inferred_dump is None:
                msg = "Could not infer dump interval from file"
                raise ValueError(msg)
            dump_interval = inferred_dump
        elif not isinstance(dump_interval, pint.Quantity):
            # Accept '15 min', '1h', etc.
            dump_interval = UREG.Quantity(dump_interval).to(UREG.seconds)

        # Build inputs for the validator model
        epoch_indices = self.get_epoch_record_batches()

        # This throws MissingEpochError automatically if inconsistent
        cast(Any, Rnxv3ObsEpochRecordCompletenessModel)(
            epoch_records_indeces=epoch_indices,
            rnx_file_dump_interval=dump_interval,
            sampling_interval=sampling_interval,
        )

    def filter_by_overlapping_groups(
        self,
        ds: xr.Dataset,
        group_preference: dict[str, str] | None = None,
    ) -> xr.Dataset:
        """Filter overlapping bands using per-group preferences.

        Parameters
        ----------
        ds : xr.Dataset
            Dataset with `sid` dimension and signal properties.
        group_preference : dict[str, str], optional
            Mapping of overlap group to preferred band.

        Returns
        -------
        xr.Dataset
            Dataset filtered to preferred overlapping bands.

        """
        if group_preference is None:
            group_preference = {
                "L1_E1_B1I": "L1",
                "L5_E5a": "L5",
                "L2_E5b_B2b": "L2",
            }

        keep = []
        for sid in ds.sid.values:
            parts = str(sid).split("|")
            band = parts[1] if len(parts) >= 2 else ""
            group = self._signal_mapper.get_overlapping_group(band)
            if group and group in group_preference:
                if band == group_preference[group]:
                    keep.append(sid)
            else:
                keep.append(sid)
        return ds.sel(sid=keep)

    def _precompute_sids_from_header(
        self,
    ) -> tuple[list[str], dict[str, dict[str, object]]]:
        """Build sorted SID list and properties from header info alone.

        Uses the header's obs_codes_per_system and static constellation
        SV lists to pre-compute the full theoretical SID set, eliminating
        the discovery pass.

        Returns
        -------
        sorted_sids : list[str]
            Sorted list of signal IDs.
        sid_properties : dict[str, dict[str, object]]
            Mapping of SID to its properties (sv, system, band, code,
            freq_center, freq_min, freq_max, bandwidth, overlapping_group).

        """
        mapper = self._signal_mapper
        signal_ids: set[str] = set()
        sid_properties: dict[str, dict[str, object]] = {}

        # Pre-compute pint arithmetic once per unique band
        band_freq_cache: dict[str, tuple[float, float, float, float]] = {}

        for system, obs_codes in self.header.obs_codes_per_system.items():
            svs = _get_constellation_svs(system)

            for obs_code in obs_codes:
                if len(obs_code) < 3:
                    continue
                band_num = obs_code[1]
                code_char = obs_code[2]

                band_name = mapper.SYSTEM_BANDS.get(system, {}).get(
                    band_num, f"UnknownBand{band_num}"
                )

                # Cache frequency arithmetic per band
                if band_name not in band_freq_cache:
                    center_frequency = mapper.get_band_frequency(band_name)
                    bandwidth = mapper.get_band_bandwidth(band_name)

                    if center_frequency is not None and bandwidth is not None:
                        bw = bandwidth[0] if isinstance(bandwidth, list) else bandwidth
                        freq_min = center_frequency - (bw / 2.0)
                        freq_max = center_frequency + (bw / 2.0)
                        band_freq_cache[band_name] = (
                            float(center_frequency),
                            float(freq_min),
                            float(freq_max),
                            float(bw),
                        )
                    else:
                        band_freq_cache[band_name] = (
                            np.nan,
                            np.nan,
                            np.nan,
                            np.nan,
                        )

                freq_center, freq_min, freq_max, bw = band_freq_cache[band_name]
                overlapping_group = mapper.get_overlapping_group(band_name)

                sid_suffix = "|" + band_name + "|" + code_char

                for sv in svs:
                    sid = sv + sid_suffix
                    if sid not in signal_ids:
                        signal_ids.add(sid)
                        sid_properties[sid] = {
                            "sv": sv,
                            "system": system,
                            "band": band_name,
                            "code": code_char,
                            "freq_center": freq_center,
                            "freq_min": freq_min,
                            "freq_max": freq_max,
                            "bandwidth": bw,
                            "overlapping_group": overlapping_group,
                        }

        sorted_sids = sorted(signal_ids)
        return sorted_sids, {s: sid_properties[s] for s in sorted_sids}

    def _create_dataset_single_pass(self) -> xr.Dataset:
        """Create xarray Dataset in a single pass over the file.

        Pre-allocates arrays using header-derived SID set and epoch count,
        then fills them by parsing observations inline without Pydantic
        models or function-call overhead.

        Returns
        -------
        xr.Dataset
            Dataset with dimensions (epoch, sid) and standard variables.

        """
        lines = self._load_file()
        epoch_batches = self.get_epoch_record_batches()
        n_epochs = len(epoch_batches)

        sorted_sids, sid_properties = self._precompute_sids_from_header()
        n_sids = len(sorted_sids)
        sid_to_idx = {sid: i for i, sid in enumerate(sorted_sids)}

        # Pre-allocate arrays
        timestamps = np.empty(n_epochs, dtype="datetime64[ns]")
        snr = np.full((n_epochs, n_sids), np.nan, dtype=DTYPES["SNR"])
        pseudo = np.full((n_epochs, n_sids), np.nan, dtype=DTYPES["Pseudorange"])
        phase = np.full((n_epochs, n_sids), np.nan, dtype=DTYPES["Phase"])
        doppler = np.full((n_epochs, n_sids), np.nan, dtype=DTYPES["Doppler"])
        lli = np.full((n_epochs, n_sids), -1, dtype=DTYPES["LLI"])
        ssi = np.full((n_epochs, n_sids), -1, dtype=DTYPES["SSI"])

        # Build obs_code → (obs_type, sid_suffix) lookup per system
        mapper = self._signal_mapper
        system_obs_lut: dict[str, list[tuple[str, str]]] = {}
        for system, obs_codes in self.header.obs_codes_per_system.items():
            lut: list[tuple[str, str]] = []
            for obs_code in obs_codes:
                if len(obs_code) < 3:
                    lut.append(("", ""))
                    continue
                band_num = obs_code[1]
                code_char = obs_code[2]
                band_name = mapper.SYSTEM_BANDS.get(system, {}).get(
                    band_num, f"UnknownBand{band_num}"
                )
                obs_type = obs_code[0]
                lut.append((obs_type, "|" + band_name + "|" + code_char))
            system_obs_lut[system] = lut

        # Single pass over all epochs — skip unparseable epoch lines
        valid_mask = np.ones(n_epochs, dtype=bool)
        for t_idx, (start, end) in enumerate(epoch_batches):
            epoch_line = lines[start]

            # Inline epoch parsing (no Pydantic model)
            m = _EPOCH_RE.match(epoch_line)
            if m is None:
                valid_mask[t_idx] = False
                continue

            year, month, day = int(m[1]), int(m[2]), int(m[3])
            hour, minute = int(m[4]), int(m[5])
            seconds = float(m[6])
            sec_int = int(seconds)
            usec = int((seconds - sec_int) * 1_000_000)
            ts = np.datetime64(
                f"{year:04d}-{month:02d}-{day:02d}"
                f"T{hour:02d}:{minute:02d}:{sec_int:02d}",
                "ns",
            )
            ts += np.timedelta64(usec, "us")
            timestamps[t_idx] = ts

            # Parse satellite data lines inline
            for line_idx in range(start + 1, end):
                sat_line = lines[line_idx]
                if len(sat_line) < 3:
                    continue
                sv = sat_line[:3].strip()
                if not sv:
                    continue
                system = sv[0]
                lut_list = system_obs_lut.get(system)
                if lut_list is None:
                    continue

                data_part = sat_line[3:]
                data_part_len = len(data_part)

                for i, (obs_type, sid_suffix) in enumerate(lut_list):
                    if not obs_type:
                        continue

                    col_start = i * 16
                    if col_start >= data_part_len:
                        break

                    sid_key = sv + sid_suffix
                    s_idx = sid_to_idx.get(sid_key)
                    if s_idx is None:
                        continue

                    col_end = col_start + 16
                    slice_text = data_part[col_start:col_end]

                    value, obs_lli, obs_ssi = _parse_obs_fast(slice_text)
                    if value is None:
                        continue

                    if obs_type == "S":
                        if value != 0:
                            snr[t_idx, s_idx] = value
                    elif obs_type == "C":
                        pseudo[t_idx, s_idx] = value
                    elif obs_type == "L":
                        phase[t_idx, s_idx] = value
                    elif obs_type == "D":
                        doppler[t_idx, s_idx] = value

                    if obs_lli is not None:
                        lli[t_idx, s_idx] = obs_lli
                    if obs_ssi is not None:
                        ssi[t_idx, s_idx] = obs_ssi

        # Drop epochs that failed to parse
        if not valid_mask.all():
            timestamps = timestamps[valid_mask]
            snr = snr[valid_mask]
            pseudo = pseudo[valid_mask]
            phase = phase[valid_mask]
            doppler = doppler[valid_mask]
            lli = lli[valid_mask]
            ssi = ssi[valid_mask]

        # Build coordinate arrays from pre-computed properties
        sv_list = np.array(
            [sid_properties[sid]["sv"] for sid in sorted_sids], dtype=object
        )
        constellation_list = np.array(
            [sid_properties[sid]["system"] for sid in sorted_sids], dtype=object
        )
        band_list = np.array(
            [sid_properties[sid]["band"] for sid in sorted_sids], dtype=object
        )
        code_list = np.array(
            [sid_properties[sid]["code"] for sid in sorted_sids], dtype=object
        )
        freq_center_list = [sid_properties[sid]["freq_center"] for sid in sorted_sids]
        freq_min_list = [sid_properties[sid]["freq_min"] for sid in sorted_sids]
        freq_max_list = [sid_properties[sid]["freq_max"] for sid in sorted_sids]

        signal_id_coord = xr.DataArray(
            np.array(sorted_sids, dtype=object),
            dims=["sid"],
            attrs=COORDS_METADATA["sid"],
        )
        coords = {
            "epoch": ("epoch", timestamps, COORDS_METADATA["epoch"]),
            "sid": signal_id_coord,
            "sv": ("sid", sv_list, COORDS_METADATA["sv"]),
            "system": ("sid", constellation_list, COORDS_METADATA["system"]),
            "band": ("sid", band_list, COORDS_METADATA["band"]),
            "code": ("sid", code_list, COORDS_METADATA["code"]),
            "freq_center": (
                "sid",
                np.asarray(freq_center_list, dtype=DTYPES["freq_center"]),
                COORDS_METADATA["freq_center"],
            ),
            "freq_min": (
                "sid",
                np.asarray(freq_min_list, dtype=DTYPES["freq_min"]),
                COORDS_METADATA["freq_min"],
            ),
            "freq_max": (
                "sid",
                np.asarray(freq_max_list, dtype=DTYPES["freq_max"]),
                COORDS_METADATA["freq_max"],
            ),
        }

        if self.header.signal_strength_unit == UREG.dBHz:
            snr_meta = CN0_METADATA
        else:
            snr_meta = SNR_METADATA

        ds = xr.Dataset(
            data_vars={
                "SNR": (["epoch", "sid"], snr, snr_meta),
                "Pseudorange": (
                    ["epoch", "sid"],
                    pseudo,
                    OBSERVABLES_METADATA["Pseudorange"],
                ),
                "Phase": (
                    ["epoch", "sid"],
                    phase,
                    OBSERVABLES_METADATA["Phase"],
                ),
                "Doppler": (
                    ["epoch", "sid"],
                    doppler,
                    OBSERVABLES_METADATA["Doppler"],
                ),
                "LLI": (
                    ["epoch", "sid"],
                    lli,
                    OBSERVABLES_METADATA["LLI"],
                ),
                "SSI": (
                    ["epoch", "sid"],
                    ssi,
                    OBSERVABLES_METADATA["SSI"],
                ),
            },
            coords=coords,
            attrs={**self._build_attrs()},
        )

        if self.apply_overlap_filter:
            ds = self.filter_by_overlapping_groups(ds, self.overlap_preferences)

        return ds

    def create_rinex_netcdf_with_signal_id(
        self,
        start: datetime | None = None,
        end: datetime | None = None,
    ) -> xr.Dataset:
        """Create a NetCDF dataset with signal IDs.

        Always uses the fast single-pass path.  Optionally restricts to
        epochs within a datetime range via post-filtering.

        Parameters
        ----------
        start : datetime, optional
            Start of time range (inclusive).
        end : datetime, optional
            End of time range (inclusive).

        Returns
        -------
        xr.Dataset
            Dataset with dimensions (epoch, sid).

        """
        ds = self._create_dataset_single_pass()

        if start or end:
            ds = ds.sel(epoch=slice(start, end))

        return ds

    def to_ds(
        self,
        keep_data_vars: list[str] | None = None,
        **kwargs: object,
    ) -> xr.Dataset:
        """Convert RINEX observations to xarray.Dataset with signal ID structure.

        Parameters
        ----------
        outname : Path or str, optional
            If provided, saves dataset to this file path
        keep_data_vars : list of str or None, optional
            Data variables to include in dataset. Defaults to config value.
        write_global_attrs : bool, default False
            If True, adds comprehensive global attributes
        pad_global_sid : bool, default True
            If True, pads to global signal ID space
        strip_fillval : bool, default True
            If True, removes fill values
        add_future_datavars : bool, default True
            If True, adds placeholder variables for future data
        keep_sids : list of str or None, default None
            If provided, filters/pads dataset to these specific SIDs.
            If None and pad_global_sid=True, pads to all possible SIDs.

        Returns
        -------
        xr.Dataset
            Dataset with dimensions (epoch, sid) and requested data variables

        """
        outname = cast(Path | str | None, kwargs.pop("outname", None))
        write_global_attrs = bool(kwargs.pop("write_global_attrs", False))
        pad_global_sid = bool(kwargs.pop("pad_global_sid", True))
        strip_fillval = bool(kwargs.pop("strip_fillval", True))
        add_future_datavars = bool(kwargs.pop("add_future_datavars", True))
        keep_sids = cast(list[str] | None, kwargs.pop("keep_sids", None))

        if keep_data_vars is None:
            from canvod.utils.config import load_config

            keep_data_vars = load_config().processing.processing.keep_rnx_vars

        ds = self.create_rinex_netcdf_with_signal_id()

        # drop unwanted vars
        for var in list(ds.data_vars):
            if var not in keep_data_vars:
                ds = ds.drop_vars([var])

        if pad_global_sid:
            from canvod.auxiliary.preprocessing import pad_to_global_sid

            # Pad/filter to specified sids or all possible sids
            ds = pad_to_global_sid(ds, keep_sids=keep_sids)

        if strip_fillval:
            from canvod.auxiliary.preprocessing import strip_fillvalue

            ds = strip_fillvalue(ds)

        if add_future_datavars:
            pass

        if write_global_attrs:
            ds.attrs.update(self._create_comprehensive_attrs())

        ds.attrs.update(self._build_attrs())

        if outname:
            from canvod.utils.config import load_config as _load_config

            comp = _load_config().processing.compression
            encoding = {
                var: {"zlib": comp.zlib, "complevel": comp.complevel}
                for var in ds.data_vars
            }
            ds.to_netcdf(str(outname), encoding=encoding)

        # Validate output structure for pipeline compatibility
        validate_dataset(ds, required_vars=keep_data_vars)

        return ds

    def validate_rinex_304_compliance(
        self,
        ds: xr.Dataset | None = None,
        strict: bool = False,
        print_report: bool = True,
    ) -> dict[str, list[str]]:
        """Run enhanced RINEX 3.04 specification validation.

        Validates:
        1. System-specific observation codes
        2. GLONASS mandatory fields (slot/frequency, biases)
        3. Phase shift records (RINEX 3.01+)
        4. Observation value ranges

        Parameters
        ----------
        ds : xr.Dataset, optional
            Dataset to validate. If None, creates one from current file.
        strict : bool
            If True, raise ValueError on validation failures
        print_report : bool
            If True, print validation report to console

        Returns
        -------
        dict[str, list[str]]
            Validation results by category

        Examples
        --------
        >>> reader = Rnxv3Obs(fpath="station.24o")
        >>> results = reader.validate_rinex_304_compliance()
        >>> # Or validate a specific dataset
        >>> ds = reader.to_ds()
        >>> results = reader.validate_rinex_304_compliance(ds=ds)

        """
        if ds is None:
            ds = self.to_ds(write_global_attrs=False)

        # Prepare header dict for validators
        header_dict: dict[str, Any] = {
            "obs_codes_per_system": self.header.obs_codes_per_system,
        }

        # Add GLONASS-specific headers if available
        if hasattr(self.header, "glonass_slot_frq"):
            header_dict["GLONASS SLOT / FRQ #"] = self.header.glonass_slot_frq

        if hasattr(self.header, "glonass_cod_phs_bis"):
            header_dict["GLONASS COD/PHS/BIS"] = self.header.glonass_cod_phs_bis

        if hasattr(self.header, "phase_shift"):
            header_dict["SYS / PHASE SHIFT"] = self.header.phase_shift

        # Run validation
        results = RINEX304ComplianceValidator.validate_all(
            ds=ds, header_dict=header_dict, strict=strict
        )

        if print_report:
            RINEX304ComplianceValidator.print_validation_report(results)

        return results

    def _create_comprehensive_attrs(self) -> dict[str, object]:
        attrs: dict[str, object] = {
            "File Path": str(self.fpath),
            "File Type": self.header.filetype,
            "RINEX Version": self.header.version,
            "RINEX Type": self.header.rinextype,
            "Observer": self.header.observer,
            "Agency": self.header.agency,
            "Date": self.header.date.isoformat(),
            "Marker Name": self.header.marker_name,
            "Marker Number": self.header.marker_number,
            "Marker Type": self.header.marker_type,
            "Approximate Position": (
                f"(X = {self.header.approx_position[0].magnitude} "
                f"{self.header.approx_position[0].units:~}, "
                f"Y = {self.header.approx_position[1].magnitude} "
                f"{self.header.approx_position[1].units:~}, "
                f"Z = {self.header.approx_position[2].magnitude} "
                f"{self.header.approx_position[2].units:~})"
            ),
            "Receiver Type": self.header.receiver_type,
            "Receiver Version": self.header.receiver_version,
            "Receiver Number": self.header.receiver_number,
            "Antenna Type": self.header.antenna_type,
            "Antenna Number": self.header.antenna_number,
            "Antenna Position": (
                f"(X = {self.header.antenna_position[0].magnitude} "
                f"{self.header.antenna_position[0].units:~}, "
                f"Y = {self.header.antenna_position[1].magnitude} "
                f"{self.header.antenna_position[1].units:~}, "
                f"Z = {self.header.antenna_position[2].magnitude} "
                f"{self.header.antenna_position[2].units:~})"
            ),
            "Program": self.header.pgm,
            "Run By": self.header.run_by,
            "Time of First Observation": json.dumps(
                {k: v.isoformat() for k, v in self.header.t0.items()}
            ),
            "GLONASS COD": self.header.glonass_cod,
            "GLONASS PHS": self.header.glonass_phs,
            "GLONASS BIS": self.header.glonass_bis,
            "GLONASS Slot Frequency Dict": json.dumps(
                self.header.glonass_slot_freq_dict
            ),
            "Leap Seconds": f"{self.header.leap_seconds:~}",
        }
        return attrs