Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Clip negative FCI radiances #3013

Open
wants to merge 14 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion doc/source/config.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
12 changes: 11 additions & 1 deletion satpy/readers/fci_l1c_nc.py
Original file line number Diff line number Diff line change
Expand Up @@ -208,31 +208,33 @@
"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,
cache_var_size=0,
cache_handle=True)
logger.debug("Reading: {}".format(self.filename))
logger.debug("Start: {}".format(self.start_time))
logger.debug("End: {}".format(self.end_time))

if self.filename_info["coverage"] == "Q4":
# change the chunk number so that padding gets activated correctly for Q4, which corresponds to the upper
# quarter of the disc
self.filename_info["count_in_repeat_cycle"] += 28

if self.filename_info["coverage"] == "AF":
# change number of chunk from 0 to 1 so that the padding is not activated (chunk 1 is present and only 1
# chunk is expected), as the African dissemination products come in one file per full disk.
self.filename_info["count_in_repeat_cycle"] = 1

if self.filename_info["facility_or_tool"] == "IQTI":
self.is_iqt = True
else:
self.is_iqt = False

self.clip_negative_radiances = clip_negative_radiances

Check warning on line 237 in satpy/readers/fci_l1c_nc.py

View check run for this annotation

CodeScene Delta Analysis / CodeScene Cloud Delta Analysis (main)

❌ New issue: Excess Number of Function Arguments

FCIL1cNCFileHandler.__init__ has 5 arguments, threshold = 4. This function has too many arguments, indicating a lack of encapsulation. Avoid adding more arguments.
self._cache = {}

@property
Expand Down Expand Up @@ -661,6 +663,8 @@

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) +
Expand All @@ -677,6 +681,12 @@
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.notnull())|(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
Expand Down
16 changes: 9 additions & 7 deletions satpy/tests/behave/features/image_comparison.feature
Original file line number Diff line number Diff line change
@@ -1,13 +1,15 @@
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> with clipping <clip>
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 | 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 |
27 changes: 19 additions & 8 deletions satpy/tests/behave/features/steps/image_comparison.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,9 @@
# satpy. If not, see <http://www.gnu.org/licenses/>.
"""Image comparison tests."""

import hdf5plugin # noqa: F401 isort:skip
import os
import os.path
import warnings
from datetime import datetime
from glob import glob
Expand Down Expand Up @@ -51,34 +53,43 @@
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"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
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} 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"
warnings.simplefilter("ignore")
dask.config.set(scheduler="threads", num_workers=4)

# 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)
reader_kwargs = {}
if clip != "null":
reader_kwargs["clip_negative_radiances"] = clip
scn = Scene(reader=reader, filenames=filenames, reader_kwargs=reader_kwargs)

scn.load([composite])

if area == "null":
ls = scn
else:
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",
f"generated_{context.satellite}_{context.composite}.png")
scn.save_datasets(writer="simple_image", filename=generated_image_path)
f"generated_{context.satellite}_{context.composite}_{context.area}.png")
ls.save_datasets(writer="simple_image", filename=generated_image_path)

Check warning on line 92 in satpy/tests/behave/features/steps/image_comparison.py

View check run for this annotation

CodeScene Delta Analysis / CodeScene Cloud Delta Analysis (main)

❌ New issue: Excess Number of Function Arguments

step_when_generate_image has 6 arguments, threshold = 4. This function has too many arguments, indicating a lack of encapsulation. Avoid adding more arguments.

# Save the generated image in the context
context.generated_image = cv2.imread(generated_image_path)
Expand Down
41 changes: 32 additions & 9 deletions satpy/tests/reader_tests/test_fci_l1c_nc.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
#!/usr/bin/env python

Check notice on line 1 in satpy/tests/reader_tests/test_fci_l1c_nc.py

View check run for this annotation

CodeScene Delta Analysis / CodeScene Cloud Delta Analysis (main)

ℹ Getting worse: Lines of Code in a Single File

The lines of code increases from 632 to 648, improve code health by reducing it to 600. The number of Lines of Code in a single file. More Lines of Code lowers the code health.
# -*- coding: utf-8 -*-
# Copyright (c) 2019 Satpy developers
#
Expand Down Expand Up @@ -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)
Expand All @@ -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",
},
Expand All @@ -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",
},
Expand Down Expand Up @@ -293,16 +296,17 @@

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],
Expand All @@ -313,7 +317,7 @@
)
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],
Expand Down Expand Up @@ -440,9 +444,14 @@
"""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:
Expand Down Expand Up @@ -542,11 +551,11 @@
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

Expand Down Expand Up @@ -738,7 +747,8 @@
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"])
Expand Down Expand Up @@ -860,6 +870,19 @@
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_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

@pytest.mark.parametrize(("calibration", "channel", "resolution"), [
(calibration, channel, resolution)
for calibration in ["counts", "radiance", "brightness_temperature", "reflectance"]
Expand Down
95 changes: 77 additions & 18 deletions utils/create_reference.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,34 +18,93 @@

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 sys
from glob import glob
import argparse
import os
import pathlib

from dask.diagnostics import ProgressBar
import hdf5plugin # noqa: F401

from satpy import Scene

ext_data_path = sys.argv[1]
outdir = sys.argv[2]
satellite = sys.argv[3]

filenames = glob(f"{ext_data_path}/satellite_data/{satellite}/*.nc")
def generate_images(props):
"""Generate reference images for testing purposes.

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

from dask.diagnostics import ProgressBar
with ProgressBar():
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."""
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(
"-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(
"-o", "--outdir", action="store", type=pathlib.Path,
default=pathlib.Path("."),
help="Directory where to write resulting images.")

parser.add_argument(
"-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()

scn = Scene(reader="abi_l1b", filenames=filenames)
generate_images(parsed)

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()
Loading