Add doppy

pyproject.toml

@@ -21,7 +21,7 @@ classifiers = [
 dependencies = [
   "cftime",
   "cloudnetpy[extras]>=1.58.0",
-  "halo-reader==0.1.8",
+  "doppy",
   "influxdb-client",
   "requests",
   "rpgpy",

src/data_processing/instrument_process.py

@@ -2,10 +2,10 @@
 import gzip
 import logging
 import os
-import pathlib
 import shutil
 import tempfile
 
+import doppy
 from cloudnetpy.instruments import (
     basta2nc,
     ceilo2nc,
@@ -23,11 +23,10 @@
     ws2nc,
 )
 from cloudnetpy.utils import is_timestamp
-from haloreader.read import read as read_halo
-from haloreader.read import read_bg as read_halobg
 from requests.exceptions import HTTPError
 
 from data_processing import concat_wrapper, nc_header_augmenter, utils
+from data_processing.netcdf_utils import Dataset
 from data_processing.processing_tools import Uuid
 from data_processing.utils import MiscError, RawDataMissingError, fetch_calibration
 
@@ -121,40 +120,20 @@ def process_halo_doppler_lidar(self):
             include_pattern=r"Stare.*\.hpl",
             exclude_tag_subset={"cross"},
         )
+        measurement_date = datetime.date.fromisoformat(self.base.date_str)
         full_paths_bg, _ = self.base.download_instrument(
             self.instrument_pid,
             include_pattern=r"Background.*\.txt",
             exclude_tag_subset={"cross"},
-            date_from=str(
-                datetime.date.fromisoformat(self.base.date_str)
-                - datetime.timedelta(days=30)
-            ),
+            date_from=str(measurement_date - datetime.timedelta(days=1)),
         )
-        full_paths_ = [pathlib.Path(path) for path in full_paths]
-        full_paths_bg_ = [pathlib.Path(path) for path in full_paths_bg]
-        halo = read_halo(full_paths_)
-        halobg = read_halobg(full_paths_bg_)
-        if halo is None:
-            raise RawDataMissingError
-        if halobg is None:
-            raise RawDataMissingError
-        halo.correct_background(halobg)
-        halo.compute_beta()
-        screen = halo.compute_noise_screen()
-        halo.compute_beta_screened(screen)
-        halo.compute_doppler_velocity_screened(screen)
-        halo.convert_time_unit2cloudnet_time()
-        nc_buf = halo.to_nc(
-            nc_map={
-                "variables": {
-                    "beta": "beta_raw",
-                    "beta_screened": "beta",
-                    "doppler_velocity_screened": "v",
-                }
-            },
-            nc_exclude={"variables": {"beta_raw"}},
+        stare = doppy.product.Stare.from_halo_data(
+            data=full_paths,
+            data_bg=full_paths_bg,
+            bg_correction_method=doppy.options.BgCorrectionMethod.FIT,
         )
-        self.temp_file.write(nc_buf)
+
+        _doppy_stare_to_nc(stare, self.temp_file.name)
         data = self._get_payload_for_nc_file_augmenter(self.temp_file.name)
         self.uuid.product = nc_header_augmenter.harmonize_halo_file(data)
 
@@ -474,3 +453,66 @@ def _fix_cl51_timestamps(filename: str, hours: int) -> None:
             lines[ind] = f"-{date_time_utc}\n"
     with open(filename, "w") as file:
         file.writelines(lines)
+
+
+def _doppy_stare_to_nc(stare: doppy.product.Stare, filename: str) -> None:
+    return (
+        Dataset(filename)
+        .add_dimension("time")
+        .add_dimension("range")
+        .add_time(
+            name="time",
+            dimensions=("time",),
+            standard_name="time",
+            long_name="Time UTC",
+            data=stare.time,
+            dtype="f8",
+        )
+        .add_variable(
+            name="range",
+            dimensions=("range",),
+            units="m",
+            data=stare.radial_distance,
+            dtype="f4",
+        )
+        .add_variable(
+            name="elevation",
+            dimensions=("time",),
+            units="degrees",
+            data=stare.elevation,
+            dtype="f4",
+            long_name="elevation from horizontal",
+        )
+        .add_variable(
+            name="beta_raw",
+            dimensions=("time", "range"),
+            units="sr-1 m-1",
+            data=stare.beta,
+            dtype="f4",
+        )
+        .add_variable(
+            name="beta",
+            dimensions=("time", "range"),
+            units="sr-1 m-1",
+            data=stare.beta,
+            dtype="f4",
+            mask=stare.mask,
+        )
+        .add_variable(
+            name="v",
+            dimensions=("time", "range"),
+            units="m s-1",
+            long_name="Doppler velocity",
+            data=stare.radial_velocity,
+            dtype="f4",
+            mask=stare.mask,
+        )
+        .add_scalar_variable(
+            name="wavelength",
+            units="m",
+            standard_name="radiation_wavelength",
+            data=stare.wavelength,
+            dtype="f4",
+        )
+        .add_atribute("doppy_version", doppy.__version__)
+    ).close()

src/data_processing/netcdf_utils.py

@@ -0,0 +1,108 @@
+from __future__ import annotations
+
+from typing import Literal, TypeAlias
+
+import netCDF4
+import numpy as np
+import numpy.typing as npt
+
+NetCDFDataType: TypeAlias = Literal["f4", "f8", "i4", "i8", "u4", "u8"]
+
+
+class Dataset:
+    def __init__(self, filename: str) -> None:
+        self.nc = netCDF4.Dataset(filename, mode="w")
+
+    def add_dimension(self, dim: str) -> Dataset:
+        self.nc.createDimension(dim, None)
+        return self
+
+    def add_atribute(self, key: str, val: str) -> Dataset:
+        setattr(self.nc, key, val)
+        return self
+
+    def add_time(
+        self,
+        name: str,
+        dimensions: tuple[str, ...],
+        data: npt.NDArray[np.datetime64],
+        dtype: NetCDFDataType,
+        standard_name: str | None = None,
+        long_name: str | None = None,
+    ) -> Dataset:
+        time, units, calendar = _convert_time(data)
+        var = self.nc.createVariable(name, dtype, dimensions)
+        var.units = units
+        var.calendar = calendar
+        var[:] = time
+        if standard_name is not None:
+            var.standard_name = standard_name
+        if long_name is not None:
+            var.long_name = long_name
+        return self
+
+    def add_variable(
+        self,
+        name: str,
+        dimensions: tuple[str, ...],
+        units: str,
+        data: npt.NDArray[np.float64],
+        dtype: NetCDFDataType,
+        standard_name: str | None = None,
+        long_name: str | None = None,
+        mask: npt.NDArray[np.bool_] | None = None,
+    ) -> Dataset:
+        var = self.nc.createVariable(name, dtype, dimensions)
+        var.units = units
+        if mask is not None:
+            var[:] = np.ma.masked_array(data, mask)  # type: ignore
+        else:
+            var[:] = data
+        if standard_name is not None:
+            var.standard_name = standard_name
+        if long_name is not None:
+            var.long_name = long_name
+        return self
+
+    def add_scalar_variable(
+        self,
+        name: str,
+        units: str,
+        data: np.float64 | np.int64 | float | int,
+        dtype: NetCDFDataType,
+        standard_name: str | None = None,
+        long_name: str | None = None,
+        mask: npt.NDArray[np.bool_] | None = None,
+    ) -> Dataset:
+        var = self.nc.createVariable(name, dtype)
+        var.units = units
+        var[:] = data
+        if standard_name is not None:
+            var.standard_name = standard_name
+        if long_name is not None:
+            var.long_name = long_name
+        return self
+
+    def close(self) -> None:
+        self.nc.close()
+
+
+def _convert_time(
+    time: npt.NDArray[np.datetime64],
+) -> tuple[npt.NDArray[np.float64], str, str]:
+    """
+    Parameters
+    ----------
+    time : npt.NDArray[np.datetime64["us"]]
+        Must be represented in UTC
+    """
+    if time.dtype != "<M8[us]":
+        raise TypeError("time must be datetime64[us]")
+    MICROSECONDS_TO_HOURS = 1 / (1e6 * 3600)
+    start_of_day = time.min().astype("datetime64[D]")
+    hours_since_start_of_day = (time - start_of_day).astype(
+        np.float64
+    ) * MICROSECONDS_TO_HOURS
+    units = f"hours since {np.datetime_as_string(start_of_day)} 00:00:00 +0000"
+    calendar = "standard"
+    return hours_since_start_of_day, units, calendar

src/data_processing/subcmds/process_cloudnet.py

@@ -11,7 +11,7 @@
 from cloudnetpy.categorize import generate_categorize
 from cloudnetpy.exceptions import CloudnetException, ModelDataError
 from cloudnetpy.utils import date_range
-from haloreader.exceptions import HaloException
+from requests.exceptions import RequestException
 
 from data_processing import instrument_process, processing_tools, utils
 from data_processing.processing_tools import ProcessBase, Uuid
@@ -89,7 +89,6 @@ def process_product(self, product: str):
             RawDataMissingError,
             MiscError,
             NotImplementedError,
-            HaloException,
         ) as err:
             logging.warning(err)
         except CloudnetException as err:

src/data_processing/utils.py

@@ -62,8 +62,8 @@ def create_product_put_payload(
             payload["software"]["cloudnetpy"] = version
         if version := getattr(nc, "mwrpy_version", None):
             payload["software"]["mwrpy"] = version
-        if version := getattr(nc, "haloreader_version", None):
-            payload["software"]["halo-reader"] = version
+        if version := getattr(nc, "doppy_version", None):
+            payload["software"]["doppy"] = version
         if version := getattr(nc, "voodoonet_version", None):
             payload["software"]["voodoonet"] = version
     return payload