Implement epsilon-radar product

cloudnetpy/categorize/model.py

@@ -11,7 +11,6 @@
 
 from cloudnetpy import utils
 from cloudnetpy.cloudnetarray import CloudnetArray
-from cloudnetpy.constants import G
 from cloudnetpy.datasource import DataSource
 from cloudnetpy.exceptions import ModelDataError
 
@@ -146,27 +145,9 @@ def _get_model_heights(self, alt_site: float) -> npt.NDArray:
         except KeyError as err:
             msg = "No 'height' variable in the model file."
             raise ModelDataError(msg) from err
-        surface_altitude = self._get_model_surface_altitude(alt_site)
+        surface_altitude = utils.get_model_surface_altitude(self.dataset, alt_site)
         return self.to_m(model_heights) + surface_altitude
 
-    def _get_model_surface_altitude(self, alt_site: float) -> float | npt.NDArray:
-        """Returns model surface altitude from geopotential if available.
-
-        For sites in complex terrain (e.g. mountains), the model grid cell
-        surface height can differ significantly from the actual site altitude.
-        Using the model's own surface height ensures that thermodynamic fields
-        are placed at their physically correct absolute heights.
-
-        Note: Model surface altitude might be higher than the site altitude.
-        """
-        if "sfc_height" in self.dataset.variables:
-            sfc_height = self.dataset.variables["sfc_height"][:]
-            return sfc_height[:, np.newaxis]
-        if "sfc_geopotential" in self.dataset.variables:
-            geopotential = self.dataset.variables["sfc_geopotential"][:]
-            return geopotential[:, np.newaxis] / G
-        return alt_site
-
     def calc_attenuations(self, frequency: float) -> None:
         temperature = self.getvar("temperature")
         pressure = self.getvar("pressure")

cloudnetpy/cli.py

@@ -90,6 +90,33 @@ def run(args: argparse.Namespace, tmpdir: str, client: APIClient) -> None:
         l2_filename = _process_mwrpy_product(product, mwrpy_filepath, args)
         _plot(l2_filename, product, args)
 
+    # Radar + model based products (e.g. epsilon-radar)
+    if "epsilon-radar" in args.products:
+        epsilon_filepath = _process_epsilon_radar(cat_files, args, client)
+        if epsilon_filepath is not None:
+            _plot(epsilon_filepath, "epsilon-radar", args)
+
+
+def _process_epsilon_radar(
+    cat_files: dict, args: argparse.Namespace, client: APIClient
+) -> str | None:
+    radar_filepath = cat_files.get("radar") or _fetch_product(args, "radar", client)
+    if radar_filepath is None:
+        logging.info("No radar data available for epsilon-radar")
+        return None
+    model_filepath = _fetch_model(args, client)
+    if model_filepath is None:
+        logging.info("No model data available for epsilon-radar")
+        return None
+    filename = f"{args.date.replace('-', '')}_{args.site}_epsilon-radar.nc"
+    output_file = _create_output_folder("geophysical", args) / filename
+    products = importlib.import_module("cloudnetpy.products")
+    products.generate_epsilon_from_radar(
+        str(radar_filepath), model_filepath, str(output_file)
+    )
+    logging.info("Processed epsilon-radar: %s", output_file)
+    return str(output_file)
+
 
 def _process_categorize(
     input_files: dict,

cloudnetpy/output.py

@@ -74,6 +74,7 @@ def save_product_file(
     file_name: str | PathLike,
     uuid: UUID,
     copy_from_cat: tuple = (),
+    extra_sources: tuple[DataSource, ...] = (),
 ) -> None:
     """Saves a standard Cloudnet product file.
 
@@ -83,13 +84,18 @@ def save_product_file(
         file_name: Name of the output file to be generated.
         uuid: Set specific UUID for the file.
         copy_from_cat: Variables to be copied from the categorize file.
+        extra_sources: Additional input DataSources whose ``file_uuid`` should
+            also be listed in ``source_file_uuids`` (e.g. a model file used
+            alongside the primary L1b input).
 
     """
     human_readable_file_type = _get_identifier(short_id)
-    dimensions = {
-        "time": len(obj.time),
-        "height": len(obj.dataset.variables["height"]),
-    }
+    height_size = (
+        len(obj.data["height"][:])
+        if "height" in obj.data
+        else len(obj.dataset.variables["height"])
+    )
+    dimensions = {"time": len(obj.time), "height": height_size}
     with init_file(file_name, dimensions, obj.data, uuid) as nc:
         nc.cloudnet_file_type = short_id
         vars_from_source = (
@@ -105,7 +111,7 @@ def save_product_file(
             f"{human_readable_file_type.capitalize()} products from"
             f" {obj.dataset.location}"
         )
-        nc.source_file_uuids = get_source_uuids([nc, obj])
+        nc.source_file_uuids = get_source_uuids([nc, obj, *extra_sources])
         copy_global(
             obj.dataset,
             nc,
@@ -116,13 +122,28 @@ def save_product_file(
                 "year",
                 "source",
                 "source_instrument_pids",
+                "instrument_pid",
                 "voodoonet_version",
             ),
         )
-        merge_history(nc, human_readable_file_type, obj)
+        _append_extra_sources(nc, extra_sources)
+        merge_history(nc, human_readable_file_type, obj, extra_sources=extra_sources)
         nc.references = get_references(short_id)
 
 
+def _append_extra_sources(
+    nc: netCDF4.Dataset, extra_sources: tuple[DataSource, ...]
+) -> None:
+    """Merges ``source`` strings from auxiliary input files into ``nc.source``."""
+    if not extra_sources:
+        return
+    existing = nc.source.split("\n") if "source" in nc.ncattrs() else []
+    extras = [src.dataset.source for src in extra_sources if src.dataset.source]
+    merged = list(dict.fromkeys([*existing, *extras]))
+    if merged:
+        nc.source = "\n".join(merged)
+
+
 def get_l1b_source(instrument: Instrument) -> str:
     """Returns level 1b file source."""
     parts = [
@@ -173,6 +194,11 @@ 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 "epsilon-radar":
+            references += (
+                ", https://doi.org/10.5194/amt-13-5335-2020"
+                ", https://doi.org/10.1002/2015JD024543"
+            )
     if extra is not None:
         for reference in extra:
             references += f", {reference}"
@@ -203,14 +229,21 @@ def get_source_uuids(data: Observations | list[netCDF4.Dataset | DataSource]) ->
 
 
 def merge_history(
-    nc: netCDF4.Dataset, file_type: str, data: Observations | DataSource
+    nc: netCDF4.Dataset,
+    file_type: str,
+    data: Observations | DataSource,
+    extra_sources: tuple[DataSource, ...] = (),
 ) -> None:
     """Merges history fields from one or several files and creates a new record."""
 
     def extract_history(obj: DataSource | Observations) -> list[str]:
         if hasattr(obj, "dataset") and hasattr(obj.dataset, "history"):
             history = obj.dataset.history
-            if isinstance(obj, Model):
+            is_model_file = (
+                isinstance(obj, Model)
+                or getattr(obj.dataset, "cloudnet_file_type", "") == "model"
+            )
+            if is_model_file:
                 return [history.split("\n")[-1]]
             return history.split("\n")
         return []
@@ -221,6 +254,8 @@ def extract_history(obj: DataSource | Observations) -> list[str]:
     elif isinstance(data, Observations):
         for field in fields(data):
             histories.extend(extract_history(getattr(data, field.name)))
+    for src in extra_sources:
+        histories.extend(extract_history(src))
 
     # Remove duplicates
     histories = list(dict.fromkeys(histories))
@@ -298,7 +333,7 @@ def copy_variables(
 
     """
     for key in keys:
-        if key in source.variables:
+        if key in source.variables and key not in target.variables:
             fill_value = getattr(source.variables[key], "_FillValue", False)
             variable = source.variables[key]
             var_out = target.createVariable(
@@ -436,6 +471,7 @@ def _get_identifier(short_id: str) -> str:
         "der",
         "ier",
         "classification-voodoo",
+        "epsilon-radar",
     )
     if short_id not in valid_ids:
         msg = f"Invalid file identifier: {short_id}"
@@ -448,6 +484,8 @@ def _get_identifier(short_id: str) -> str:
         return "ice effective radius"
     if short_id == "der":
         return "droplet effective radius"
+    if short_id == "epsilon-radar":
+        return "dissipation rate of turbulent kinetic energy"
     return short_id
 
 

cloudnetpy/plotting/plot_meta.py

@@ -667,5 +667,10 @@ class PlotMeta(NamedTuple):
             plot_range=(1e-7, 1e-1),
             log_scale=True,
         ),
+        "epsilon_error": PlotMeta(
+            cmap="inferno",
+            plot_range=(1e-7, 1e-1),
+            log_scale=True,
+        ),
     },
 }

cloudnetpy/products/__init__.py

@@ -1,6 +1,7 @@
 from .classification import generate_classification
 from .der import generate_der
 from .drizzle import generate_drizzle
+from .epsilon_radar import generate_epsilon_from_radar
 from .ier import generate_ier
 from .iwc import generate_iwc
 from .lwc import generate_lwc