canvod.vod API Reference¶

class TauOmegaZerothOrder(VODCalculator):
    """Calculate VOD using the zeroth-order Tau-Omega approximation.

    Based on Humphrey, V., & Frankenberg, C. (2022).
    """

    def get_delta_snr(self) -> xr.DataArray:
        """Calculate delta SNR = SNR_canopy - SNR_sky.

        Returns
        -------
        xr.DataArray
            Delta SNR in decibels.
        """
        return self.canopy_ds["SNR"] - self.sky_ds["SNR"]

    def decibel2linear(self, delta_snr_db: xr.DataArray) -> xr.DataArray:
        """Convert decibel values to linear values.

        Parameters
        ----------
        delta_snr_db : xr.DataArray
            Delta SNR in decibels.

        Returns
        -------
        xr.DataArray
            Linear-scale values.
        """
        return 10 ** (delta_snr_db / 10)

    def calculate_vod(self) -> xr.Dataset:
        """Calculate VOD using the zeroth-order approximation.

        Returns a lazy xarray.Dataset when the inputs are dask-backed so that
        the caller can rechunk and let a single write trigger computation.
        For in-memory (numpy) inputs the result is computed immediately.

        Returns
        -------
        xr.Dataset
            Dataset containing VOD and angular coordinates.

        Raises
        ------
        ValueError
            If all delta SNR values are NaN (eager arrays only).
        """
        start_time = time.time()
        log.info(
            "vod_calculation_started",
            canopy_epochs=len(self.canopy_ds.epoch),
            sky_epochs=len(self.sky_ds.epoch),
            sids=len(self.canopy_ds.sid),
        )

        delta_snr = self.get_delta_snr()

        # Detect lazy (dask-backed) arrays — avoid .item()/.any()/.all() on
        # large dask arrays, as each call triggers a full compute pass.
        _lazy = delta_snr.chunks is not None

        if not _lazy and delta_snr.isnull().all():
            log.error("vod_calculation_failed", reason="all_delta_snr_nan")
            raise ValueError(
                "All delta_snr values are NaN - check data alignment",
            )

        canopy_transmissivity = self.decibel2linear(delta_snr)

        if not _lazy and (canopy_transmissivity <= 0).any():
            n_invalid = int((canopy_transmissivity <= 0).sum())
            total = canopy_transmissivity.size
            log.warning(
                "invalid_transmissivity",
                invalid_count=n_invalid,
                total_count=total,
                percent=round(100 * n_invalid / total, 2),
            )

        theta = self.canopy_ds["theta"]
        vod = -np.log(canopy_transmissivity) * np.cos(theta)

        vod_ds = xr.Dataset(
            {
                "VOD": vod,
                "delta_snr": delta_snr,
                "phi": self.canopy_ds["phi"],
                "theta": self.canopy_ds["theta"],
            },
            coords=self.canopy_ds.coords,
        )

        duration = time.time() - start_time
        if not _lazy:
            n_valid = int((~vod.isnull()).sum())
            log.info(
                "vod_calculation_complete",
                duration_seconds=round(duration, 2),
                vod_values=vod.size,
                valid_values=n_valid,
                valid_percent=round(100 * n_valid / vod.size, 2),
            )
        else:
            log.info(
                "vod_calculation_complete",
                duration_seconds=round(duration, 2),
                vod_values=vod.size,
                note="lazy graph built; compute deferred to write",
            )

        return vod_ds

class VODCalculator(ABC, BaseModel):
    """Abstract base class for VOD calculation from RINEX store data.

    Notes
    -----
    This is an abstract base class (ABC) and a Pydantic model.
    """

    model_config = ConfigDict(arbitrary_types_allowed=True)

    canopy_ds: xr.Dataset
    sky_ds: xr.Dataset

    @field_validator("canopy_ds", "sky_ds")
    @classmethod
    def validate_datasets(cls, v: xr.Dataset) -> xr.Dataset:
        """Validate dataset inputs.

        Parameters
        ----------
        v : xr.Dataset
            Dataset to validate.

        Returns
        -------
        xr.Dataset
            Validated dataset.

        Raises
        ------
        ValueError
            If the input is not an xarray Dataset or lacks SNR.
        """
        if not isinstance(v, xr.Dataset):
            raise ValueError("Must be xr.Dataset")
        if "SNR" not in v.data_vars:
            raise ValueError("Dataset must contain 'SNR' variable")
        return v

    @abstractmethod
    def calculate_vod(self) -> xr.Dataset:
        """Calculate VOD and return a dataset with VOD, phi, theta.

        Returns
        -------
        xr.Dataset
            Dataset containing VOD and angular coordinates.
        """
        raise NotImplementedError

    @classmethod
    def from_icechunkstore(
        cls,
        icechunk_store_pth: Path,
        canopy_group: str = "canopy_01",
        sky_group: str = "reference_01",
        **open_kwargs: Any,
    ) -> xr.Dataset:
        """Convenience method to calculate VOD directly from an IcechunkStore.

        Parameters
        ----------
        icechunk_store_pth : Path
            Path to Icechunk store.
        canopy_group : str
            Canopy receiver group name.
        sky_group : str
            Sky/reference receiver group name.
        open_kwargs : dict[str, Any]
            Additional keyword arguments for IcechunkStore.open().
            Currently unused.

        Returns
        -------
        xr.Dataset
            VOD dataset.

        Notes
        -----
        Requires canvod-store to be installed.
        """
        try:
            from canvod.store import MyIcechunkStore
        except ImportError as e:
            raise ImportError(
                "canvod-store package required for from_icechunkstore(). "
                "Install with: pip install canvod-store"
            ) from e

        store = MyIcechunkStore(icechunk_store_pth)

        with store.readonly_session() as session:
            canopy_ds = xr.open_zarr(store=session.store, group=canopy_group)
            sky_ds = xr.open_zarr(store=session.store, group=sky_group)

        return cls.from_datasets(
            canopy_ds=canopy_ds,
            sky_ds=sky_ds,
            align=True,
        )

    @classmethod
    def from_datasets(
        cls,
        canopy_ds: xr.Dataset,
        sky_ds: xr.Dataset,
        align: bool = True,
    ) -> xr.Dataset:
        """Convenience method to calculate VOD directly from datasets.

        Parameters
        ----------
        canopy_ds : xr.Dataset
            Canopy receiver dataset.
        sky_ds : xr.Dataset
            Sky/reference receiver dataset.
        align : bool
            Whether to align datasets on common coordinates.

        Returns
        -------
        xr.Dataset
            VOD dataset.
        """
        if align:
            canopy_ds, sky_ds = xr.align(canopy_ds, sky_ds, join="inner")

        calculator = cls(canopy_ds=canopy_ds, sky_ds=sky_ds)
        return calculator.calculate_vod()

class VODCalculator(ABC, BaseModel):
    """Abstract base class for VOD calculation from RINEX store data.

    Notes
    -----
    This is an abstract base class (ABC) and a Pydantic model.
    """

    model_config = ConfigDict(arbitrary_types_allowed=True)

    canopy_ds: xr.Dataset
    sky_ds: xr.Dataset

    @field_validator("canopy_ds", "sky_ds")
    @classmethod
    def validate_datasets(cls, v: xr.Dataset) -> xr.Dataset:
        """Validate dataset inputs.

        Parameters
        ----------
        v : xr.Dataset
            Dataset to validate.

        Returns
        -------
        xr.Dataset
            Validated dataset.

        Raises
        ------
        ValueError
            If the input is not an xarray Dataset or lacks SNR.
        """
        if not isinstance(v, xr.Dataset):
            raise ValueError("Must be xr.Dataset")
        if "SNR" not in v.data_vars:
            raise ValueError("Dataset must contain 'SNR' variable")
        return v

    @abstractmethod
    def calculate_vod(self) -> xr.Dataset:
        """Calculate VOD and return a dataset with VOD, phi, theta.

        Returns
        -------
        xr.Dataset
            Dataset containing VOD and angular coordinates.
        """
        raise NotImplementedError

    @classmethod
    def from_icechunkstore(
        cls,
        icechunk_store_pth: Path,
        canopy_group: str = "canopy_01",
        sky_group: str = "reference_01",
        **open_kwargs: Any,
    ) -> xr.Dataset:
        """Convenience method to calculate VOD directly from an IcechunkStore.

        Parameters
        ----------
        icechunk_store_pth : Path
            Path to Icechunk store.
        canopy_group : str
            Canopy receiver group name.
        sky_group : str
            Sky/reference receiver group name.
        open_kwargs : dict[str, Any]
            Additional keyword arguments for IcechunkStore.open().
            Currently unused.

        Returns
        -------
        xr.Dataset
            VOD dataset.

        Notes
        -----
        Requires canvod-store to be installed.
        """
        try:
            from canvod.store import MyIcechunkStore
        except ImportError as e:
            raise ImportError(
                "canvod-store package required for from_icechunkstore(). "
                "Install with: pip install canvod-store"
            ) from e

        store = MyIcechunkStore(icechunk_store_pth)

        with store.readonly_session() as session:
            canopy_ds = xr.open_zarr(store=session.store, group=canopy_group)
            sky_ds = xr.open_zarr(store=session.store, group=sky_group)

        return cls.from_datasets(
            canopy_ds=canopy_ds,
            sky_ds=sky_ds,
            align=True,
        )

    @classmethod
    def from_datasets(
        cls,
        canopy_ds: xr.Dataset,
        sky_ds: xr.Dataset,
        align: bool = True,
    ) -> xr.Dataset:
        """Convenience method to calculate VOD directly from datasets.

        Parameters
        ----------
        canopy_ds : xr.Dataset
            Canopy receiver dataset.
        sky_ds : xr.Dataset
            Sky/reference receiver dataset.
        align : bool
            Whether to align datasets on common coordinates.

        Returns
        -------
        xr.Dataset
            VOD dataset.
        """
        if align:
            canopy_ds, sky_ds = xr.align(canopy_ds, sky_ds, join="inner")

        calculator = cls(canopy_ds=canopy_ds, sky_ds=sky_ds)
        return calculator.calculate_vod()

class TauOmegaZerothOrder(VODCalculator):
    """Calculate VOD using the zeroth-order Tau-Omega approximation.

    Based on Humphrey, V., & Frankenberg, C. (2022).
    """

    def get_delta_snr(self) -> xr.DataArray:
        """Calculate delta SNR = SNR_canopy - SNR_sky.

        Returns
        -------
        xr.DataArray
            Delta SNR in decibels.
        """
        return self.canopy_ds["SNR"] - self.sky_ds["SNR"]

    def decibel2linear(self, delta_snr_db: xr.DataArray) -> xr.DataArray:
        """Convert decibel values to linear values.

        Parameters
        ----------
        delta_snr_db : xr.DataArray
            Delta SNR in decibels.

        Returns
        -------
        xr.DataArray
            Linear-scale values.
        """
        return 10 ** (delta_snr_db / 10)

    def calculate_vod(self) -> xr.Dataset:
        """Calculate VOD using the zeroth-order approximation.

        Returns a lazy xarray.Dataset when the inputs are dask-backed so that
        the caller can rechunk and let a single write trigger computation.
        For in-memory (numpy) inputs the result is computed immediately.

        Returns
        -------
        xr.Dataset
            Dataset containing VOD and angular coordinates.

        Raises
        ------
        ValueError
            If all delta SNR values are NaN (eager arrays only).
        """
        start_time = time.time()
        log.info(
            "vod_calculation_started",
            canopy_epochs=len(self.canopy_ds.epoch),
            sky_epochs=len(self.sky_ds.epoch),
            sids=len(self.canopy_ds.sid),
        )

        delta_snr = self.get_delta_snr()

        # Detect lazy (dask-backed) arrays — avoid .item()/.any()/.all() on
        # large dask arrays, as each call triggers a full compute pass.
        _lazy = delta_snr.chunks is not None

        if not _lazy and delta_snr.isnull().all():
            log.error("vod_calculation_failed", reason="all_delta_snr_nan")
            raise ValueError(
                "All delta_snr values are NaN - check data alignment",
            )

        canopy_transmissivity = self.decibel2linear(delta_snr)

        if not _lazy and (canopy_transmissivity <= 0).any():
            n_invalid = int((canopy_transmissivity <= 0).sum())
            total = canopy_transmissivity.size
            log.warning(
                "invalid_transmissivity",
                invalid_count=n_invalid,
                total_count=total,
                percent=round(100 * n_invalid / total, 2),
            )

        theta = self.canopy_ds["theta"]
        vod = -np.log(canopy_transmissivity) * np.cos(theta)

        vod_ds = xr.Dataset(
            {
                "VOD": vod,
                "delta_snr": delta_snr,
                "phi": self.canopy_ds["phi"],
                "theta": self.canopy_ds["theta"],
            },
            coords=self.canopy_ds.coords,
        )

        duration = time.time() - start_time
        if not _lazy:
            n_valid = int((~vod.isnull()).sum())
            log.info(
                "vod_calculation_complete",
                duration_seconds=round(duration, 2),
                vod_values=vod.size,
                valid_values=n_valid,
                valid_percent=round(100 * n_valid / vod.size, 2),
            )
        else:
            log.info(
                "vod_calculation_complete",
                duration_seconds=round(duration, 2),
                vod_values=vod.size,
                note="lazy graph built; compute deferred to write",
            )

        return vod_ds