From c58b1eb04849016a652ef2cc27876fe714156deb Mon Sep 17 00:00:00 2001 From: Simo Tukiainen Date: Thu, 10 Sep 2026 09:26:59 +0300 Subject: [PATCH 1/4] Add cloud optical depth product --- cloudnetpy/cli.py | 4 +- cloudnetpy/output.py | 11 + cloudnetpy/plotting/plot_meta.py | 44 +++ cloudnetpy/plotting/plotting.py | 20 +- cloudnetpy/products/__init__.py | 1 + cloudnetpy/products/optical_depth.py | 484 +++++++++++++++++++++++++++ docs/source/api.rst | 2 + tests/unit/test_optical_depth.py | 396 ++++++++++++++++++++++ 8 files changed, 960 insertions(+), 2 deletions(-) create mode 100644 cloudnetpy/products/optical_depth.py create mode 100644 tests/unit/test_optical_depth.py diff --git a/cloudnetpy/cli.py b/cloudnetpy/cli.py index 429ef855..0311c986 100644 --- a/cloudnetpy/cli.py +++ b/cloudnetpy/cli.py @@ -822,7 +822,9 @@ def _plot_l3( def _process_cat_product(product: str, categorize_file: str) -> str: output_file = categorize_file.replace("categorize", product) module = importlib.import_module("cloudnetpy.products") - getattr(module, f"generate_{product}")(categorize_file, output_file) + getattr(module, f"generate_{product.replace('-', '_')}")( + categorize_file, output_file + ) logging.info("Processed %s: %s", product, output_file) return output_file diff --git a/cloudnetpy/output.py b/cloudnetpy/output.py index 52157e30..f5d63133 100644 --- a/cloudnetpy/output.py +++ b/cloudnetpy/output.py @@ -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 "optical-depth": + references += ( + ", https://doi.org/10.1175/BAMS-88-6-883" + ", https://doi.org/10.1175/1520-0426(2002)019<0835:TROSCD>2.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", + "optical-depth", "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 == "optical-depth": + return "cloud optical depth" if short_id == "epsilon-radar": return "dissipation rate of turbulent kinetic energy" return short_id diff --git a/cloudnetpy/plotting/plot_meta.py b/cloudnetpy/plotting/plot_meta.py index a22b75e8..4e8221dd 100644 --- a/cloudnetpy/plotting/plot_meta.py +++ b/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, + ), }, } diff --git a/cloudnetpy/plotting/plotting.py b/cloudnetpy/plotting/plotting.py index 36cb6d5f..217b615a 100644 --- a/cloudnetpy/plotting/plotting.py +++ b/cloudnetpy/plotting/plotting.py @@ -904,6 +904,8 @@ 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() + if self._is_log: + self._mask_non_positive() self._mark_gaps(figure_data) self._ax.plot( figure_data.time_including_gaps, @@ -914,9 +916,11 @@ 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) + if self._is_log: + self._ax.set_yscale("log") self.sub_plot.set_yax(ylabel=units, y_limits=self._get_y_limits()) pos = self._ax.get_position() self._ax.set_position((pos.x0, pos.y0, pos.width * 0.965, pos.height)) @@ -1021,9 +1025,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 +1079,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 +1113,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 diff --git a/cloudnetpy/products/__init__.py b/cloudnetpy/products/__init__.py index bff43d97..c89e9121 100644 --- a/cloudnetpy/products/__init__.py +++ b/cloudnetpy/products/__init__.py @@ -7,3 +7,4 @@ from .iwc import generate_iwc from .lwc import generate_lwc from .mwr_tools import generate_mwr_lhumpro, generate_mwr_multi, generate_mwr_single +from .optical_depth import generate_optical_depth diff --git a/cloudnetpy/products/optical_depth.py b/cloudnetpy/products/optical_depth.py new file mode 100644 index 00000000..d1079968 --- /dev/null +++ b/cloudnetpy/products/optical_depth.py @@ -0,0 +1,484 @@ +"""Module for creating Cloudnet cloud optical depth product.""" + +from os import PathLike +from uuid import UUID + +import numpy as np +import numpy.typing as npt +from numpy import ma + +from cloudnetpy import constants, output, utils +from cloudnetpy.datasource import DataSource +from cloudnetpy.metadata import MetaData +from cloudnetpy.products.der import DerSource +from cloudnetpy.products.ier import IerSource +from cloudnetpy.products.iwc import IwcSource +from cloudnetpy.products.lwc import CloudAdjustor, Lwc, LwcError, LwcSource +from cloudnetpy.products.product_tools import IceClassification + +# Density of liquid water (kg m-3) +RHO_WATER = 1000 + +# Droplet effective radius (m) assumed when it cannot be retrieved from radar +DEFAULT_ASSUMED_DER = 10e-6 + +# Radar-retrieved droplet effective radius (m) outside this range is rejected +DER_VALID_RANGE = (2e-6, 50e-6) + +# Measured LWP above this multiple of the adiabatic LWP of the detected liquid +# layers indicates liquid in layers not classified as droplets +LWP_RATIO_LIMIT = 2.0 + +# Relative error assumed for the droplet effective radius when not retrieved +ASSUMED_DER_REL_ERROR = 0.4 + + +def generate_optical_depth( + categorize_file: str | PathLike, + output_file: str | PathLike, + uuid: str | UUID | None = None, + assumed_der: float = DEFAULT_ASSUMED_DER, +) -> UUID: + """Generates Cloudnet cloud optical depth product. + + This function calculates the visible extinction coefficient of liquid and + ice clouds and integrates it over the profile to give cloud optical depth. + Liquid extinction is derived from the adiabatic-scaled liquid water content + (see :func:`generate_lwc`) and the LWP-scaled droplet effective radius + (see :func:`generate_der`). Ice extinction is derived from ice water content + and ice effective radius, both retrieved from radar reflectivity and model + temperature (see :func:`generate_iwc` and :func:`generate_ier`). Extinction + is calculated in the geometric optics limit, i.e. as 3/2 times the mass + content divided by the particle density and effective radius. The results + are written in a netCDF file. + + Args: + categorize_file: Categorize file name. + output_file: Output file name. + uuid: Set specific UUID for the file. + assumed_der: Droplet effective radius (m) used in liquid layers where + the radar-based retrieval is not available. + + Returns: + UUID of the generated file. + + Examples: + >>> from cloudnetpy.products import generate_optical_depth + >>> generate_optical_depth('categorize.nc', 'optical_depth.nc') + + References: + Frisch, S., Shupe, M., Djalalova, I., Feingold, G., & Poellot, M. (2002). + The Retrieval of Stratus Cloud Droplet Effective Radius with Cloud Radars, + Journal of Atmospheric and Oceanic Technology, 19(6), 835-842. + https://doi.org/10.1175/1520-0426(2002)019%3C0835:TROSCD%3E2.0.CO;2 + + Hogan, R. J., Mittermaier, M. P., & Illingworth, A. J. (2006). The + Retrieval of Ice Water Content from Radar Reflectivity Factor and + Temperature and Its Use in Evaluating a Mesoscale Model, Journal of + Applied Meteorology and Climatology, 45(2), 301-317. + https://doi.org/10.1175/JAM2340.1 + + Delanoë, J., Protat, A., Bouniol, D., Heymsfield, A., Bansemer, A., & + Brown, P. (2007). The Characterization of Ice Cloud Properties from + Doppler Radar Measurements, Journal of Applied Meteorology and + Climatology, 46(10), 1682-1698. https://doi.org/10.1175/JAM2543.1 + + Griesche, H. J., Seifert, P., Ansmann, A., Baars, H., Barrientos + Velasco, C., Bühl, J., Engelmann, R., Radenz, M., Zhenping, Y., & + Macke, A. (2020): Application of the shipborne remote sensing supersite + OCEANET for profiling of Arctic aerosols and clouds during Polarstern + cruise PS106, Atmos. Meas. Tech., 13, 5335–5358. + https://doi.org/10.5194/amt-13-5335-2020 + + """ + uuid = utils.get_uuid(uuid) + with OpticalDepthSource(categorize_file, assumed_der) as od_source: + od_source.append_liquid_extinction() + od_source.append_ice_extinction() + od_source.append_extinction_status() + od_source.append_optical_depths() + od_source.append_optical_depth_error() + od_source.append_optical_depth_status() + date = od_source.get_date() + attributes = output.add_time_attribute(dict(OPTICAL_DEPTH_ATTRIBUTES), date) + attributes = _add_extinction_comments(attributes, od_source) + output.update_attributes(od_source.data, attributes) + output.save_product_file( + "optical-depth", od_source, output_file, uuid, copy_from_cat=("lwp",) + ) + return uuid + + +class OpticalDepthSource(DataSource): + """Data container for cloud optical depth calculations.""" + + def __init__(self, categorize_file: str | PathLike, assumed_der: float) -> None: + super().__init__(categorize_file) + self.categorize_file = categorize_file + self.assumed_der = assumed_der + self.height_agl: npt.NDArray + self.path_lengths = utils.path_lengths_from_ground(self.height_agl) + self.ice_classification = IceClassification(categorize_file) + self.is_rain = self.ice_classification.is_rain.astype(bool) + self.is_liquid = self.ice_classification.category_bits.droplet + self.has_lwp = "lwp" in self.dataset.variables + self.is_der_assumed = np.zeros(self.is_liquid.shape, dtype=bool) + self.is_lidar_only_ice = np.zeros(self.is_liquid.shape, dtype=bool) + self.is_lwp_inconsistent = np.zeros(self.is_liquid.shape[0], dtype=bool) + self._rel_error: dict[str, ma.MaskedArray] = {} + + def append_liquid_extinction(self) -> None: + """Calculates liquid extinction from LWC and droplet effective radius.""" + lwc, lwc_rel_error, lwp_adiabatic = self._get_lwc() + self.is_lwp_inconsistent = self._find_inconsistent_lwp(lwp_adiabatic) + der, der_rel_error = self._get_der() + der = ma.masked_outside(der, *DER_VALID_RANGE) + self.is_der_assumed = ~ma.getmaskarray(lwc) & ma.getmaskarray(der) + der_filled = ma.filled(der, self.assumed_der) + der_rel_error_filled = ma.filled(der_rel_error, ASSUMED_DER_REL_ERROR).copy() + der_rel_error_filled[self.is_der_assumed] = ASSUMED_DER_REL_ERROR + # Includes pixels added at lidar-only cloud tops by the lwc retrieval + extinction = 3 * lwc / (2 * RHO_WATER * der_filled) + rel_error = utils.l2norm(ma.filled(lwc_rel_error, 0), der_rel_error_filled) + self._rel_error["liquid"] = ma.masked_where( + ma.getmaskarray(extinction), rel_error + ) + self.append_data(extinction, "extinction_liquid") + self.append_data( + _relative_to_db(self._rel_error["liquid"]), "extinction_liquid_error" + ) + + def append_ice_extinction(self) -> None: + """Calculates ice extinction from IWC and ice effective radius.""" + iwc, iwc_error = self._get_iwc() + ier = self._get_ier() + self.is_lidar_only_ice = self.ice_classification.is_ice & ma.getmaskarray(iwc) + extinction = 3 * iwc / (2 * constants.RHO_ICE * ier) + extinction[~self.ice_classification.is_ice] = ma.masked + extinction[self.is_rain, :] = ma.masked + rel_error = utils.db2lin(ma.array(iwc_error, copy=True)) - 1 + self._rel_error["ice"] = ma.masked_where(ma.getmaskarray(extinction), rel_error) + error = ma.array(iwc_error, copy=True) + error[ma.getmaskarray(extinction)] = ma.masked + self.append_data(extinction, "extinction_ice") + self.append_data(error, "extinction_ice_error") + + def append_optical_depth_error(self) -> None: + """Estimates the error of the column optical depth. + + Errors are assumed fully correlated within a profile for each phase + (they are dominated by LWP, effective radius and Z-T relation + uncertainties) and independent between liquid and ice. + """ + tau = self.data["optical_depth"][:] + abs_errors = [ + self._integrate(ma.filled(self._rel_error[phase], 0) * ma.filled(ext, 0)) + for phase, ext in ( + ("liquid", self.data["extinction_liquid"][:]), + ("ice", self.data["extinction_ice"][:]), + ) + ] + abs_error = utils.l2norm(*abs_errors) + error = _relative_to_db(abs_error / ma.masked_less_equal(tau, 0)) + self.append_data(error, "optical_depth_error") + + def append_extinction_status(self) -> None: + """Adds pixel-wise retrieval status.""" + is_liquid = ~ma.getmaskarray(self.data["extinction_liquid"][:]) + is_ice = ~ma.getmaskarray(self.data["extinction_ice"][:]) + status = np.zeros(is_liquid.shape, dtype=int) + status[is_liquid] = 1 + status[is_liquid & self.is_der_assumed] = 2 + status[is_ice] = 3 + status[is_ice & self.ice_classification.corrected_ice] = 4 + status[is_liquid & is_ice] = 5 + status[self._find_cloud_without_retrieval()] = 6 + self.append_data(status, "extinction_retrieval_status") + + def append_optical_depths(self) -> None: + """Integrates extinction over the profile.""" + tau_liquid = self._integrate(self.data["extinction_liquid"][:]) + tau_ice = self._integrate(self.data["extinction_ice"][:]) + tau = tau_liquid + tau_ice + no_retrieval = self._find_profiles_without_retrieval() + for array in (tau_liquid, tau_ice, tau): + array[no_retrieval] = ma.masked + self.append_data(tau_liquid, "optical_depth_liquid") + self.append_data(tau_ice, "optical_depth_ice") + self.append_data(tau, "optical_depth") + + def append_optical_depth_status(self) -> None: + """Adds profile-wise retrieval status.""" + is_cloud = self.is_liquid | self.ice_classification.is_ice + status = np.zeros(is_cloud.shape[0], dtype=int) + status[np.any(is_cloud, axis=1)] = 1 + status[np.any(self.is_der_assumed, axis=1)] = 2 + status[np.any(self.is_lidar_only_ice, axis=1)] = 3 + status[self.is_lwp_inconsistent] = 4 + status[self._find_profiles_without_retrieval()] = 5 + self.append_data(status, "optical_depth_retrieval_status") + + def _integrate(self, extinction: npt.NDArray) -> ma.MaskedArray: + return ma.sum(ma.filled(extinction, 0) * self.path_lengths, axis=1) + + def _find_cloud_without_retrieval(self) -> npt.NDArray: + missing_ice = self.ice_classification.is_ice & self._is_missing( + "extinction_ice" + ) + return self._find_missing_liquid() | missing_ice + + def _find_missing_liquid(self) -> npt.NDArray: + return self.is_liquid & self._is_missing("extinction_liquid") + + def _is_missing(self, key: str) -> npt.NDArray: + return ma.getmaskarray(self.data[key][:]) + + def _find_profiles_without_retrieval(self) -> npt.NDArray: + missing_liquid = self._find_missing_liquid() + uncorrected_ice = self.ice_classification.uncorrected_ice + return ( + self.is_rain + | np.any(missing_liquid, axis=1) + | np.any(uncorrected_ice, axis=1) + ) + + def _find_inconsistent_lwp(self, lwp_adiabatic: npt.NDArray) -> npt.NDArray: + """Finds profiles where LWP exceeds what the detected layers can hold.""" + if not self.has_lwp: + return np.zeros(len(lwp_adiabatic), dtype=bool) + lwp = ma.filled(self.getvar("lwp"), 0) + lwp_adiabatic = ma.filled(lwp_adiabatic, 0) + has_liquid = lwp_adiabatic > 0 + return has_liquid & (lwp > LWP_RATIO_LIMIT * lwp_adiabatic) + + def _get_lwc(self) -> tuple[ma.MaskedArray, npt.NDArray, npt.NDArray]: + """Returns LWC (kg m-3), its relative error and adiabatic LWP (kg m-2). + + Without a microwave radiometer, LWC is masked everywhere and the + liquid layers get no retrieval. + """ + shape = self.is_liquid.shape + if not self.has_lwp: + return ma.masked_all(shape), np.zeros(shape), np.zeros(shape[0]) + with LwcSource(self.categorize_file) as lwc_source: + lwc = Lwc(lwc_source) + status = CloudAdjustor(lwc_source, lwc).status + rel_error = LwcError(lwc_source, lwc).error + lwp_positive = lwc_source.lwp > 0 + lwp_adiabatic = ma.sum(lwc.lwc_adiabatic * self.path_lengths, axis=1) + valid = np.isin(status, (1, 2, 3)) & utils.transpose(lwp_positive) + return ma.masked_where(~valid, lwc.lwc), rel_error, lwp_adiabatic + + def _get_der(self) -> tuple[ma.MaskedArray, ma.MaskedArray]: + """Returns LWP-scaled droplet effective radius (m) and its relative error.""" + if not self.has_lwp: + return ma.masked_all(self.is_liquid.shape), ma.masked_all( + self.is_liquid.shape + ) + with DerSource(self.categorize_file) as der_source: + der_source.append_der() + der = der_source.data["der_scaled"][:] + error = der_source.data["der_scaled_error"][:] + return der, error / der + + def _get_iwc(self) -> tuple[ma.MaskedArray, ma.MaskedArray]: + """Returns IWC (kg m-3) and its random error (dB).""" + with IwcSource(self.categorize_file, "iwc") as iwc_source: + iwc_source.append_icy_data(self.ice_classification) + iwc_source.append_error(self.ice_classification) + return iwc_source.data["iwc"][:], iwc_source.data["iwc_error"][:] + + def _get_ier(self) -> ma.MaskedArray: + with IerSource(self.categorize_file, "ier") as ier_source: + ier_source.append_icy_data(self.ice_classification) + ier_source.convert_units() + return ier_source.data["ier"][:] + + +def _relative_to_db(rel_error: npt.NDArray) -> ma.MaskedArray: + """Converts relative error to dB, i.e. 10 log10(1 + error).""" + return ma.array(utils.lin2db(1 + ma.array(rel_error))) + + +def _add_extinction_comments(attributes: dict, od_source: OpticalDepthSource) -> dict: + comment = attributes["extinction_liquid"].comment.format( + der=od_source.assumed_der * 1e6, + der_min=DER_VALID_RANGE[0] * 1e6, + der_max=DER_VALID_RANGE[1] * 1e6, + ) + attributes["extinction_liquid"] = attributes["extinction_liquid"]._replace( + comment=comment + ) + comment = attributes["extinction_liquid_error"].comment.format( + der_error=ASSUMED_DER_REL_ERROR * 100 + ) + attributes["extinction_liquid_error"] = attributes[ + "extinction_liquid_error" + ]._replace(comment=comment) + return attributes + + +COMMENTS = { + "extinction_liquid": ( + "This variable was calculated for the pixels where the categorization\n" + "data has diagnosed liquid droplets and a reliable liquid water path was\n" + "available from a coincident microwave radiometer. Where the liquid\n" + "layer was not detected by the radar, or the retrieved droplet\n" + "effective radius was outside the range {der_min:.0f}-{der_max:.0f} um\n" + "(e.g. due to drizzle), an assumed effective radius of {der:.0f} um\n" + "was used. Missing values indicate that liquid water path was\n" + "unavailable, unreliable or zero, or that rain was present in the\n" + "profile.\n" + "Note that the liquid water path is distributed over the detected\n" + "liquid layers only. If liquid is present in layers not classified as\n" + "droplets (e.g. mixed-phase cloud above the lidar-detected liquid\n" + "base), the extinction of the detected layers is overestimated. Such\n" + "profiles are flagged in the optical_depth_retrieval_status variable." + ), + "extinction_ice": ( + "This variable was calculated for the pixels where the categorization\n" + "data has diagnosed that the radar echo is due to ice. Missing values\n" + "indicate that ice was detected only by the lidar, or that the radar\n" + "reflectivity was affected by uncorrected liquid, rain or melting\n" + "attenuation." + ), + "extinction_liquid_error": ( + "Random error in liquid extinction, one standard deviation, expressed\n" + "as 10 log10(1 + relative error). It combines the liquid water content\n" + "error (liquid water path error and cloud boundary uncertainty, see the\n" + "lwc product) with the droplet effective radius error from the Frisch\n" + "method, or an assumed {der_error:.0f} % where the effective radius\n" + "was assumed." + ), + "extinction_ice_error": ( + "Random error in ice extinction, one standard deviation, expressed as\n" + "10 log10(1 + relative error). The error of the ice water content\n" + "retrieval (see the iwc product) is used, as the empirical\n" + "reflectivity-temperature relation for extinction has a similar\n" + "uncertainty. It includes the additional uncertainty of the liquid\n" + "attenuation correction where applied." + ), + "optical_depth_error": ( + "Random error in cloud optical depth, one standard deviation, expressed\n" + "as 10 log10(1 + relative error). The extinction errors are assumed\n" + "fully correlated within a profile for each phase, as they are\n" + "dominated by the liquid water path, effective radius and\n" + "reflectivity-temperature relation uncertainties, and independent\n" + "between liquid and ice. Systematic errors, such as liquid in layers\n" + "not classified as droplets or ice detected only by the lidar, are\n" + "not included; see the retrieval status." + ), + "optical_depth": ( + "Vertical integral of extinction over the profile. The value is zero\n" + "in profiles where no cloud was detected. The retrieval is not\n" + "performed, and the value is missing, if rain is present, if liquid\n" + "cloud is present but liquid water path is unavailable or unreliable,\n" + "or if ice is present but affected by uncorrected radar attenuation.\n" + "Pixels where ice was detected only by the lidar do not contribute,\n" + "so the value may be a lower bound; see the retrieval status." + ), +} + +DEFINITIONS = { + "extinction_retrieval_status": utils.status_field_definition( + { + 0: """No cloud detected.""", + 1: """Liquid: reliable retrieval with radar-based droplet + effective radius.""", + 2: """Liquid: droplet effective radius not retrieved from radar, + assumed value used.""", + 3: """Ice: reliable retrieval.""", + 4: """Ice: retrieval performed with radar corrected for liquid, + rain or melting attenuation.""", + 5: """Mixed phase: both liquid and ice extinction retrieved + in the same pixel.""", + 6: """Cloud detected but no retrieval: rain, missing or + non-positive liquid water path, ice detected only by lidar, + or uncorrected radar attenuation.""", + } + ), + "optical_depth_retrieval_status": utils.status_field_definition( + { + 0: """Clear sky: no cloud detected, optical depth is zero.""", + 1: """Reliable retrieval.""", + 2: """Assumed droplet effective radius used in liquid layers + not detected by the radar.""", + 3: """Ice detected only by lidar in part of the profile: ice + optical depth is a lower bound.""", + 4: """Liquid water path exceeds the adiabatic liquid water path of + the detected liquid layers by more than a factor of two: + liquid is probably present in layers not classified as + droplets, and the vertical distribution of liquid + extinction is unreliable.""", + 5: """No retrieval: rain, missing or unreliable liquid water path, + or uncorrected radar attenuation in ice.""", + } + ), +} + +OPTICAL_DEPTH_ATTRIBUTES = { + "extinction_liquid": MetaData( + long_name="Visible extinction coefficient of liquid cloud", + units="m-1", + ancillary_variables="extinction_liquid_error", + comment=COMMENTS["extinction_liquid"], + dimensions=("time", "height"), + ), + "extinction_liquid_error": MetaData( + long_name="Random error in liquid extinction coefficient", + units="dB", + comment=COMMENTS["extinction_liquid_error"], + dimensions=("time", "height"), + ), + "extinction_ice": MetaData( + long_name="Visible extinction coefficient of ice cloud", + units="m-1", + ancillary_variables="extinction_ice_error", + comment=COMMENTS["extinction_ice"], + dimensions=("time", "height"), + ), + "extinction_ice_error": MetaData( + long_name="Random error in ice extinction coefficient", + units="dB", + comment=COMMENTS["extinction_ice_error"], + dimensions=("time", "height"), + ), + "optical_depth_error": MetaData( + long_name="Random error in cloud optical depth", + units="dB", + comment=COMMENTS["optical_depth_error"], + dimensions=("time",), + ), + "extinction_retrieval_status": MetaData( + long_name="Extinction coefficient retrieval status", + definition=DEFINITIONS["extinction_retrieval_status"], + units="1", + dimensions=("time", "height"), + ), + "optical_depth_liquid": MetaData( + long_name="Liquid cloud optical depth", + units="1", + comment=COMMENTS["optical_depth"], + dimensions=("time",), + ), + "optical_depth_ice": MetaData( + long_name="Ice cloud optical depth", + units="1", + comment=COMMENTS["optical_depth"], + dimensions=("time",), + ), + "optical_depth": MetaData( + long_name="Cloud optical depth", + standard_name="atmosphere_optical_thickness_due_to_cloud", + units="1", + ancillary_variables="optical_depth_error optical_depth_retrieval_status", + comment=COMMENTS["optical_depth"], + dimensions=("time",), + ), + "optical_depth_retrieval_status": MetaData( + long_name="Cloud optical depth retrieval status", + definition=DEFINITIONS["optical_depth_retrieval_status"], + units="1", + dimensions=("time",), + ), +} diff --git a/docs/source/api.rst b/docs/source/api.rst index 412ac1c9..975ae893 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -112,6 +112,8 @@ generated. .. autofunction:: products.generate_ier +.. autofunction:: products.generate_optical_depth + .. autofunction:: products.generate_mwr_single .. autofunction:: products.generate_mwr_multi diff --git a/tests/unit/test_optical_depth.py b/tests/unit/test_optical_depth.py new file mode 100644 index 00000000..f09e00fb --- /dev/null +++ b/tests/unit/test_optical_depth.py @@ -0,0 +1,396 @@ +import netCDF4 +import numpy as np +import pytest +from numpy import ma +from numpy.testing import assert_array_almost_equal, assert_array_equal + +from cloudnetpy import constants +from cloudnetpy.products import optical_depth +from cloudnetpy.products.optical_depth import OpticalDepthSource, generate_optical_depth + +# category bits: droplet=1, falling=2, freezing=4 +DROPLET, FALLING, FREEZING = 1, 2, 4 +# quality bits: radar=1, lidar=2 +RADAR, LIDAR = 1, 2 +# adiabatic LWP consistent with the fixture LWP of 0.1 kg m-2 +LWP_ADIABATIC = np.array([0, 0.1, 0, 0.1]) +LWC_REL_ERROR = np.full((4, 4), 0.2) +DER_REL_ERROR = np.full((4, 4), 0.3) +IWC_ERROR_DB = np.full((4, 4), 1.7) + + +@pytest.fixture(scope="session") +def categorize_file(tmpdir_factory): + file_name = tmpdir_factory.mktemp("data").join("categorize.nc") + _create_categorize_file(file_name, with_lwp=True) + return str(file_name) + + +@pytest.fixture(scope="session") +def categorize_file_no_mwr(tmpdir_factory): + file_name = tmpdir_factory.mktemp("data").join("categorize_no_mwr.nc") + _create_categorize_file(file_name, with_lwp=False) + return str(file_name) + + +def _create_categorize_file(file_name, *, with_lwp: bool) -> None: + n_time, n_height = 4, 4 + with netCDF4.Dataset(file_name, "w", format="NETCDF4_CLASSIC") as nc: + for name, n in ( + ("time", n_time), + ("height", n_height), + ("model_time", n_time), + ("model_height", n_height), + ): + nc.createDimension(name, n) + var = nc.createVariable(name, "f8", name) + var[:] = np.arange(n) + nc.variables["height"][:] = [100, 200, 300, 400] + nc.variables["height"].units = "m" + nc.variables["model_height"][:] = [100, 200, 300, 400] + var = nc.createVariable("altitude", "f8") + var[:] = 0 + var.units = "m" + nc.createVariable("radar_frequency", "f8")[:] = 35.5 + if with_lwp: + nc.createVariable("lwp", "f8", "time")[:] = [0.1, 0.1, 0.1, 0.1] + nc.createVariable("lwp_error", "f8", "time")[:] = [0.01, 0.01, 0.01, 0.01] + # profile 0: clear, 1: liquid + ice, 2: ice only, 3: liquid in rain + nc.createVariable("rainfall_rate", "f8", "time")[:] = [0, 0, 0, 1] + cat = np.zeros((n_time, n_height), dtype=int) + cat[1, 0:2] = DROPLET + cat[1, 2:4] = FALLING | FREEZING + cat[2, 1:4] = FALLING | FREEZING + cat[3, 0:2] = DROPLET + nc.createVariable("category_bits", "i4", ("time", "height"))[:] = cat + qual = np.full((n_time, n_height), RADAR | LIDAR, dtype=int) + nc.createVariable("quality_bits", "i4", ("time", "height"))[:] = qual + nc.createVariable("temperature", "f8", ("model_time", "model_height"))[:] = ( + np.full((n_time, n_height), 260.0) + ) + nc.createVariable("pressure", "f8", ("model_time", "model_height"))[:] = ( + np.full((n_time, n_height), 90000.0) + ) + nc.createVariable("Z", "f8", ("time", "height"))[:] = np.full( + (n_time, n_height), -10.0 + ) + nc.createVariable("Z_error", "f8", ("time", "height"))[:] = np.full( + (n_time, n_height), 1.0 + ) + nc.createVariable("is_rain", "i4", "time")[:] = [0, 0, 0, 1] + nc.year, nc.month, nc.day = "2025", "06", "10" + nc.location = "Kumpula" + nc.file_uuid = "b7d3e2f0-1234-5678-9abc-def012345678" + nc.cloudnet_file_type = "categorize" + + +class TestOpticalDepthSource: + @pytest.fixture(autouse=True) + def run_before_tests(self, categorize_file): + self.obj = OpticalDepthSource(categorize_file, assumed_der=10e-6) + yield + self.obj.close() + + def test_liquid_extinction_formula(self, monkeypatch): + lwc = ma.array(np.full((4, 4), 1e-3), mask=~self.obj.is_liquid) + der = ma.array(np.full((4, 4), 20e-6), mask=np.zeros((4, 4), dtype=bool)) + monkeypatch.setattr( + self.obj, "_get_lwc", lambda: (lwc, LWC_REL_ERROR, LWP_ADIABATIC) + ) + monkeypatch.setattr( + self.obj, "_get_der", lambda: (der, ma.array(DER_REL_ERROR)) + ) + self.obj.append_liquid_extinction() + ext = self.obj.data["extinction_liquid"][:] + expected = 3 * 1e-3 / (2 * 1000 * 20e-6) + assert_array_almost_equal(ext[1, 0:2], expected) + assert ext[0, :].mask.all() + assert ext[1, 2:4].mask.all() + assert not self.obj.is_der_assumed.any() + + def test_liquid_extinction_assumed_der(self, monkeypatch): + lwc = ma.array(np.full((4, 4), 1e-3), mask=~self.obj.is_liquid) + der = ma.masked_all((4, 4)) + der[1, 0] = 20e-6 + der[1, 1] = 100e-6 # outside the valid range + monkeypatch.setattr( + self.obj, "_get_lwc", lambda: (lwc, LWC_REL_ERROR, LWP_ADIABATIC) + ) + monkeypatch.setattr( + self.obj, "_get_der", lambda: (der, ma.array(DER_REL_ERROR)) + ) + self.obj.append_liquid_extinction() + ext = self.obj.data["extinction_liquid"][:] + assert_array_almost_equal(ext[1, 0], 3 * 1e-3 / (2 * 1000 * 20e-6)) + assert_array_almost_equal(ext[1, 1], 3 * 1e-3 / (2 * 1000 * 10e-6)) + assert_array_equal(self.obj.is_der_assumed[1, :], [False, True, False, False]) + assert_array_equal(self.obj.is_der_assumed[3, :], [True, True, False, False]) + + def test_extended_cloud_top_pixel_is_kept(self, monkeypatch): + # lwc retrieval placed liquid in pixel (1, 2), which has no droplet bit + lwc = ma.array(np.full((4, 4), 1e-3), mask=~self.obj.is_liquid) + lwc[1, 2] = 1e-3 + der = ma.masked_all((4, 4)) + monkeypatch.setattr( + self.obj, "_get_lwc", lambda: (lwc, LWC_REL_ERROR, LWP_ADIABATIC) + ) + monkeypatch.setattr( + self.obj, "_get_der", lambda: (der, ma.array(DER_REL_ERROR)) + ) + self.obj.append_liquid_extinction() + ext = self.obj.data["extinction_liquid"][:] + assert ext[1, 2] > 0 + assert ext[1, 3] is ma.masked + assert (DER_REL_ERROR == 0.3).all() # input constant not mutated + + def test_ice_extinction_formula(self, monkeypatch): + is_ice = self.obj.ice_classification.is_ice + iwc = ma.array(np.full((4, 4), 1e-4), mask=~is_ice) + ier = ma.array(np.full((4, 4), 50e-6), mask=~is_ice) + iwc[2, 1] = ma.masked # lidar-only ice pixel + monkeypatch.setattr(self.obj, "_get_iwc", lambda: (iwc, ma.array(IWC_ERROR_DB))) + monkeypatch.setattr(self.obj, "_get_ier", lambda: ier) + self.obj.append_ice_extinction() + ext = self.obj.data["extinction_ice"][:] + expected = 3 * 1e-4 / (2 * constants.RHO_ICE * 50e-6) + assert_array_almost_equal(ext[1, 2:4], expected) + assert_array_almost_equal(ext[2, 2:4], expected) + assert ext[2, 1] is ma.masked + assert ext[0, :].mask.all() + assert_array_equal( + self.obj.is_lidar_only_ice[2, :], [False, True, False, False] + ) + + def test_ice_masked_in_rain(self, monkeypatch): + is_ice = self.obj.ice_classification.is_ice + iwc = ma.array(np.full((4, 4), 1e-4), mask=~is_ice) + ier = ma.array(np.full((4, 4), 50e-6), mask=~is_ice) + self.obj.is_rain = np.array([False, False, True, False]) + monkeypatch.setattr(self.obj, "_get_iwc", lambda: (iwc, ma.array(IWC_ERROR_DB))) + monkeypatch.setattr(self.obj, "_get_ier", lambda: ier) + self.obj.append_ice_extinction() + lwc = ma.masked_all((4, 4)) + monkeypatch.setattr( + self.obj, "_get_lwc", lambda: (lwc, LWC_REL_ERROR, LWP_ADIABATIC) + ) + monkeypatch.setattr(self.obj, "_get_der", lambda: (lwc, lwc)) + self.obj.append_liquid_extinction() + self.obj.append_extinction_status() + assert self.obj.data["extinction_ice"][:][2, :].mask.all() + assert_array_equal(self.obj.data["extinction_retrieval_status"][:][2, 1:], 6) + + def test_optical_depths_and_status(self, monkeypatch): + is_ice = self.obj.ice_classification.is_ice + lwc = ma.array(np.full((4, 4), 1e-3), mask=~self.obj.is_liquid) + lwc[3, :] = ma.masked # rain profile has no lwc + der = ma.array(np.full((4, 4), 20e-6), mask=~self.obj.is_liquid) + der[1, 1] = ma.masked + iwc = ma.array(np.full((4, 4), 1e-4), mask=~is_ice) + iwc[2, 1] = ma.masked + ier = ma.array(np.full((4, 4), 50e-6), mask=~is_ice) + monkeypatch.setattr( + self.obj, "_get_lwc", lambda: (lwc, LWC_REL_ERROR, LWP_ADIABATIC) + ) + monkeypatch.setattr( + self.obj, "_get_der", lambda: (der, ma.array(DER_REL_ERROR)) + ) + monkeypatch.setattr(self.obj, "_get_iwc", lambda: (iwc, ma.array(IWC_ERROR_DB))) + monkeypatch.setattr(self.obj, "_get_ier", lambda: ier) + self.obj.append_liquid_extinction() + self.obj.append_ice_extinction() + self.obj.append_extinction_status() + self.obj.append_optical_depths() + self.obj.append_optical_depth_error() + self.obj.append_optical_depth_status() + + ext_liq = 3 * 1e-3 / (2 * 1000 * 20e-6) + ext_liq_assumed = 3 * 1e-3 / (2 * 1000 * 10e-6) + ext_ice = 3 * 1e-4 / (2 * constants.RHO_ICE * 50e-6) + dz = 100 + + tau_liq = self.obj.data["optical_depth_liquid"][:] + tau_ice = self.obj.data["optical_depth_ice"][:] + tau = self.obj.data["optical_depth"][:] + assert_array_almost_equal( + tau_liq[0:3], [0, (ext_liq + ext_liq_assumed) * dz, 0] + ) + assert_array_almost_equal(tau_ice[0:3], [0, 2 * ext_ice * dz, 2 * ext_ice * dz]) + assert_array_almost_equal(tau[0:3], tau_liq[0:3] + tau_ice[0:3]) + assert tau_liq[3] is ma.masked + assert tau_ice[3] is ma.masked + assert tau[3] is ma.masked + + pixel_status = self.obj.data["extinction_retrieval_status"][:] + assert_array_equal(pixel_status[0, :], [0, 0, 0, 0]) + assert_array_equal(pixel_status[1, :], [1, 2, 3, 3]) + assert_array_equal(pixel_status[2, :], [0, 6, 3, 3]) + assert_array_equal(pixel_status[3, :], [6, 6, 0, 0]) + + status = self.obj.data["optical_depth_retrieval_status"][:] + assert_array_equal(status, [0, 2, 3, 5]) + + # errors: 10 log10(1 + relative error) + err_liq = self.obj.data["extinction_liquid_error"][:] + assert_array_almost_equal(err_liq[1, 0], 10 * np.log10(1 + np.hypot(0.2, 0.3))) + assert_array_almost_equal(err_liq[1, 1], 10 * np.log10(1 + np.hypot(0.2, 0.4))) + assert err_liq[0, :].mask.all() + err_ice = self.obj.data["extinction_ice_error"][:] + assert_array_almost_equal(err_ice[2, 2:4], 1.7) + assert err_ice[2, 1] is ma.masked + err = self.obj.data["optical_depth_error"][:] + rel_ice = 10 ** (1.7 / 10) - 1 + rel_liq = np.array([np.hypot(0.2, 0.3), np.hypot(0.2, 0.4)]) + w = np.array([ext_liq, ext_liq_assumed]) + rel_liq_col = np.sum(rel_liq * w) / np.sum(w) + abs_err = np.hypot(rel_liq_col * tau_liq[1], rel_ice * tau_ice[1]) + assert_array_almost_equal(err[1], 10 * np.log10(1 + abs_err / tau[1])) + assert_array_almost_equal(err[2], 1.7) + assert err[0] is ma.masked + assert err[3] is ma.masked + + def test_corrected_ice_status(self, monkeypatch): + is_ice = self.obj.ice_classification.is_ice + lwc = ma.masked_all((4, 4)) + der = ma.masked_all((4, 4)) + iwc = ma.array(np.full((4, 4), 1e-4), mask=~is_ice) + ier = ma.array(np.full((4, 4), 50e-6), mask=~is_ice) + corrected = np.zeros((4, 4), dtype=bool) + corrected[2, 2] = True + monkeypatch.setattr(self.obj.ice_classification, "corrected_ice", corrected) + monkeypatch.setattr( + self.obj, "_get_lwc", lambda: (lwc, LWC_REL_ERROR, LWP_ADIABATIC) + ) + monkeypatch.setattr( + self.obj, "_get_der", lambda: (der, ma.array(DER_REL_ERROR)) + ) + monkeypatch.setattr(self.obj, "_get_iwc", lambda: (iwc, ma.array(IWC_ERROR_DB))) + monkeypatch.setattr(self.obj, "_get_ier", lambda: ier) + self.obj.append_liquid_extinction() + self.obj.append_ice_extinction() + self.obj.append_extinction_status() + self.obj.append_optical_depths() + self.obj.append_optical_depth_status() + assert self.obj.data["extinction_retrieval_status"][:][2, 2] == 4 + assert self.obj.data["optical_depth_retrieval_status"][:][2] == 1 + # liquid present without lwc -> no retrieval + assert self.obj.data["optical_depth_retrieval_status"][:][1] == 5 + assert self.obj.data["optical_depth"][:][1] is ma.masked + + def test_inconsistent_lwp_status(self, monkeypatch): + is_ice = self.obj.ice_classification.is_ice + lwc = ma.array(np.full((4, 4), 1e-3), mask=~self.obj.is_liquid) + der = ma.array(np.full((4, 4), 20e-6), mask=~self.obj.is_liquid) + iwc = ma.array(np.full((4, 4), 1e-4), mask=~is_ice) + ier = ma.array(np.full((4, 4), 50e-6), mask=~is_ice) + # profile 1 can hold only 0.02 kg m-2 but LWP is 0.1 kg m-2 + lwp_adiabatic = np.array([0, 0.02, 0, 0.1]) + monkeypatch.setattr( + self.obj, "_get_lwc", lambda: (lwc, LWC_REL_ERROR, lwp_adiabatic) + ) + monkeypatch.setattr( + self.obj, "_get_der", lambda: (der, ma.array(DER_REL_ERROR)) + ) + monkeypatch.setattr(self.obj, "_get_iwc", lambda: (iwc, ma.array(IWC_ERROR_DB))) + monkeypatch.setattr(self.obj, "_get_ier", lambda: ier) + self.obj.append_liquid_extinction() + self.obj.append_ice_extinction() + self.obj.append_extinction_status() + self.obj.append_optical_depths() + self.obj.append_optical_depth_status() + assert_array_equal(self.obj.is_lwp_inconsistent, [False, True, False, False]) + status = self.obj.data["optical_depth_retrieval_status"][:] + assert_array_equal(status, [0, 4, 1, 5]) + # value is kept, only flagged + assert self.obj.data["optical_depth"][:][1] > 0 + + +def test_generate_optical_depth(categorize_file, tmp_path): + output_file = tmp_path / "optical_depth.nc" + uuid = generate_optical_depth(categorize_file, output_file) + with netCDF4.Dataset(output_file) as nc: + assert nc.file_uuid == str(uuid) + assert nc.cloudnet_file_type == "optical-depth" + for key in ( + "extinction_liquid", + "extinction_ice", + "extinction_retrieval_status", + "optical_depth_liquid", + "optical_depth_ice", + "optical_depth", + "optical_depth_retrieval_status", + "extinction_liquid_error", + "extinction_ice_error", + "optical_depth_error", + "lwp", + ): + assert key in nc.variables + assert nc.variables["extinction_liquid"].dimensions == ("time", "height") + assert nc.variables["optical_depth"].dimensions == ("time",) + assert nc.variables["optical_depth"].units == "1" + assert "10 um" in nc.variables["extinction_liquid"].comment + tau = nc.variables["optical_depth"][:] + status = nc.variables["optical_depth_retrieval_status"][:] + assert tau[0] == 0 + assert status[0] == 0 + assert tau[1] > 0 + assert tau[2] > 0 + assert tau[3] is ma.masked + assert status[3] == 5 + # thin fixture layer cannot hold 0.1 kg m-2 adiabatically + assert status[1] == 4 + ext_liq = nc.variables["extinction_liquid"][:] + assert ext_liq[1, 0:2].count() == 2 + ext_ice = nc.variables["extinction_ice"][:] + assert ext_ice[1, 2:4].count() == 2 + assert ext_ice[2, 1:4].count() == 3 + + +def test_non_positive_lwp_gives_no_liquid_retrieval(categorize_file, tmp_path): + src = OpticalDepthSource(categorize_file, assumed_der=10e-6) + src.dataset.variables["lwp"][:] # LWP is 0.1 in the fixture + lwc, _, _ = src._get_lwc() + assert not lwc[1, 0:2].mask.any() + src.close() + with netCDF4.Dataset(categorize_file, "a") as nc: + nc.variables["lwp"][1] = 0.0 + try: + output_file = tmp_path / "optical_depth.nc" + generate_optical_depth(categorize_file, output_file) + with netCDF4.Dataset(output_file) as nc: + assert nc.variables["optical_depth"][:][1] is ma.masked + assert nc.variables["optical_depth_retrieval_status"][:][1] == 5 + assert_array_equal( + nc.variables["extinction_retrieval_status"][:][1, 0:2], 6 + ) + finally: + with netCDF4.Dataset(categorize_file, "a") as nc: + nc.variables["lwp"][1] = 0.1 + + +def test_generate_without_mwr(categorize_file_no_mwr, tmp_path): + output_file = tmp_path / "optical_depth.nc" + generate_optical_depth(categorize_file_no_mwr, output_file) + with netCDF4.Dataset(output_file) as nc: + assert "lwp" not in nc.variables + assert nc.variables["extinction_liquid"][:].mask.all() + tau = nc.variables["optical_depth"][:] + status = nc.variables["optical_depth_retrieval_status"][:] + assert tau[2] > 0 # ice-only profile still retrieved + assert status[2] == 1 + assert tau[1] is ma.masked # liquid present, no LWP + assert status[1] == 5 + + +def test_generate_optical_depth_custom_der(categorize_file, tmp_path): + output_file = tmp_path / "optical_depth.nc" + generate_optical_depth(categorize_file, output_file, assumed_der=5e-6) + with netCDF4.Dataset(output_file) as nc: + assert "5 um" in nc.variables["extinction_liquid"].comment + + +def test_attributes_cover_all_variables(categorize_file, tmp_path): + output_file = tmp_path / "optical_depth.nc" + generate_optical_depth(categorize_file, output_file) + with netCDF4.Dataset(output_file) as nc: + for key in optical_depth.OPTICAL_DEPTH_ATTRIBUTES: + assert nc.variables[key].long_name From f91939bf66f7622b9b4e0a919bfeb641c78fa116 Mon Sep 17 00:00:00 2001 From: Simo Tukiainen Date: Thu, 10 Sep 2026 12:21:42 +0300 Subject: [PATCH 2/4] Address review comments --- cloudnetpy/output.py | 4 ++-- cloudnetpy/products/optical_depth.py | 13 ++++++++++--- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/cloudnetpy/output.py b/cloudnetpy/output.py index f5d63133..f2fd8504 100644 --- a/cloudnetpy/output.py +++ b/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 += ( @@ -198,7 +198,7 @@ def get_references(identifier: str | None = None, extra: list | None = None) -> case "optical-depth": references += ( ", https://doi.org/10.1175/BAMS-88-6-883" - ", 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" ", https://doi.org/10.1175/JAM2340.1" ", https://doi.org/10.1175/JAM2543.1" ", https://doi.org/10.5194/amt-13-5335-2020" diff --git a/cloudnetpy/products/optical_depth.py b/cloudnetpy/products/optical_depth.py index d1079968..8280bece 100644 --- a/cloudnetpy/products/optical_depth.py +++ b/cloudnetpy/products/optical_depth.py @@ -104,7 +104,11 @@ def generate_optical_depth( attributes = _add_extinction_comments(attributes, od_source) output.update_attributes(od_source.data, attributes) output.save_product_file( - "optical-depth", od_source, output_file, uuid, copy_from_cat=("lwp",) + "optical-depth", + od_source, + output_file, + uuid, + copy_from_cat=("lwp", "lwp_error"), ) return uuid @@ -113,6 +117,9 @@ class OpticalDepthSource(DataSource): """Data container for cloud optical depth calculations.""" def __init__(self, categorize_file: str | PathLike, assumed_der: float) -> None: + if not np.isfinite(assumed_der) or assumed_der <= 0: + msg = "Assumed droplet effective radius must be finite and positive." + raise ValueError(msg) super().__init__(categorize_file) self.categorize_file = categorize_file self.assumed_der = assumed_der @@ -401,8 +408,8 @@ def _add_extinction_comments(attributes: dict, od_source: OpticalDepthSource) -> { 0: """Clear sky: no cloud detected, optical depth is zero.""", 1: """Reliable retrieval.""", - 2: """Assumed droplet effective radius used in liquid layers - not detected by the radar.""", + 2: """Assumed droplet effective radius used where the radar-based + retrieval was unavailable or outside the valid range.""", 3: """Ice detected only by lidar in part of the profile: ice optical depth is a lower bound.""", 4: """Liquid water path exceeds the adiabatic liquid water path of From 380edb676489af1c9b9c30ee4b8c011ceb2b56a2 Mon Sep 17 00:00:00 2001 From: Simo Tukiainen Date: Thu, 10 Sep 2026 14:03:16 +0300 Subject: [PATCH 3/4] Rename product id to cod --- cloudnetpy/cli.py | 4 +--- cloudnetpy/output.py | 6 ++--- cloudnetpy/products/__init__.py | 2 +- .../products/{optical_depth.py => cod.py} | 8 +++---- docs/source/api.rst | 2 +- .../{test_optical_depth.py => test_cod.py} | 22 +++++++++---------- 6 files changed, 21 insertions(+), 23 deletions(-) rename cloudnetpy/products/{optical_depth.py => cod.py} (99%) rename tests/unit/{test_optical_depth.py => test_cod.py} (96%) diff --git a/cloudnetpy/cli.py b/cloudnetpy/cli.py index 0311c986..429ef855 100644 --- a/cloudnetpy/cli.py +++ b/cloudnetpy/cli.py @@ -822,9 +822,7 @@ def _plot_l3( def _process_cat_product(product: str, categorize_file: str) -> str: output_file = categorize_file.replace("categorize", product) module = importlib.import_module("cloudnetpy.products") - getattr(module, f"generate_{product.replace('-', '_')}")( - categorize_file, output_file - ) + getattr(module, f"generate_{product}")(categorize_file, output_file) logging.info("Processed %s: %s", product, output_file) return output_file diff --git a/cloudnetpy/output.py b/cloudnetpy/output.py index f2fd8504..084e42c7 100644 --- a/cloudnetpy/output.py +++ b/cloudnetpy/output.py @@ -195,7 +195,7 @@ 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 "optical-depth": + 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" @@ -480,7 +480,7 @@ def _get_identifier(short_id: str) -> str: "classification", "der", "ier", - "optical-depth", + "cod", "classification-voodoo", "epsilon-radar", ) @@ -495,7 +495,7 @@ def _get_identifier(short_id: str) -> str: return "ice effective radius" if short_id == "der": return "droplet effective radius" - if short_id == "optical-depth": + if short_id == "cod": return "cloud optical depth" if short_id == "epsilon-radar": return "dissipation rate of turbulent kinetic energy" diff --git a/cloudnetpy/products/__init__.py b/cloudnetpy/products/__init__.py index c89e9121..b88a8617 100644 --- a/cloudnetpy/products/__init__.py +++ b/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 @@ -7,4 +8,3 @@ from .iwc import generate_iwc from .lwc import generate_lwc from .mwr_tools import generate_mwr_lhumpro, generate_mwr_multi, generate_mwr_single -from .optical_depth import generate_optical_depth diff --git a/cloudnetpy/products/optical_depth.py b/cloudnetpy/products/cod.py similarity index 99% rename from cloudnetpy/products/optical_depth.py rename to cloudnetpy/products/cod.py index 8280bece..71107411 100644 --- a/cloudnetpy/products/optical_depth.py +++ b/cloudnetpy/products/cod.py @@ -33,7 +33,7 @@ ASSUMED_DER_REL_ERROR = 0.4 -def generate_optical_depth( +def generate_cod( categorize_file: str | PathLike, output_file: str | PathLike, uuid: str | UUID | None = None, @@ -63,8 +63,8 @@ def generate_optical_depth( UUID of the generated file. Examples: - >>> from cloudnetpy.products import generate_optical_depth - >>> generate_optical_depth('categorize.nc', 'optical_depth.nc') + >>> from cloudnetpy.products import generate_cod + >>> generate_cod('categorize.nc', 'cod.nc') References: Frisch, S., Shupe, M., Djalalova, I., Feingold, G., & Poellot, M. (2002). @@ -104,7 +104,7 @@ def generate_optical_depth( attributes = _add_extinction_comments(attributes, od_source) output.update_attributes(od_source.data, attributes) output.save_product_file( - "optical-depth", + "cod", od_source, output_file, uuid, diff --git a/docs/source/api.rst b/docs/source/api.rst index 975ae893..8cb230d2 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -112,7 +112,7 @@ generated. .. autofunction:: products.generate_ier -.. autofunction:: products.generate_optical_depth +.. autofunction:: products.generate_cod .. autofunction:: products.generate_mwr_single diff --git a/tests/unit/test_optical_depth.py b/tests/unit/test_cod.py similarity index 96% rename from tests/unit/test_optical_depth.py rename to tests/unit/test_cod.py index f09e00fb..4b52c6e3 100644 --- a/tests/unit/test_optical_depth.py +++ b/tests/unit/test_cod.py @@ -5,8 +5,8 @@ from numpy.testing import assert_array_almost_equal, assert_array_equal from cloudnetpy import constants -from cloudnetpy.products import optical_depth -from cloudnetpy.products.optical_depth import OpticalDepthSource, generate_optical_depth +from cloudnetpy.products import cod +from cloudnetpy.products.cod import OpticalDepthSource, generate_cod # category bits: droplet=1, falling=2, freezing=4 DROPLET, FALLING, FREEZING = 1, 2, 4 @@ -304,12 +304,12 @@ def test_inconsistent_lwp_status(self, monkeypatch): assert self.obj.data["optical_depth"][:][1] > 0 -def test_generate_optical_depth(categorize_file, tmp_path): +def test_generate_cod(categorize_file, tmp_path): output_file = tmp_path / "optical_depth.nc" - uuid = generate_optical_depth(categorize_file, output_file) + uuid = generate_cod(categorize_file, output_file) with netCDF4.Dataset(output_file) as nc: assert nc.file_uuid == str(uuid) - assert nc.cloudnet_file_type == "optical-depth" + assert nc.cloudnet_file_type == "cod" for key in ( "extinction_liquid", "extinction_ice", @@ -355,7 +355,7 @@ def test_non_positive_lwp_gives_no_liquid_retrieval(categorize_file, tmp_path): nc.variables["lwp"][1] = 0.0 try: output_file = tmp_path / "optical_depth.nc" - generate_optical_depth(categorize_file, output_file) + generate_cod(categorize_file, output_file) with netCDF4.Dataset(output_file) as nc: assert nc.variables["optical_depth"][:][1] is ma.masked assert nc.variables["optical_depth_retrieval_status"][:][1] == 5 @@ -369,7 +369,7 @@ def test_non_positive_lwp_gives_no_liquid_retrieval(categorize_file, tmp_path): def test_generate_without_mwr(categorize_file_no_mwr, tmp_path): output_file = tmp_path / "optical_depth.nc" - generate_optical_depth(categorize_file_no_mwr, output_file) + generate_cod(categorize_file_no_mwr, output_file) with netCDF4.Dataset(output_file) as nc: assert "lwp" not in nc.variables assert nc.variables["extinction_liquid"][:].mask.all() @@ -381,16 +381,16 @@ def test_generate_without_mwr(categorize_file_no_mwr, tmp_path): assert status[1] == 5 -def test_generate_optical_depth_custom_der(categorize_file, tmp_path): +def test_generate_cod_custom_der(categorize_file, tmp_path): output_file = tmp_path / "optical_depth.nc" - generate_optical_depth(categorize_file, output_file, assumed_der=5e-6) + generate_cod(categorize_file, output_file, assumed_der=5e-6) with netCDF4.Dataset(output_file) as nc: assert "5 um" in nc.variables["extinction_liquid"].comment def test_attributes_cover_all_variables(categorize_file, tmp_path): output_file = tmp_path / "optical_depth.nc" - generate_optical_depth(categorize_file, output_file) + generate_cod(categorize_file, output_file) with netCDF4.Dataset(output_file) as nc: - for key in optical_depth.OPTICAL_DEPTH_ATTRIBUTES: + for key in cod.OPTICAL_DEPTH_ATTRIBUTES: assert nc.variables[key].long_name From 323488ffd804f8764aa2132fa2e96fa49df2371f Mon Sep 17 00:00:00 2001 From: Simo Tukiainen Date: Thu, 10 Sep 2026 14:21:32 +0300 Subject: [PATCH 4/4] Mark clear sky in log-scale plots --- cloudnetpy/plotting/plotting.py | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/cloudnetpy/plotting/plotting.py b/cloudnetpy/plotting/plotting.py index 217b615a..b0390b8f 100644 --- a/cloudnetpy/plotting/plotting.py +++ b/cloudnetpy/plotting/plotting.py @@ -904,6 +904,7 @@ 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) @@ -921,11 +922,29 @@ def plot(self, figure_data: FigureData, hacky_freq_ind: int | None = None) -> No self._fill_between_data_gaps(figure_data) if self._is_log: self._ax.set_yscale("log") - self.sub_plot.set_yax(ylabel=units, y_limits=self._get_y_limits()) + 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)