Add noise screening for CT25K

Author: siiptuo

README.md

@@ -3,7 +3,10 @@
 [![Run tests](https://github.com/actris-cloudnet/ceilopyter/actions/workflows/test.yml/badge.svg)](https://github.com/actris-cloudnet/ceilopyter/actions/workflows/test.yml)
 [![PyPI version](https://badge.fury.io/py/ceilopyter.svg)](https://badge.fury.io/py/ceilopyter)
 
-Python package for reading ceilometer data.
+Python package for reading ceilometer data and doing some post-processing like
+screening background noise. This package is used in
+[CloudnetPy](https://github.com/actris-cloudnet/cloudnetpy) for generating the
+[Cloudnet lidar product](https://cloudnet.fmi.fi/product/lidar).
 
 ## Supported ceilometers
 

ceilopyter/ceilo.py

@@ -1,12 +1,32 @@
-from .ceilo_raw import CeiloRaw, concatenate_raw
+import numpy as np
+import numpy.typing as npt
+
+from .ceilo_raw import CeiloRaw
 
 
 class Ceilo:
-    def __init__(self, raw: list[CeiloRaw], calibration_factor: float):
-        concat = concatenate_raw(raw)
-        self.time = concat.time
-        self.range = concat.range
-        self.beta_raw = concat.beta * calibration_factor
+    """Raw ceilometer data.
+
+    Attributes:
+        time: Time
+        range: Range (m)
+        beta_raw: Non-screened range-corrected backscatter coefficient (sr-1 m-1)
+        beta: Screened range-corrected backscatter coefficient (sr-1 m-1)
+        wavelength: Wavelength (nm)
+        zenith_angle: Zenith angle (deg)
+    """
+
+    def __init__(
+        self,
+        raw: CeiloRaw,
+        beta_raw: npt.NDArray[np.floating],
+        beta: npt.NDArray[np.floating] | None,
+        calibration_factor: float,
+    ):
+        self.time = raw.time
+        self.range = raw.range
+        self.beta = beta
+        self.beta_raw = beta_raw
         self.calibration_factor = calibration_factor
-        self.wavelength = concat.wavelength
-        self.zenith_angle = concat.zenith_angle
+        self.wavelength = raw.wavelength
+        self.zenith_angle = raw.zenith_angle

ceilopyter/instruments/chm15k.py

@@ -8,7 +8,7 @@
 from numpy import ma
 
 from ..ceilo import Ceilo
-from ..ceilo_raw import CeiloRaw
+from ..ceilo_raw import CeiloRaw, concatenate_raw
 
 
 def read_chm15k(
@@ -23,7 +23,9 @@ def read_chm15k(
     raw = []
     for file in files:
         raw.append(_read_file(file))
-    return Ceilo(raw, calibration_factor)
+    concat = concatenate_raw(raw)
+    beta_raw = concat.beta * calibration_factor
+    return Ceilo(concat, beta_raw, None, calibration_factor)
 
 
 def _read_file(file: str | PathLike) -> CeiloRaw:

ceilopyter/instruments/cl31.py

@@ -2,6 +2,7 @@
 from os import PathLike
 
 from ..ceilo import Ceilo
+from ..ceilo_raw import concatenate_raw
 from ..common import read_msgs
 from ..readers.read_cl import read_cl_file
 
@@ -13,5 +14,7 @@ def read_cl31(
     if calibration_factor is None:
         calibration_factor = 1.0
         logging.warning("Using default calibration factor: %s", calibration_factor)
-    raw = read_msgs(files, read_cl_file, 910.0)
-    return Ceilo(raw, calibration_factor)
+    raw = read_msgs(files, read_cl_file, wavelength=910.0)
+    concat = concatenate_raw(raw)
+    beta_raw = concat.beta * calibration_factor
+    return Ceilo(concat, beta_raw, None, calibration_factor)

ceilopyter/instruments/cl51.py

@@ -2,6 +2,7 @@
 from os import PathLike
 
 from ..ceilo import Ceilo
+from ..ceilo_raw import concatenate_raw
 from ..common import read_msgs
 from ..readers.read_cl import read_cl_file
 
@@ -13,5 +14,7 @@ def read_cl51(
     if calibration_factor is None:
         calibration_factor = 1.0
         logging.warning("Using default calibration factor: %s", calibration_factor)
-    raw = read_msgs(files, read_cl_file, 910.0)
-    return Ceilo(raw, calibration_factor)
+    raw = read_msgs(files, read_cl_file, wavelength=910.0)
+    concat = concatenate_raw(raw)
+    beta_raw = concat.beta * calibration_factor
+    return Ceilo(concat, beta_raw, None, calibration_factor)

ceilopyter/instruments/cl61.py

@@ -6,7 +6,7 @@
 from cftime import num2pydate
 
 from ..ceilo import Ceilo
-from ..ceilo_raw import CeiloRaw
+from ..ceilo_raw import CeiloRaw, concatenate_raw
 
 
 def read_cl61(
@@ -21,7 +21,9 @@ def read_cl61(
             raw = list(executor.map(_read_file, files))
     else:
         raw = [_read_file(files)]
-    return Ceilo(raw, calibration_factor)
+    concat = concatenate_raw(raw)
+    beta_raw = concat.beta * calibration_factor
+    return Ceilo(concat, beta_raw, None, calibration_factor)
 
 
 def _read_file(file: str | PathLike) -> CeiloRaw:

ceilopyter/instruments/cs135.py

@@ -1,6 +1,8 @@
 import logging
 from os import PathLike
 
+from ceilopyter.ceilo_raw import concatenate_raw
+
 from ..ceilo import Ceilo
 from ..common import read_msgs
 from ..readers.read_cs import read_cs_file
@@ -13,5 +15,7 @@ def read_cs135(
     if calibration_factor is None:
         calibration_factor = 1.0
         logging.warning("Using default calibration factor: %s", calibration_factor)
-    raw = read_msgs(files, read_cs_file, 905.0)
-    return Ceilo(raw, calibration_factor)
+    raw = read_msgs(files, read_cs_file, wavelength=905.0)
+    concat = concatenate_raw(raw)
+    beta_raw = concat.beta * calibration_factor
+    return Ceilo(concat, beta_raw, None, calibration_factor)

ceilopyter/instruments/ct25k.py

@@ -1,17 +1,44 @@
 import logging
 from os import PathLike
 
+import numpy as np
+import numpy.typing as npt
+from numpy import ma
+
 from ..ceilo import Ceilo
+from ..ceilo_raw import concatenate_raw
 from ..common import read_msgs
+from ..noise import remove_noise
 from ..readers.read_ct import read_ct_file
 
 
 def read_ct25k(
     files: str | PathLike | list[str | PathLike],
     calibration_factor: float | None = None,
+    noise_h2: bool = True,
 ) -> Ceilo:
     if calibration_factor is None:
         calibration_factor = 1.0
         logging.warning("Using default calibration factor: %s", calibration_factor)
-    raw = read_msgs(files, read_ct_file, 905.0)
-    return Ceilo(raw, calibration_factor)
+
+    raw = read_msgs(files, read_ct_file, wavelength=905.0)
+    concat = concatenate_raw(raw)
+
+    r2 = (concat.range * 1e-3) ** 2
+    beta_calib = concat.beta * calibration_factor
+    beta_uncorr = beta_calib / r2 if noise_h2 else _fix_beta(r2, beta_calib)
+
+    is_noise = remove_noise(beta_uncorr, noise_floor=6e-8)
+    beta_raw = beta_uncorr * r2
+    beta = ma.masked_where(is_noise, beta_raw)
+
+    return Ceilo(concat, beta_raw, beta, calibration_factor)
+
+
+def _fix_beta(
+    r2: npt.NDArray[np.floating], beta_raw: npt.NDArray[np.floating]
+) -> npt.NDArray[np.floating]:
+    range_broad = np.broadcast_to(r2, beta_raw.shape)
+    is_strong = beta_raw > 1e-7
+    beta_raw[is_strong] /= range_broad[is_strong]
+    return beta_raw

ceilopyter/instruments/ld40.py

@@ -5,7 +5,7 @@
 import numpy as np
 
 from ..ceilo import Ceilo
-from ..ceilo_raw import CeiloRaw
+from ..ceilo_raw import CeiloRaw, concatenate_raw
 from ..common import InvalidMessageError
 
 
@@ -21,7 +21,9 @@ def read_ld40(
     raw = []
     for file in files:
         raw.append(_read_file(file))
-    return Ceilo(raw, calibration_factor)
+    concat = concatenate_raw(raw)
+    beta_raw = concat.beta * calibration_factor
+    return Ceilo(concat, beta_raw, None, calibration_factor)
 
 
 def _read_file(filename: str | PathLike) -> CeiloRaw:

ceilopyter/noise.py

@@ -0,0 +1,22 @@
+import numpy as np
+import numpy.typing as npt
+from numpy import ma
+
+
+def remove_noise(
+    beta_uncorr: npt.NDArray[np.floating], noise_floor: float, snr_limit: float = 5
+) -> npt.NDArray[np.bool]:
+    zero_ranges = np.all(beta_uncorr == 0, axis=0)
+    n_zeros = np.argmax(~zero_ranges[::-1])
+
+    fraction = 0.1
+    n_top_gates = round(beta_uncorr.shape[1] * fraction)
+    beta_top = beta_uncorr[:, -n_top_gates - n_zeros : -n_zeros]
+    noise = ma.std(beta_top, axis=1)
+
+    noise = np.maximum(noise, noise_floor)
+    snr = beta_uncorr / noise[:, np.newaxis]
+
+    is_noise = snr < snr_limit
+
+    return is_noise