From 576e94180b584a33706bc84c9f633016bdcb4130 Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Mon, 9 Dec 2024 16:28:17 +0100
Subject: [PATCH 01/14] Started behave test for FCI cloudtop

---
 .../behave/features/image_comparison.feature      | 15 ++++++++-------
 .../behave/features/steps/image_comparison.py     | 15 ++++++++-------
 2 files changed, 16 insertions(+), 14 deletions(-)

diff --git a/satpy/tests/behave/features/image_comparison.feature b/satpy/tests/behave/features/image_comparison.feature
index 3352db70b4..03064fe824 100755
--- a/satpy/tests/behave/features/image_comparison.feature
+++ b/satpy/tests/behave/features/image_comparison.feature
@@ -1,13 +1,14 @@
 Feature: Image Comparison
 
   Scenario Outline: Compare generated image with reference image
-    Given I have a <composite> reference image file from <satellite>
-    When I generate a new <composite> image file from <satellite>
+    Given I have a <composite> reference image file from <satellite> resampled to <area>
+    When I generate a new <composite> image file from <satellite> with <reader> for <area>
     Then the generated image should be the same as the reference image
 
     Examples:
-      |satellite |composite  |
-      |GOES17   |airmass  |
-      |GOES16   |airmass  |
-      |GOES16   |ash      |
-      |GOES17   |ash      |
+      |satellite |composite  | reader | area
+      |GOES17   |airmass  | abi_l1b | null
+      |GOES16   |airmass  | abi_l1b | null
+      |GOES16   |ash      | abi_l1b | null
+      |GOES17   |ash      | abi_l1b | null
+      |METEOSAT12 | cloudtop | fci_l1b_nc | north_atlantic
diff --git a/satpy/tests/behave/features/steps/image_comparison.py b/satpy/tests/behave/features/steps/image_comparison.py
index 17a17cdeeb..6c658ceaac 100644
--- a/satpy/tests/behave/features/steps/image_comparison.py
+++ b/satpy/tests/behave/features/steps/image_comparison.py
@@ -51,17 +51,18 @@ def setup_hooks():
     use_fixture(before_all, Context)
 
 setup_hooks()
-@given("I have a {composite} reference image file from {satellite}")
-def step_given_reference_image(context, composite, satellite):
+@given("I have a {composite} reference image file from {satellite} resampled to <area>")
+def step_given_reference_image(context, composite, satellite, area):
     """Prepare a reference image."""
-    reference_image = f"reference_image_{satellite}_{composite}.png"
+    reference_image = f"reference_image_{satellite}_{composite}_{area}.png"
     context.reference_image = cv2.imread(f"{ext_data_path}/reference_images/{reference_image}")
     context.satellite = satellite
     context.composite = composite
+    context.area = area
 
 
-@when("I generate a new {composite} image file from {satellite}")
-def step_when_generate_image(context, composite, satellite):
+@when("I generate a new {composite} image file from {satellite} with {reader} for {area}")
+def step_when_generate_image(context, composite, satellite, reader, area):
     """Generate test images."""
     os.environ["OMP_NUM_THREADS"] = os.environ["MKL_NUM_THREADS"] = "2"
     os.environ["PYTROLL_CHUNK_SIZE"] = "1024"
@@ -71,13 +72,13 @@ def step_when_generate_image(context, composite, satellite):
     # Get the list of satellite files to open
     filenames = glob(f"{ext_data_path}/satellite_data/{satellite}/*.nc")
 
-    scn = Scene(reader="abi_l1b", filenames=filenames)
+    scn = Scene(reader=reader, filenames=filenames)
 
     scn.load([composite])
 
     # Save the generated image in the generated folder
     generated_image_path = os.path.join(context.test_results_dir, "generated",
-                                        f"generated_{context.satellite}_{context.composite}.png")
+                                        f"generated_{context.satellite}_{context.composite}_{context.area}.png")
     scn.save_datasets(writer="simple_image", filename=generated_image_path)
 
     # Save the generated image in the context

From d4a1e20adbde3d2701bb77ddff2de4ef0aaeee6a Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Mon, 9 Dec 2024 18:57:10 +0100
Subject: [PATCH 02/14] Continue development of FCI imagery tests

---
 .../behave/features/steps/image_comparison.py | 18 ++++-
 utils/create_reference.py                     | 74 +++++++++++++++----
 2 files changed, 76 insertions(+), 16 deletions(-)

diff --git a/satpy/tests/behave/features/steps/image_comparison.py b/satpy/tests/behave/features/steps/image_comparison.py
index 6c658ceaac..3891e8f5dc 100644
--- a/satpy/tests/behave/features/steps/image_comparison.py
+++ b/satpy/tests/behave/features/steps/image_comparison.py
@@ -24,12 +24,21 @@
 import dask
 import numpy as np
 from behave import given, then, when
+from pyresample.area_config import create_area_def
 
 from satpy import Scene
 
 ext_data_path = "/app/ext_data"
 threshold = 2000
 
+# test areas used only for testing
+test_areas = {
+        "north_atlantic": create_area_def(
+            "ofz", 4087, description="oceanographer fracture zone",
+            area_extent=[-4230000, 4675000, -3562000, 5232000],
+            resolution=750)
+        }
+
 def before_all(context):
     """Define a before_all hook to create the timestamp and test results directory."""
     tm = datetime.now()
@@ -51,7 +60,7 @@ def setup_hooks():
     use_fixture(before_all, Context)
 
 setup_hooks()
-@given("I have a {composite} reference image file from {satellite} resampled to <area>")
+@given("I have a {composite} reference image file from {satellite} resampled to {area}")
 def step_given_reference_image(context, composite, satellite, area):
     """Prepare a reference image."""
     reference_image = f"reference_image_{satellite}_{composite}_{area}.png"
@@ -76,10 +85,15 @@ def step_when_generate_image(context, composite, satellite, reader, area):
 
     scn.load([composite])
 
+    if area == "null":
+        ls = scn
+    else:
+        ls = scn.resample(test_areas.get(area, area))
+
     # Save the generated image in the generated folder
     generated_image_path = os.path.join(context.test_results_dir, "generated",
                                         f"generated_{context.satellite}_{context.composite}_{context.area}.png")
-    scn.save_datasets(writer="simple_image", filename=generated_image_path)
+    ls.save_datasets(writer="simple_image", filename=generated_image_path)
 
     # Save the generated image in the context
     context.generated_image = cv2.imread(generated_image_path)
diff --git a/utils/create_reference.py b/utils/create_reference.py
index 8604cb45c0..489a37cddf 100644
--- a/utils/create_reference.py
+++ b/utils/create_reference.py
@@ -28,24 +28,70 @@
 imagery.
 """
 
-import sys
+import argparse
+import pathlib
 from glob import glob
 
-from dask.diagnostics import ProgressBar
-
 from satpy import Scene
+from satpy.tests.behave.features.steps.image_comparison import test_areas
+
+
+def generate_images(reader, filenames, area, composites, outdir):
+    """Generate reference images for testing purposes."""
+    from dask.diagnostics import ProgressBar
+    scn = Scene(reader="abi_l1b", filenames=filenames)
+
+    composites = ["ash", "airmass"]
+    scn.load(composites)
+    if area is None:
+        ls = scn
+    elif area == "native":
+        ls = scn.resample(resampler="native")
+    else:
+        ls = scn.resample(test_areas.get(area, area))
+
+    with ProgressBar():
+        ls.save_datasets(writer="simple_image", filename=outdir +
+                          "/satpy-reference-image-{platform_name}-{sensor}-{start_time:%Y%m%d%H%M}-{area.area_id}-{name}.png")
+
+def get_parser():
+    """Return argument parser."""
+    parser = argparse.ArgumentParser(description=__doc__)
+
+    parser.add_argument(
+            "satellite", action="store", type=str,
+            help="Satellite name.")
+
+    parser.add_argument(
+            "reader", action="store", type=str,
+            help="Reader name.")
+
+    parser.add_argument(
+            "area", action="store", type=str,
+            help="Area name, 'null' (no resampling) or 'native' (native resampling)")
+
+    parser.add_argument(
+            "basedir", action="store", type=pathlib.Path,
+            help="Root directory where reference input data are contained.")
+
+    parser.add_argument(
+            "outdir", action="store", type=pathlib.Path,
+            help="Directory where to write resulting images.")
 
-ext_data_path = sys.argv[1]
-outdir = sys.argv[2]
-satellite = sys.argv[3]
+    return parser
 
-filenames = glob(f"{ext_data_path}/satellite_data/{satellite}/*.nc")
+def main():
+    """Main function."""
+    parsed = get_parser().parse_args()
+    ext_data_path = parsed.basedir
+    reader = parsed.reader
+    area = parsed.area
+    outdir = parsed.outdir
+    satellite = parsed.satellite
 
-scn = Scene(reader="abi_l1b", filenames=filenames)
+    filenames = glob(f"{ext_data_path}/satellite_data/{satellite}/*")
+    generate_images(reader, filenames, None if area.lower() == "null" else
+                    area, ["airmass", "ash"], outdir)
 
-composites = ["ash", "airmass"]
-scn.load(composites)
-ls = scn.resample(resampler="native")
-with ProgressBar():
-    ls.save_datasets(writer="simple_image", filename=outdir +
-                      "/satpy-reference-image-{platform_name}-{sensor}-{start_time:%Y%m%d%H%M}-{area.area_id}-{name}.png")
+if __name__ == "__main__":
+    main()

From 1f562ae413bbf0871d98be6680be7fe28870f4db Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Tue, 10 Dec 2024 10:01:08 +0100
Subject: [PATCH 03/14] Change reference data and area for ABI test

---
 satpy/tests/behave/features/image_comparison.feature  |  2 +-
 satpy/tests/behave/features/steps/image_comparison.py | 11 +----------
 utils/create_reference.py                             |  3 +--
 3 files changed, 3 insertions(+), 13 deletions(-)

diff --git a/satpy/tests/behave/features/image_comparison.feature b/satpy/tests/behave/features/image_comparison.feature
index 03064fe824..5bb6c9c2c9 100755
--- a/satpy/tests/behave/features/image_comparison.feature
+++ b/satpy/tests/behave/features/image_comparison.feature
@@ -11,4 +11,4 @@ Feature: Image Comparison
       |GOES16   |airmass  | abi_l1b | null
       |GOES16   |ash      | abi_l1b | null
       |GOES17   |ash      | abi_l1b | null
-      |METEOSAT12 | cloudtop | fci_l1b_nc | north_atlantic
+      |METEOSAT12 | cloudtop | fci_l1b_nc | sve
diff --git a/satpy/tests/behave/features/steps/image_comparison.py b/satpy/tests/behave/features/steps/image_comparison.py
index 3891e8f5dc..74f0af3229 100644
--- a/satpy/tests/behave/features/steps/image_comparison.py
+++ b/satpy/tests/behave/features/steps/image_comparison.py
@@ -24,21 +24,12 @@
 import dask
 import numpy as np
 from behave import given, then, when
-from pyresample.area_config import create_area_def
 
 from satpy import Scene
 
 ext_data_path = "/app/ext_data"
 threshold = 2000
 
-# test areas used only for testing
-test_areas = {
-        "north_atlantic": create_area_def(
-            "ofz", 4087, description="oceanographer fracture zone",
-            area_extent=[-4230000, 4675000, -3562000, 5232000],
-            resolution=750)
-        }
-
 def before_all(context):
     """Define a before_all hook to create the timestamp and test results directory."""
     tm = datetime.now()
@@ -88,7 +79,7 @@ def step_when_generate_image(context, composite, satellite, reader, area):
     if area == "null":
         ls = scn
     else:
-        ls = scn.resample(test_areas.get(area, area))
+        ls = scn.resample(area)
 
     # Save the generated image in the generated folder
     generated_image_path = os.path.join(context.test_results_dir, "generated",
diff --git a/utils/create_reference.py b/utils/create_reference.py
index 489a37cddf..ca0d0a180d 100644
--- a/utils/create_reference.py
+++ b/utils/create_reference.py
@@ -33,7 +33,6 @@
 from glob import glob
 
 from satpy import Scene
-from satpy.tests.behave.features.steps.image_comparison import test_areas
 
 
 def generate_images(reader, filenames, area, composites, outdir):
@@ -48,7 +47,7 @@ def generate_images(reader, filenames, area, composites, outdir):
     elif area == "native":
         ls = scn.resample(resampler="native")
     else:
-        ls = scn.resample(test_areas.get(area, area))
+        ls = scn.resample(area)
 
     with ProgressBar():
         ls.save_datasets(writer="simple_image", filename=outdir +

From 2297ee7c5d40f19598ba1a77d2c531b308bb9668 Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Fri, 13 Dec 2024 09:35:46 +0100
Subject: [PATCH 04/14] Improve reference data generation script.

---
 utils/create_reference.py | 78 +++++++++++++++++++++++----------------
 1 file changed, 46 insertions(+), 32 deletions(-)

diff --git a/utils/create_reference.py b/utils/create_reference.py
index ca0d0a180d..5510054099 100644
--- a/utils/create_reference.py
+++ b/utils/create_reference.py
@@ -18,40 +18,50 @@
 
 Script to create reference images for the automated image testing system.
 
-create_reference.py <input-data> <output-directory> <satellite-name>
-
 The input data directory must follow the data structure from the
 image-comparison-tests repository with satellite_data/<satellite-name>.
 
 This script is a work in progress and expected to change significantly.
-It is absolutely not intended for any operational production of satellite
-imagery.
+
+DO NOT USE FOR OPERATIONAL PRODUCTION!
 """
 
 import argparse
+import os
 import pathlib
-from glob import glob
+
+import hdf5plugin  # noqa: F401
 
 from satpy import Scene
 
 
-def generate_images(reader, filenames, area, composites, outdir):
-    """Generate reference images for testing purposes."""
-    from dask.diagnostics import ProgressBar
-    scn = Scene(reader="abi_l1b", filenames=filenames)
+def generate_images(props):
+    """Generate reference images for testing purposes.
 
-    composites = ["ash", "airmass"]
-    scn.load(composites)
-    if area is None:
-        ls = scn
-    elif area == "native":
+    Args:
+        props (namespace): Object with attributes corresponding to command line
+        arguments as defined by :func:get_parser.
+    """
+    filenames = (props.basedir / "satellite_data" / props.satellite).glob("*")
+
+    scn = Scene(reader=props.reader, filenames=filenames)
+
+    scn.load(props.composites)
+    if props.area == "native":
         ls = scn.resample(resampler="native")
+    elif props.area is not None:
+        ls = scn.resample(props.area, resampler="gradient_search")
     else:
-        ls = scn.resample(area)
+        ls = scn
 
+    from dask.diagnostics import ProgressBar
     with ProgressBar():
-        ls.save_datasets(writer="simple_image", filename=outdir +
-                          "/satpy-reference-image-{platform_name}-{sensor}-{start_time:%Y%m%d%H%M}-{area.area_id}-{name}.png")
+        ls.save_datasets(
+                writer="simple_image",
+                filename=os.fspath(
+                    props.basedir / "reference_images" /
+                    "satpy-reference-image-{platform_name}-{sensor}-"
+                    "{start_time:%Y%m%d%H%M}-{area.area_id}-{name}.png"))
 
 def get_parser():
     """Return argument parser."""
@@ -66,31 +76,35 @@ def get_parser():
             help="Reader name.")
 
     parser.add_argument(
-            "area", action="store", type=str,
-            help="Area name, 'null' (no resampling) or 'native' (native resampling)")
+            "-b", "--basedir", action="store", type=pathlib.Path,
+            default=pathlib.Path("."),
+            help="Base directory for reference data. "
+                 "This must contain a subdirectories satellite_data and "
+                 "reference_images.  The directory satellite_data must contain "
+                 "input data in a subdirectory for the satellite. Output images "
+                 "will be written to the subdirectory reference_images.")
 
     parser.add_argument(
-            "basedir", action="store", type=pathlib.Path,
-            help="Root directory where reference input data are contained.")
+            "-o", "--outdir", action="store", type=pathlib.Path,
+            default=pathlib.Path("."),
+            help="Directory where to write resulting images.")
 
     parser.add_argument(
-            "outdir", action="store", type=pathlib.Path,
-            help="Directory where to write resulting images.")
+            "-c", "--composites", nargs="+", help="composites to generate",
+            type=str, default=["ash", "airmass"])
+
+    parser.add_argument(
+            "-a", "--area", action="store",
+            default=None,
+            help="Area name, or 'native' (native resampling)")
 
     return parser
 
 def main():
     """Main function."""
     parsed = get_parser().parse_args()
-    ext_data_path = parsed.basedir
-    reader = parsed.reader
-    area = parsed.area
-    outdir = parsed.outdir
-    satellite = parsed.satellite
-
-    filenames = glob(f"{ext_data_path}/satellite_data/{satellite}/*")
-    generate_images(reader, filenames, None if area.lower() == "null" else
-                    area, ["airmass", "ash"], outdir)
+
+    generate_images(parsed)
 
 if __name__ == "__main__":
     main()

From 02989b9ad3cf850a66cb1577be880d050bb9f38f Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Fri, 13 Dec 2024 12:55:31 +0100
Subject: [PATCH 05/14] Fix feature syntax and filename setup

---
 .../tests/behave/features/image_comparison.feature  | 13 +++++++------
 .../tests/behave/features/steps/image_comparison.py |  3 ++-
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/satpy/tests/behave/features/image_comparison.feature b/satpy/tests/behave/features/image_comparison.feature
index 5bb6c9c2c9..d367732ae3 100755
--- a/satpy/tests/behave/features/image_comparison.feature
+++ b/satpy/tests/behave/features/image_comparison.feature
@@ -6,9 +6,10 @@ Feature: Image Comparison
     Then the generated image should be the same as the reference image
 
     Examples:
-      |satellite |composite  | reader | area
-      |GOES17   |airmass  | abi_l1b | null
-      |GOES16   |airmass  | abi_l1b | null
-      |GOES16   |ash      | abi_l1b | null
-      |GOES17   |ash      | abi_l1b | null
-      |METEOSAT12 | cloudtop | fci_l1b_nc | sve
+      |satellite |composite  | reader | area |
+      |GOES17   |airmass  | abi_l1b | null |
+      |GOES16   |airmass  | abi_l1b | null |
+      |GOES16   |ash      | abi_l1b | null |
+      |GOES17   |ash      | abi_l1b | null |
+      |Meteosat-12 | cloudtop | fci_l1c_nc | sve |
+      |Meteosat-12 | night_microphysics | fci_l1c_nc | sve |
diff --git a/satpy/tests/behave/features/steps/image_comparison.py b/satpy/tests/behave/features/steps/image_comparison.py
index 74f0af3229..6e33baeb1e 100644
--- a/satpy/tests/behave/features/steps/image_comparison.py
+++ b/satpy/tests/behave/features/steps/image_comparison.py
@@ -22,6 +22,7 @@
 
 import cv2
 import dask
+import hdf5plugin  # noqa: F401
 import numpy as np
 from behave import given, then, when
 
@@ -54,7 +55,7 @@ def setup_hooks():
 @given("I have a {composite} reference image file from {satellite} resampled to {area}")
 def step_given_reference_image(context, composite, satellite, area):
     """Prepare a reference image."""
-    reference_image = f"reference_image_{satellite}_{composite}_{area}.png"
+    reference_image = f"satpy-reference-image-{satellite}-{composite}-{area}.png"
     context.reference_image = cv2.imread(f"{ext_data_path}/reference_images/{reference_image}")
     context.satellite = satellite
     context.composite = composite

From 8ce9a62535769c0e4d4171392d025e2541f69fe4 Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Fri, 13 Dec 2024 14:27:22 +0100
Subject: [PATCH 06/14] import hdf5plugin first despite isort

---
 satpy/tests/behave/features/steps/image_comparison.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/satpy/tests/behave/features/steps/image_comparison.py b/satpy/tests/behave/features/steps/image_comparison.py
index 6e33baeb1e..1349e0b64e 100644
--- a/satpy/tests/behave/features/steps/image_comparison.py
+++ b/satpy/tests/behave/features/steps/image_comparison.py
@@ -15,6 +15,7 @@
 # satpy.  If not, see <http://www.gnu.org/licenses/>.
 """Image comparison tests."""
 
+import hdf5plugin  # noqa: F401  isort:skip
 import os
 import warnings
 from datetime import datetime
@@ -22,7 +23,6 @@
 
 import cv2
 import dask
-import hdf5plugin  # noqa: F401
 import numpy as np
 from behave import given, then, when
 
@@ -80,7 +80,7 @@ def step_when_generate_image(context, composite, satellite, reader, area):
     if area == "null":
         ls = scn
     else:
-        ls = scn.resample(area)
+        ls = scn.resample(area, resampler="gradient_search")
 
     # Save the generated image in the generated folder
     generated_image_path = os.path.join(context.test_results_dir, "generated",

From 7cd8858d6f887a259d66737ed2264bf1a73b334e Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Fri, 13 Dec 2024 14:42:16 +0100
Subject: [PATCH 07/14] Does it help if I do Meteosat-12 first?

---
 satpy/tests/behave/features/image_comparison.feature | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/satpy/tests/behave/features/image_comparison.feature b/satpy/tests/behave/features/image_comparison.feature
index d367732ae3..b794fa22ef 100755
--- a/satpy/tests/behave/features/image_comparison.feature
+++ b/satpy/tests/behave/features/image_comparison.feature
@@ -7,9 +7,9 @@ Feature: Image Comparison
 
     Examples:
       |satellite |composite  | reader | area |
+      |Meteosat-12 | cloudtop | fci_l1c_nc | sve |
+      |Meteosat-12 | night_microphysics | fci_l1c_nc | sve |
       |GOES17   |airmass  | abi_l1b | null |
       |GOES16   |airmass  | abi_l1b | null |
       |GOES16   |ash      | abi_l1b | null |
       |GOES17   |ash      | abi_l1b | null |
-      |Meteosat-12 | cloudtop | fci_l1c_nc | sve |
-      |Meteosat-12 | night_microphysics | fci_l1c_nc | sve |

From 4d13d122b0e5068c85f641d4d1a23fe801641971 Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Fri, 13 Dec 2024 15:58:21 +0100
Subject: [PATCH 08/14] workaround for failing hdf5lpguni?

Trying John Cintineos workaround.
---
 satpy/tests/behave/features/steps/image_comparison.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/satpy/tests/behave/features/steps/image_comparison.py b/satpy/tests/behave/features/steps/image_comparison.py
index 1349e0b64e..a35dd22c00 100644
--- a/satpy/tests/behave/features/steps/image_comparison.py
+++ b/satpy/tests/behave/features/steps/image_comparison.py
@@ -17,6 +17,10 @@
 
 import hdf5plugin  # noqa: F401  isort:skip
 import os
+import os.path
+
+os.environ["HDF5_PLUGIN_PATH"] = os.path.dirname(hdf5plugin.__file__) + "/plugins/"
+
 import warnings
 from datetime import datetime
 from glob import glob

From 4c912ebe7d47368451227d6e00e00eed5241a60d Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Fri, 13 Dec 2024 16:36:47 +0100
Subject: [PATCH 09/14] switch on debug logging for troubleshooting purposes

---
 satpy/tests/behave/features/steps/image_comparison.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/satpy/tests/behave/features/steps/image_comparison.py b/satpy/tests/behave/features/steps/image_comparison.py
index a35dd22c00..a7afd56008 100644
--- a/satpy/tests/behave/features/steps/image_comparison.py
+++ b/satpy/tests/behave/features/steps/image_comparison.py
@@ -32,6 +32,8 @@
 
 from satpy import Scene
 
+from satpy.utils import debug_on; debug_on()
+
 ext_data_path = "/app/ext_data"
 threshold = 2000
 

From 50a6dfc67de9069049a64ceb6cf95f0c4696d0a2 Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Mon, 16 Dec 2024 17:16:05 +0100
Subject: [PATCH 10/14] Add unit test for FCI radiance clipping

Add a unit test for FCI radiance clipping.

The unit test fails, because there is no implementation yet.
---
 satpy/tests/reader_tests/test_fci_l1c_nc.py | 34 ++++++++++++++++-----
 1 file changed, 26 insertions(+), 8 deletions(-)

diff --git a/satpy/tests/reader_tests/test_fci_l1c_nc.py b/satpy/tests/reader_tests/test_fci_l1c_nc.py
index aa98990df3..8365205779 100644
--- a/satpy/tests/reader_tests/test_fci_l1c_nc.py
+++ b/satpy/tests/reader_tests/test_fci_l1c_nc.py
@@ -121,6 +121,7 @@
 DICT_CALIBRATION = {"radiance": {"dtype": np.float32,
                                  "value_1": 15,
                                  "value_0": 9700,
+                                 "value_2": -5,
                                  "attrs_dict": {"calibration": "radiance",
                                                 "units": "mW m-2 sr-1 (cm-1)-1",
                                                 "radiance_unit_conversion_coefficient": np.float32(1234.56)
@@ -134,8 +135,9 @@
                                     },
 
                     "counts": {"dtype": np.uint16,
-                               "value_1": 1,
+                               "value_1": 5,
                                "value_0": 5000,
+                               "value_2": 1,
                                "attrs_dict": {"calibration": "counts",
                                               "units": "count",
                                               },
@@ -144,6 +146,7 @@
                     "brightness_temperature": {"dtype": np.float32,
                                                "value_1": np.float32(209.68275),
                                                "value_0": np.float32(1888.8513),
+                                               "value_2": np.float32("nan"),
                                                "attrs_dict": {"calibration": "brightness_temperature",
                                                               "units": "K",
                                                               },
@@ -293,16 +296,17 @@ def _get_test_image_data_for_channel(data, ch_str, n_rows_cols):
 
     common_attrs = {
         "scale_factor": 5,
-        "add_offset": 10,
+        "add_offset": -10,
         "long_name": "Effective Radiance",
         "units": "mW.m-2.sr-1.(cm-1)-1",
         "ancillary_variables": "pixel_quality"
     }
     if "38" in ch_path:
         fire_line = da.ones((1, n_rows_cols[1]), dtype="uint16", chunks=1024) * 5000
-        data_without_fires = da.ones((n_rows_cols[0] - 1, n_rows_cols[1]), dtype="uint16", chunks=1024)
+        data_without_fires = da.full((n_rows_cols[0] - 2, n_rows_cols[1]), 5, dtype="uint16", chunks=1024)
+        neg_rad = da.ones((1, n_rows_cols[1]), dtype="uint16", chunks=1024)
         d = FakeH5Variable(
-            da.concatenate([fire_line, data_without_fires], axis=0),
+            da.concatenate([fire_line, data_without_fires, neg_rad], axis=0),
             dims=("y", "x"),
             attrs={
                 "valid_range": [0, 8191],
@@ -313,7 +317,7 @@ def _get_test_image_data_for_channel(data, ch_str, n_rows_cols):
         )
     else:
         d = FakeH5Variable(
-            da.ones(n_rows_cols, dtype="uint16", chunks=1024),
+            da.full(n_rows_cols, 5, dtype="uint16", chunks=1024),
             dims=("y", "x"),
             attrs={
                 "valid_range": [0, 4095],
@@ -542,11 +546,11 @@ def reader_configs():
         os.path.join("readers", "fci_l1c_nc.yaml"))
 
 
-def _get_reader_with_filehandlers(filenames, reader_configs):
+def _get_reader_with_filehandlers(filenames, reader_configs, **reader_kwargs):
     from satpy.readers import load_reader
     reader = load_reader(reader_configs)
     loadables = reader.select_files_from_pathnames(filenames)
-    reader.create_filehandlers(loadables)
+    reader.create_filehandlers(loadables, fh_kwargs=reader_kwargs)
     clear_cache(reader)
     return reader
 
@@ -738,7 +742,8 @@ def _reflectance_test(tab, filenames):
     def _other_calibration_test(res, ch, dict_arg):
         """Test of other calibration test."""
         if ch == "ir_38":
-            numpy.testing.assert_array_equal(res[ch][-1], dict_arg["value_1"])
+            numpy.testing.assert_array_equal(res[ch][-1], dict_arg["value_2"])
+            numpy.testing.assert_array_equal(res[ch][-2], dict_arg["value_1"])
             numpy.testing.assert_array_equal(res[ch][0], dict_arg["value_0"])
         else:
             numpy.testing.assert_array_equal(res[ch], dict_arg["value_1"])
@@ -860,6 +865,19 @@ def test_load_calibration(self, reader_configs, fh_param,
             self._get_assert_load(res, ch, DICT_CALIBRATION[calibration],
                                   fh_param["filenames"][0])
 
+    @pytest.mark.parametrize("fh_param", [lazy_fixture("FakeFCIFileHandlerFDHSI_fixture")])
+    def test_load_calibration_negative_rad(self, reader_configs, fh_param):
+        """Test calibrating negative radiances.
+
+        See https://github.com/pytroll/satpy/issues/3009.
+        """
+        reader = _get_reader_with_filehandlers(fh_param["filenames"],
+                                               reader_configs,
+                                               clip_negative_radiance=True)
+        res = reader.load([make_dataid(name="ir_38", calibration="radiance")],
+                           pad_data=False)
+        numpy.testing.assert_array_equal(res["ir_38"][-1, :], 5)  # smallest positive radiance
+
     @pytest.mark.parametrize(("calibration", "channel", "resolution"), [
         (calibration, channel, resolution)
         for calibration in ["counts", "radiance", "brightness_temperature", "reflectance"]

From 5d0675bd2d7e05728d6e1173de906375f2a79d13 Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Mon, 16 Dec 2024 17:58:05 +0100
Subject: [PATCH 11/14] Clip negative radiances

Clip negative radiances when a keyword arguments asks for it.
---
 satpy/readers/fci_l1c_nc.py                 | 12 +++++++++++-
 satpy/readers/yaml_reader.py                |  3 +--
 satpy/tests/reader_tests/test_fci_l1c_nc.py |  9 +++++++--
 3 files changed, 19 insertions(+), 5 deletions(-)

diff --git a/satpy/readers/fci_l1c_nc.py b/satpy/readers/fci_l1c_nc.py
index fc40916699..caad045f90 100644
--- a/satpy/readers/fci_l1c_nc.py
+++ b/satpy/readers/fci_l1c_nc.py
@@ -208,7 +208,8 @@ class using the :mod:`~satpy.Scene.load` method with the reader
         "MTI3": "MTG-I3",
         "MTI4": "MTG-I4"}
 
-    def __init__(self, filename, filename_info, filetype_info):
+    def __init__(self, filename, filename_info, filetype_info,
+                 clip_negative_radiances=False, **kwargs):
         """Initialize file handler."""
         super().__init__(filename, filename_info,
                          filetype_info,
@@ -233,6 +234,7 @@ def __init__(self, filename, filename_info, filetype_info):
         else:
             self.is_iqt = False
 
+        self.clip_negative_radiances = clip_negative_radiances
         self._cache = {}
 
     @property
@@ -661,6 +663,8 @@ def calibrate_counts_to_physical_quantity(self, data, key):
 
     def calibrate_counts_to_rad(self, data, key):
         """Calibrate counts to radiances."""
+        if self.clip_negative_radiances:
+            data = self._clipneg(data)
         if key["name"] == "ir_38":
             data = xr.where(((2 ** 12 - 1 < data) & (data <= 2 ** 13 - 1)),
                             (data * data.attrs.get("warm_scale_factor", 1) +
@@ -677,6 +681,12 @@ def calibrate_counts_to_rad(self, data, key):
                                self.get_and_cache_npxr(measured + "/radiance_unit_conversion_coefficient")})
         return data
 
+    @staticmethod
+    def _clipneg(data):
+        """Clip counts to avoid negative radiances."""
+        lo = -data.attrs.get("add_offset", 0) // data.attrs.get("scale_factor", 1) + 1
+        return data.where(data>=lo, lo)
+
     def calibrate_rad_to_bt(self, radiance, key):
         """IR channel calibration."""
         # using the method from PUG section Converting from Effective Radiance to Brightness Temperature for IR Channels
diff --git a/satpy/readers/yaml_reader.py b/satpy/readers/yaml_reader.py
index 5bbaba4a6c..ef371a6284 100644
--- a/satpy/readers/yaml_reader.py
+++ b/satpy/readers/yaml_reader.py
@@ -474,8 +474,7 @@ class FileYAMLReader(GenericYAMLReader, DataDownloadMixin):
     def __init__(self,
                  config_dict,
                  filter_parameters=None,
-                 filter_filenames=True,
-                 **kwargs):
+                 filter_filenames=True):
         """Set up initial internal storage for loading file data."""
         super().__init__(config_dict, filter_parameters, filter_filenames)
 
diff --git a/satpy/tests/reader_tests/test_fci_l1c_nc.py b/satpy/tests/reader_tests/test_fci_l1c_nc.py
index 8365205779..718d51c819 100644
--- a/satpy/tests/reader_tests/test_fci_l1c_nc.py
+++ b/satpy/tests/reader_tests/test_fci_l1c_nc.py
@@ -444,9 +444,14 @@ class FakeFCIFileHandlerBase(FakeNetCDF4FileHandler):
     """Class for faking the NetCDF4 Filehandler."""
 
     cached_file_content: Dict[str, xr.DataArray] = {}
-    # overwritten by FDHSI and HRFI FIle Handlers
+    # overwritten by FDHSI and HRFI File Handlers
     chan_patterns: Dict[str, Dict[str, Union[List[int], str]]] = {}
 
+    def __init__(self, *args, **kwargs):
+        """Initiative fake file handler."""
+        kwargs.pop("clip_negative_radiances", None)
+        super().__init__(*args, **kwargs)
+
     def _get_test_content_all_channels(self):
         data = {}
         for pat in self.chan_patterns:
@@ -873,7 +878,7 @@ def test_load_calibration_negative_rad(self, reader_configs, fh_param):
         """
         reader = _get_reader_with_filehandlers(fh_param["filenames"],
                                                reader_configs,
-                                               clip_negative_radiance=True)
+                                               clip_negative_radiances=True)
         res = reader.load([make_dataid(name="ir_38", calibration="radiance")],
                            pad_data=False)
         numpy.testing.assert_array_equal(res["ir_38"][-1, :], 5)  # smallest positive radiance

From aae848b146d10b5f424acad626035b9bb12f3926 Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Mon, 16 Dec 2024 18:22:56 +0100
Subject: [PATCH 12/14] Clip in image comparison tests

---
 .../behave/features/image_comparison.feature     | 16 ++++++++--------
 .../behave/features/steps/image_comparison.py    | 14 ++++++--------
 2 files changed, 14 insertions(+), 16 deletions(-)

diff --git a/satpy/tests/behave/features/image_comparison.feature b/satpy/tests/behave/features/image_comparison.feature
index b794fa22ef..0497a96c93 100755
--- a/satpy/tests/behave/features/image_comparison.feature
+++ b/satpy/tests/behave/features/image_comparison.feature
@@ -2,14 +2,14 @@ Feature: Image Comparison
 
   Scenario Outline: Compare generated image with reference image
     Given I have a <composite> reference image file from <satellite> resampled to <area>
-    When I generate a new <composite> image file from <satellite> with <reader> for <area>
+    When I generate a new <composite> image file from <satellite> with <reader> for <area> with clipping <clip>
     Then the generated image should be the same as the reference image
 
     Examples:
-      |satellite |composite  | reader | area |
-      |Meteosat-12 | cloudtop | fci_l1c_nc | sve |
-      |Meteosat-12 | night_microphysics | fci_l1c_nc | sve |
-      |GOES17   |airmass  | abi_l1b | null |
-      |GOES16   |airmass  | abi_l1b | null |
-      |GOES16   |ash      | abi_l1b | null |
-      |GOES17   |ash      | abi_l1b | null |
+      |satellite |composite  | reader | area | clip |
+      |Meteosat-12 | cloudtop | fci_l1c_nc | sve | True |
+      |Meteosat-12 | night_microphysics | fci_l1c_nc | sve | True |
+      |GOES17   |airmass  | abi_l1b | null | null |
+      |GOES16   |airmass  | abi_l1b | null | null |
+      |GOES16   |ash      | abi_l1b | null | null |
+      |GOES17   |ash      | abi_l1b | null | null |
diff --git a/satpy/tests/behave/features/steps/image_comparison.py b/satpy/tests/behave/features/steps/image_comparison.py
index a7afd56008..92c5fa0034 100644
--- a/satpy/tests/behave/features/steps/image_comparison.py
+++ b/satpy/tests/behave/features/steps/image_comparison.py
@@ -18,9 +18,6 @@
 import hdf5plugin  # noqa: F401  isort:skip
 import os
 import os.path
-
-os.environ["HDF5_PLUGIN_PATH"] = os.path.dirname(hdf5plugin.__file__) + "/plugins/"
-
 import warnings
 from datetime import datetime
 from glob import glob
@@ -32,8 +29,6 @@
 
 from satpy import Scene
 
-from satpy.utils import debug_on; debug_on()
-
 ext_data_path = "/app/ext_data"
 threshold = 2000
 
@@ -68,8 +63,8 @@ def step_given_reference_image(context, composite, satellite, area):
     context.area = area
 
 
-@when("I generate a new {composite} image file from {satellite} with {reader} for {area}")
-def step_when_generate_image(context, composite, satellite, reader, area):
+@when("I generate a new {composite} image file from {satellite} with {reader} for {area} with clipping {clip}")
+def step_when_generate_image(context, composite, satellite, reader, area, clip):
     """Generate test images."""
     os.environ["OMP_NUM_THREADS"] = os.environ["MKL_NUM_THREADS"] = "2"
     os.environ["PYTROLL_CHUNK_SIZE"] = "1024"
@@ -79,7 +74,10 @@ def step_when_generate_image(context, composite, satellite, reader, area):
     # Get the list of satellite files to open
     filenames = glob(f"{ext_data_path}/satellite_data/{satellite}/*.nc")
 
-    scn = Scene(reader=reader, filenames=filenames)
+    reader_kwargs = {}
+    if clip != "null":
+        reader_kwargs["clip_negative_radiances"] = clip
+    scn = Scene(reader=reader, filenames=filenames, reader_kwargs=reader_kwargs)
 
     scn.load([composite])
 

From 5f6598b6d2db00add68ca7c1a207d0659c497d45 Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Mon, 16 Dec 2024 18:27:52 +0100
Subject: [PATCH 13/14] revert erroneous removal of **kwargs

---
 satpy/readers/yaml_reader.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/satpy/readers/yaml_reader.py b/satpy/readers/yaml_reader.py
index ef371a6284..5bbaba4a6c 100644
--- a/satpy/readers/yaml_reader.py
+++ b/satpy/readers/yaml_reader.py
@@ -474,7 +474,8 @@ class FileYAMLReader(GenericYAMLReader, DataDownloadMixin):
     def __init__(self,
                  config_dict,
                  filter_parameters=None,
-                 filter_filenames=True):
+                 filter_filenames=True,
+                 **kwargs):
         """Set up initial internal storage for loading file data."""
         super().__init__(config_dict, filter_parameters, filter_filenames)
 

From d404a795704b9cf628b6210a11673bdcc19f0a45 Mon Sep 17 00:00:00 2001
From: Gerrit Holl <gerrit.holl@dwd.de>
Date: Wed, 18 Dec 2024 12:06:52 +0100
Subject: [PATCH 14/14] add doc & don't clip space pxiels

---
 doc/source/config.rst       | 2 +-
 satpy/readers/fci_l1c_nc.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/source/config.rst b/doc/source/config.rst
index 9babc1abbf..1cbbbec2ed 100644
--- a/doc/source/config.rst
+++ b/doc/source/config.rst
@@ -272,7 +272,7 @@ If ``clip_negative_radiances=False``, pixels with negative radiances will have
 
 Clipping of negative radiances is currently implemented for the following readers:
 
-* ``abi_l1b``, ``ami_l1b``
+* ``abi_l1b``, ``ami_l1b``, ``fci_l1c_nc``
 
 
 Temporary Directory
diff --git a/satpy/readers/fci_l1c_nc.py b/satpy/readers/fci_l1c_nc.py
index caad045f90..c24c9b4849 100644
--- a/satpy/readers/fci_l1c_nc.py
+++ b/satpy/readers/fci_l1c_nc.py
@@ -685,7 +685,7 @@ def calibrate_counts_to_rad(self, data, key):
     def _clipneg(data):
         """Clip counts to avoid negative radiances."""
         lo = -data.attrs.get("add_offset", 0) // data.attrs.get("scale_factor", 1) + 1
-        return data.where(data>=lo, lo)
+        return data.where((~data.notnull())|(data>=lo), lo)
 
     def calibrate_rad_to_bt(self, radiance, key):
         """IR channel calibration."""