Add cloud optical depth product

cloudnetpy/output.py

@@ -181,7 +181,7 @@ def get_references(identifier: str | None = None, extra: list | None = None) ->
     match identifier:
         case "der":
             references += (
-                ", https://doi.org/10.1175/1520-0426(2002)019<0835:TROSCD>2.0.CO;2"
+                ", https://doi.org/10.1175/1520-0426(2002)019%3C0835:TROSCD%3E2.0.CO;2"
             )
         case "ier":
             references += (
@@ -195,6 +195,14 @@ def get_references(identifier: str | None = None, extra: list | None = None) ->
             references += ", https://doi.org/10.1175/JAM2340.1"
         case "drizzle":
             references += ", https://doi.org/10.1175/JAM-2181.1"
+        case "cod":
+            references += (
+                ", https://doi.org/10.1175/BAMS-88-6-883"
+                ", https://doi.org/10.1175/1520-0426(2002)019%3C0835:TROSCD%3E2.0.CO;2"
+                ", https://doi.org/10.1175/JAM2340.1"
+                ", https://doi.org/10.1175/JAM2543.1"
+                ", https://doi.org/10.5194/amt-13-5335-2020"
+            )
         case "epsilon-radar":
             references += (
                 ", https://doi.org/10.5194/amt-13-5335-2020"
@@ -472,6 +480,7 @@ def _get_identifier(short_id: str) -> str:
         "classification",
         "der",
         "ier",
+        "cod",
         "classification-voodoo",
         "epsilon-radar",
     )
@@ -486,6 +495,8 @@ def _get_identifier(short_id: str) -> str:
         return "ice effective radius"
     if short_id == "der":
         return "droplet effective radius"
+    if short_id == "cod":
+        return "cloud optical depth"
     if short_id == "epsilon-radar":
         return "dissipation rate of turbulent kinetic energy"
     return short_id

cloudnetpy/plotting/plot_meta.py

@@ -170,6 +170,15 @@ class PlotMeta(NamedTuple):
         ("Unquantifiable", _COLORS["seaweed_roll"]),
         ("Undetected", _COLORS["skyblue"]),
     ),
+    "extinction_retrieval_status": (
+        ("_No cloud", _COLORS["white"]),
+        ("Liquid", _COLORS["lightblue"]),
+        ("Liquid, assumed radius", _COLORS["yellow"]),
+        ("Ice", _COLORS["lightsteel"]),
+        ("Ice, corrected atten.", _COLORS["skyblue"]),
+        ("Liquid & ice", _COLORS["darkpurple"]),
+        ("No retrieval", _COLORS["seaweed_roll"]),
+    ),
     "dominant_hydrometeor_type": (
         ("_Clear sky", _COLORS["white"]),
         ("Ice", _COLORS["lightgray"]),
@@ -684,5 +693,40 @@ class PlotMeta(NamedTuple):
         "dominant_hydrometeor_type": PlotMeta(
             clabel=_CLABEL["dominant_hydrometeor_type"],
         ),
+        "extinction_liquid": PlotMeta(
+            cmap="Blues",
+            plot_range=(1e-4, 1e-1),
+            log_scale=True,
+        ),
+        "extinction_ice": PlotMeta(
+            plot_range=(1e-5, 1e-2),
+            log_scale=True,
+        ),
+        "extinction_retrieval_status": PlotMeta(
+            clabel=_CLABEL["extinction_retrieval_status"],
+        ),
+        "extinction_liquid_error": PlotMeta(
+            cmap="RdYlGn_r",
+            plot_range=(0, 5),
+        ),
+        "extinction_ice_error": PlotMeta(
+            cmap="RdYlGn_r",
+            plot_range=(0, 5),
+        ),
+        "optical_depth_error": PlotMeta(
+            zero_line=True,
+        ),
+        "optical_depth": PlotMeta(
+            plot_range=(0.02, 1000),
+            log_scale=True,
+        ),
+        "optical_depth_liquid": PlotMeta(
+            plot_range=(0.02, 1000),
+            log_scale=True,
+        ),
+        "optical_depth_ice": PlotMeta(
+            plot_range=(0.02, 1000),
+            log_scale=True,
+        ),
     },
 }

cloudnetpy/plotting/plotting.py

@@ -904,6 +904,9 @@ def plot(self, figure_data: FigureData, hacky_freq_ind: int | None = None) -> No
         units = self._convert_units()
         if self._plot_meta.mask_zeros:
             self._mask_zeros()
+        is_zero = ~ma.getmaskarray(self._data) & (ma.filled(self._data, 1) == 0)
+        if self._is_log:
+            self._mask_non_positive()
         self._mark_gaps(figure_data)
         self._ax.plot(
             figure_data.time_including_gaps,
@@ -914,14 +917,34 @@ def plot(self, figure_data: FigureData, hacky_freq_ind: int | None = None) -> No
         )
         if self._plot_meta.moving_average:
             self._plot_moving_average(figure_data, hacky_freq_ind)
-        if self._plot_meta.zero_line:
+        if self._plot_meta.zero_line and not self._is_log:
             self._ax.axhline(0, color="black", alpha=0.5, label="_nolegend_")
         self._fill_between_data_gaps(figure_data)
-        self.sub_plot.set_yax(ylabel=units, y_limits=self._get_y_limits())
+        if self._is_log:
+            self._ax.set_yscale("log")
+        y_limits = self._get_y_limits()
+        self.sub_plot.set_yax(ylabel=units, y_limits=y_limits)
+        if self._is_log and np.any(is_zero):
+            self._plot_zeros(figure_data.time[is_zero], y_limits[0])
         pos = self._ax.get_position()
         self._ax.set_position((pos.x0, pos.y0, pos.width * 0.965, pos.height))
         self._plot_flags(figure_data)
 
+    def _plot_zeros(self, time: ndarray, y_min: float) -> None:
+        """Marks zero values along the bottom of a logarithmic axis."""
+        is_cloud_variable = self.sub_plot.variable.name.startswith("optical_depth")
+        self._ax.plot(
+            time,
+            np.full(len(time), y_min * 1.3),
+            color="lightgrey",
+            marker=".",
+            lw=0,
+            markersize=3,
+            label="Clear sky" if is_cloud_variable else "Zero",
+            zorder=_get_zorder("data"),
+        )
+        self._ax.legend(markerscale=3, numpoints=1, frameon=False)
+
     def _plot_flags(self, figure_data: FigureData) -> None:
         if figure_data.is_mwrpy_product():
             flags = self._read_flagged_data(figure_data)
@@ -1021,9 +1044,18 @@ def _add_legend(self, name: str | tuple = ("Flagged data",)) -> None:
             frameon=False,
         )
 
+    def _mask_non_positive(self) -> None:
+        self._data = ma.masked_less_equal(self._data, 0)
+        self._data_orig = ma.masked_less_equal(self._data_orig, 0)
+
     def _get_y_limits(self) -> tuple[float, float]:
         percent_gap = 0.05
         fallback = (-percent_gap, percent_gap)
+        if self._is_log:
+            if self._plot_meta.plot_range is not None:
+                return self._plot_meta.plot_range
+            valid = self._data[~ma.getmaskarray(self._data)]
+            return (valid.min() / 2, valid.max() * 2) if valid.size else (0.1, 10)
         if ma.all(self._data.mask):
             return fallback
         min_data = self._data.min()
@@ -1066,6 +1098,9 @@ def _plot_moving_average(
     ) -> None:
         time = figure_data.time.copy()
         data = self._data_orig.copy()
+        if self._is_log:
+            # Average in log space so large values do not dominate
+            data = ma.log10(ma.masked_less_equal(data, 0))
 
         if figure_data.is_mwrpy_product() or self.sub_plot.variable.name in (
             "tb",
@@ -1097,6 +1132,8 @@ def _plot_moving_average(
                 )
             else:
                 sma = self._calculate_moving_average(data1, time1, window=5)
+                if self._is_log:
+                    sma = 10**sma
             gap_time = _get_max_gap_in_minutes(figure_data)
             gaps = self._find_time_gap_indices(time1, max_gap_min=gap_time) + 1
 

cloudnetpy/products/__init__.py

@@ -1,4 +1,5 @@
 from .classification import generate_classification
+from .cod import generate_cod
 from .der import generate_der
 from .drizzle import generate_drizzle
 from .epsilon_lidar import generate_epsilon_from_lidar