Amaroc 841 localized base height calculation (#137)

* AMAROC-835 add ceilo id query for base height calculation

* AMAROC-835 request minimum amount of hits for base height calculation

* AMAROC-841 localized cloud base height calculation

* AMAROC-842 resolve version conflict

* AMAROC-841 AMAROC-842 adapt changelog

* AMAROC-841 raise warning if we do not filter but filter is set + adapt test

* AMAROC-841 add documentation

* AMAROC-841 fix doc referencing

CHANGELOG

@@ -3,6 +3,12 @@ All notable changes to ampycloud will be documented in this file.
 The format is inspired from [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [v2.1.1]
+### Added:
+ - [regDaniel, 2024-10-28] Feature to calculate cloud base height with a subset of instruments only
+### Fixed:
+ - [regDaniel, 2024-10-28] Incompatibility of scipy and statsmodels
+
 ## [v2.0.0]
 ### Added:
  - [fpavogt, 2024-06-26] Add ability to trigger weekly tests manually.

docs/source/changes.rst

@@ -0,0 +1,29 @@
+.. include:: ./substitutions.rst
+
+.. _changes:
+
+Scientific changes since v2.0.0
+===============================
+
+The scientific documentation of Ampycloud v2.0.0 can be found
+`here <https://amt.copernicus.org/articles/17/4891/2024/>`_. In this page, we
+list all changes since v.2.0.0 that go beyond bugfixing, refactoring, patching,
+etc. and have an impact on the science of the algorithm. More detailed
+information on changes can be found in the `changelog <changelog>`.
+
+
+v2.1.0: Enable ceilometer filtering for calculation of cloud base height
+------------------------------------------------------------------------
+
+There might be situations, where it is beneficial to calculate the cloud base
+height from a subset of ceilometer hits reported by specific ceilometers. For
+example:
+- If the cloud height is supposed to be representative for a given location, but
+you still want to use as many ceilometers as possible to infer the amount.
+- If you use different ceilometer models and know that you want to calculate the
+height only from hits of a specific ceilometer model to avoid implementing
+complicated correction factors.
+To this end, the parameter ``CEILOS_FOR_BASE_HEIGHT_CALC`` was implemented in
+this version. The default value is an empty list. In order to activate the
+filtering, it is sufficient to enter the ceilometer IDs of the ceilos to keep
+for the base height calculation.

docs/source/index.rst

@@ -40,7 +40,8 @@ Welcome to the ampycloud documentation
 * **Who:** ampycloud is being developed at `MeteoSwiss <https://www.meteoswiss.admin.ch>`__.
   See also the code's :ref:`license & copyright <license:License & Copyright>` information.
 
-* **How:** a scientific article describing the ampycloud **algorithm** is currently in preparation.
+* **How:** a scientific article describing the ampycloud **algorithm** is published
+  `here <https://amt.copernicus.org/articles/17/4891/2024/>`_.
   This article will be complemented by these webpages, that contain the official documentation of
   the ampycloud **Python package**.
 
@@ -64,6 +65,7 @@ Table of contents
     troubleshooting
     acknowledge
     license
+    changes
     changelog
     Contributing <https://github.com/MeteoSwiss/ampycloud/blob/develop/CONTRIBUTING.md>
     Github repository <https://github.com/MeteoSwiss/ampycloud>

setup.py

@@ -46,7 +46,7 @@
         "matplotlib >= 3.7.2",
         "numpy >= 1.20.3",
         "scikit-learn >= 1.2.0",
-        "scipy >= 1.7.3",
+        "scipy >= 1.7.3, < 1.14.1", # avoid conflict with statsmodels 0.14.1
         "statsmodels",
         "pandas >= 1.5",
         "pyyaml",

src/ampycloud/data.py

@@ -13,11 +13,12 @@
 import logging
 import copy
 from abc import ABC, abstractmethod
+import warnings
 import numpy as np
 import pandas as pd
 
 # Import from this package
-from .errors import AmpycloudError
+from .errors import AmpycloudError, AmpycloudWarning
 from .logger import log_func_call
 from . import scaler, cluster, layer, fluffer
 from . import wmo, icao
@@ -528,6 +529,20 @@ def _calculate_sligrolay_base_height(
         for ind, cid in enumerate(cluster_ids):
             # Which hits are in this sli/gro/lay ?
             in_sligrolay = self.data[which[:-1]+'_id'] == cid
+            if self.prms['CEILOS_FOR_BASE_HEIGHT_CALC'] != []:
+                in_sligrolay_filtered = in_sligrolay * self.data['ceilo'].apply(
+                    lambda x: x in self.prms['CEILOS_FOR_BASE_HEIGHT_CALC']
+                )
+                # We require a minimum of hits by the filtered ceilos that belong to the layer
+                # of interest. Otherwise fall back to using all ceilos for the calculation.
+                if in_sligrolay_filtered.sum() > self.prms['MAX_HITS_OKTA0']:
+                    in_sligrolay = in_sligrolay_filtered
+                else:
+                    warnings.warn(
+                        'Not enough data after filtering to calculate cloud base height, '
+                        f'will fall back to use all data in group/ slice/ layer {cid}',
+                        AmpycloudWarning
+                    )
             # Compute the base height
             pdf.iloc[ind, pdf.columns.get_loc('height_base')] = \
                 self._calculate_base_height_for_selection(in_sligrolay)

src/ampycloud/prms/ampycloud_default_prms.yml

@@ -38,6 +38,13 @@ BASE_LVL_HEIGHT_PERC: 5
 # Set to 100 to use all the points.
 BASE_LVL_LOOKBACK_PERC: 100
 
+# If ceilo Ids are set, for any given layer/ group/ slice, only ceilos with the corresponding ID
+# will be used for base height calculation. If the selection results in an empty layer/ group/
+# slice, we go back to using all ceilo hits and ignore the selection. Use this switch if you want
+# to localize the height caclulation or if you only want ceilometers from a specific type in the
+# height calculation.
+CEILOS_FOR_BASE_HEIGHT_CALC: []
+
 # LOWESS parameters used to fit Slice/Groups/Layers height trends
 LOWESS:
     # Fraction of the slice/group/layer points to consider when deriving the LOWESS smooth fit.

test/ampycloud/test_data.py

@@ -13,6 +13,7 @@
 import numpy as np
 import pandas as pd
 from pytest import raises, warns, mark, param
+import pytest
 
 # Import from the module to test
 from ampycloud.errors import AmpycloudError, AmpycloudWarning
@@ -413,3 +414,56 @@ def test_layering_dualeval():
         warnings.simplefilter("error")
         warnings.simplefilter("default", category=FutureWarning)  # Fixes #87
         chunk.find_layers()
+
+
+def _mock_layer_from_two_ceilos_with_narrow_stds_and_close_hits(n1: int, n2: int):
+    dts1 = np.linspace(-900., 0., n1)
+    layer_ceilo1 = mocker.flat_layer(dts1, 310., 2., 1.0)
+    layer_ceilo1['ceilo'] = np.array(['ceilo1'] * len(layer_ceilo1), dtype=str)
+    dts2 = np.linspace(-908., -8., n2)
+    layer_ceilo2 = mocker.flat_layer(dts2, 460., 2., 1.0)
+    layer_ceilo2['ceilo'] = np.array(['ceilo2'] * len(layer_ceilo2), dtype=str)
+    combined_layer = pd.concat([layer_ceilo1, layer_ceilo2])
+    combined_layer['type'] = np.ones_like(combined_layer['height'], dtype=int)
+
+    return combined_layer.sort_values('dt')
+
+
+def test_localized_base_height_calculation():
+    """Test the localized base height calculation. """
+
+    dynamic.AMPYCLOUD_PRMS['BASE_LVL_HEIGHT_PERC'] = 50
+    dynamic.AMPYCLOUD_PRMS['MSA'] = 10000
+    dynamic.AMPYCLOUD_PRMS['MIN_SEP_VALS'] = [250, 1000]  # make sure there is only 1 layer.
+
+    mock_layer = _mock_layer_from_two_ceilos_with_narrow_stds_and_close_hits(100, 200)
+    chunk_filtered = CeiloChunk(
+        mock_layer, prms={'CEILOS_FOR_BASE_HEIGHT_CALC': 'ceilo1'}
+    )
+    chunk_filtered.find_slices()
+    chunk_filtered.find_groups()
+    chunk_filtered.find_layers()
+
+    # Since the SD of the two layers is 2ft, it's highly unlikely that the following condition
+    # 1) is not met if the filtering works (310+/-2 ft)
+    # 2) is met if the filtering doesn't (in total twice as many point at 460 +/-2 ft)
+    assert chunk_filtered.layers.loc[0, 'height_base'] < 330.
+
+
+def test_localized_base_hgt_calc_fallback():
+    """Test if the localized base height calculation falls back if there are too
+    few hits after filtering by ceilo."""
+
+    dynamic.AMPYCLOUD_PRMS['MAX_HITS_OKTA0'] = 3
+    dynamic.AMPYCLOUD_PRMS['BASE_LVL_HEIGHT_PERC'] = 1
+    dynamic.AMPYCLOUD_PRMS['MSA'] = 10000
+    dynamic.AMPYCLOUD_PRMS['MIN_SEP_VALS'] = [250, 1000]  # make sure there is only 1 layer.
+
+    mock_layer = _mock_layer_from_two_ceilos_with_narrow_stds_and_close_hits(2, 200)
+    chunk_filtered = CeiloChunk(
+        mock_layer, prms={'CEILOS_FOR_BASE_HEIGHT_CALC': 'ceilo1'}
+    )
+    chunk_filtered.find_slices()
+    chunk_filtered.find_groups()
+    with pytest.warns(AmpycloudWarning, match="will fall back to use all data"):
+        chunk_filtered.find_layers()