Merge pull request #2687 from yufeizhu600/olci_l1b_quality_flags

add support of masking olci l1b products by using quality flags

satpy/etc/readers/olci_l1b.yaml

@@ -36,6 +36,13 @@ file_types:
         file_patterns:
          - '{mission_id:3s}_OL_1_{datatype_id:_<6s}_{start_time:%Y%m%dT%H%M%S}_{end_time:%Y%m%dT%H%M%S}_{creation_time:%Y%m%dT%H%M%S}_{duration:4d}_{cycle:3d}_{relative_orbit:3d}_{frame:4d}_{centre:3s}_{platform_mode:1s}_{timeliness:2s}_{collection:3s}.SEN3/tie_meteo.nc'
          - '{mission_id:3s}_OL_1_{datatype_id:_<6s}_{start_time:%Y%m%dT%H%M%S}_{end_time:%Y%m%dT%H%M%S}_{creation_time:%Y%m%dT%H%M%S}_{duration:4d}_{cycle:3d}_{relative_orbit:3d}______{centre:3s}_{platform_mode:1s}_{timeliness:2s}_{collection:3s}.SEN3/tie_meteo.nc'
+    esa_quality_flags:
+        file_reader: !!python/name:satpy.readers.olci_nc.NCOLCI1B
+        file_patterns:
+         - '{mission_id:3s}_OL_1_{datatype_id:_<6s}_{start_time:%Y%m%dT%H%M%S}_{end_time:%Y%m%dT%H%M%S}_{creation_time:%Y%m%dT%H%M%S}_{duration:4d}_{cycle:3d}_{relative_orbit:3d}_{frame:4d}_{centre:3s}_{platform_mode:1s}_{timeliness:2s}_{collection:3s}.SEN3/qualityFlags.nc'
+         - '{mission_id:3s}_OL_1_{datatype_id:_<6s}_{start_time:%Y%m%dT%H%M%S}_{end_time:%Y%m%dT%H%M%S}_{creation_time:%Y%m%dT%H%M%S}_{duration:4d}_{cycle:3d}_{relative_orbit:3d}______{centre:3s}_{platform_mode:1s}_{timeliness:2s}_{collection:3s}.SEN3/qualityFlags.nc'
+
+
 datasets:
   longitude:
     name: longitude
@@ -428,3 +435,18 @@ datasets:
     resolution: 300
     coordinates: [longitude, latitude]
     file_type: esa_meteo
+
+  quality_flags:
+    name: quality_flags
+    sensor: olci
+    resolution: 300
+    coordinates: [longitude, latitude]
+    file_type: esa_quality_flags
+
+  mask:
+    name: mask
+    sensor: olci
+    resolution: 300
+    coordinates: [longitude, latitude]
+    file_type: esa_quality_flags
+    nc_key: quality_flags

satpy/readers/olci_nc.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-# Copyright (c) 2016 Satpy developers
+# Copyright (c) 2016-2023 Satpy developers
 #
 # This file is part of satpy.
 #
@@ -51,9 +51,27 @@
 from satpy.readers.file_handlers import BaseFileHandler
 from satpy.utils import angle2xyz, get_legacy_chunk_size, xyz2angle
 
-DEFAULT_MASK_ITEMS = ["INVALID", "SNOW_ICE", "INLAND_WATER", "SUSPECT",
-                      "AC_FAIL", "CLOUD", "HISOLZEN", "OCNN_FAIL",
-                      "CLOUD_MARGIN", "CLOUD_AMBIGUOUS", "LOWRW", "LAND"]
+# the order of the L1B quality flags are from highest 32nd bit to the lowest 1 bit
+# https://sentinel.esa.int/documents/247904/1872756/Sentinel-3-OLCI-Product-Data-Format-Specification-OLCI-Level-1
+L1B_QUALITY_FLAGS = ["saturated@Oa21", "saturated@Oa20", "saturated@Oa19", "saturated@Oa18",
+                     "saturated@Oa17", "saturated@Oa16", "saturated@Oa15", "saturated@Oa14",
+                     "saturated@Oa13", "saturated@Oa12", "saturated@Oa11", "saturated@Oa10",
+                     "saturated@Oa09", "saturated@Oa08", "saturated@Oa07", "saturated@Oa06",
+                     "saturated@Oa05", "saturated@Oa04", "saturated@Oa03", "saturated@Oa02",
+                     "saturated@Oa01", "dubious", "sun-glint_risk", "duplicated",
+                     "cosmetic", "invalid", "straylight_risk", "bright",
+                     "tidal_region", "fresh_inland_water", "coastline", "land"]
+
+DEFAULT_L1B_MASK_ITEMS = ["dubious", "sun-glint_risk", "duplicated", "cosmetic", "invalid",
+    "straylight_risk", "bright", "tidal_region", "coastline", "land"]
+
+WQSF_FLAG_LIST = ["INVALID", "WATER", "LAND", "CLOUD", "SNOW_ICE", "INLAND_WATER", "TIDAL",
+    "COSMETIC", "SUSPECT", "HISOLZEN", "SATURATED", "MEGLINT", "HIGHGLINT", "WHITECAPS", "ADJAC",
+    "WV_FAIL", "PAR_FAIL", "AC_FAIL", "OC4ME_FAIL", "OCNN_FAIL", "Extra_1", "KDM_FAIL", "Extra_2",
+    "CLOUD_AMBIGUOUS", "CLOUD_MARGIN", "BPAC_ON", "WHITE_SCATT", "LOWRW", "HIGHRW"]
+
+DEFAULT_WQSF_MASK_ITEMS = ["INVALID", "SNOW_ICE", "INLAND_WATER", "SUSPECT", "AC_FAIL", "CLOUD",
+    "HISOLZEN", "OCNN_FAIL", "CLOUD_MARGIN", "CLOUD_AMBIGUOUS", "LOWRW", "LAND"]
 
 logger = logging.getLogger(__name__)
 
@@ -76,16 +94,7 @@ def __init__(self, value, flag_list=None):
                 meanings = value.attrs["flag_meanings"].split()
                 masks = value.attrs["flag_masks"]
             except (AttributeError, KeyError):
-                meanings = ["INVALID", "WATER", "LAND", "CLOUD", "SNOW_ICE",
-                            "INLAND_WATER", "TIDAL", "COSMETIC", "SUSPECT",
-                            "HISOLZEN", "SATURATED", "MEGLINT", "HIGHGLINT",
-                            "WHITECAPS", "ADJAC", "WV_FAIL", "PAR_FAIL",
-                            "AC_FAIL", "OC4ME_FAIL", "OCNN_FAIL",
-                            "Extra_1",
-                            "KDM_FAIL",
-                            "Extra_2",
-                            "CLOUD_AMBIGUOUS", "CLOUD_MARGIN", "BPAC_ON", "WHITE_SCATT",
-                            "LOWRW", "HIGHRW"]
+                meanings = WQSF_FLAG_LIST
                 self.meaning = {meaning: mask for mask, meaning in enumerate(meanings)}
             else:
                 self.meaning = {meaning: int(np.log2(mask)) for meaning, mask in zip(meanings, masks)}
@@ -99,8 +108,7 @@ def __getitem__(self, item):
         if isinstance(data, xr.DataArray):
             data = data.data
             res = ((data >> pos) % 2).astype(bool)
-            res = xr.DataArray(res, coords=self._value.coords,
-                               attrs=self._value.attrs,
+            res = xr.DataArray(res, coords=self._value.coords, attrs=self._value.attrs,
                                dims=self._value.dims)
         else:
             res = ((data >> pos) % 2).astype(bool)
@@ -113,8 +121,7 @@ class NCOLCIBase(BaseFileHandler):
     rows_name = "rows"
     cols_name = "columns"
 
-    def __init__(self, filename, filename_info, filetype_info,
-                 engine=None, **kwargs):
+    def __init__(self, filename, filename_info, filetype_info, engine=None, **kwargs):
         """Init the olci reader base."""
         super().__init__(filename, filename_info, filetype_info)
         self._engine = engine
@@ -176,10 +183,12 @@ def __init__(self, filename, filename_info, filetype_info, engine=None):
 class NCOLCI1B(NCOLCIChannelBase):
     """File handler for OLCI l1b."""
 
-    def __init__(self, filename, filename_info, filetype_info, cal, engine=None):
+    def __init__(self, filename, filename_info, filetype_info, cal=None, engine=None, mask_items=None):
         """Init the file handler."""
         super().__init__(filename, filename_info, filetype_info, engine)
-        self.cal = cal.nc
+        if cal is not None:
+            self.cal = cal.nc
+        self.mask_items = mask_items
 
     @staticmethod
     def _get_items(idx, solar_flux):
@@ -194,24 +203,40 @@ def _get_solar_flux(self, band):
         return da.map_blocks(self._get_items, d_index.data,
                              solar_flux=solar_flux, dtype=solar_flux.dtype)
 
+    def getbitmask(self, quality_flags, items=None):
+        """Get the quality flags bitmask."""
+        if items is None:
+            items = DEFAULT_L1B_MASK_ITEMS
+
+        bflags = BitFlags(quality_flags, flag_list=L1B_QUALITY_FLAGS)
+
+        return reduce(np.logical_or, [bflags[item] for item in items])
+
     def get_dataset(self, key, info):
         """Load a dataset."""
-        if self.channel != key["name"]:
+        if self.channel is not None and self.channel != key["name"]:
             return
+
         logger.debug("Reading %s.", key["name"])
 
-        radiances = self.nc[self.channel + "_radiance"]
+        if key["name"] == "quality_flags":
+            dataset = self.nc["quality_flags"]
+        elif key["name"] == "mask":
+            dataset = self.getbitmask(self.nc["quality_flags"], self.mask_items)
+        else:
+            dataset = self.nc[self.channel + "_radiance"]
 
-        if key["calibration"] == "reflectance":
-            idx = int(key["name"][2:]) - 1
-            sflux = self._get_solar_flux(idx)
-            radiances = radiances / sflux * np.pi * 100
-            radiances.attrs["units"] = "%"
+            if key["calibration"] == "reflectance":
+                idx = int(key["name"][2:]) - 1
+                sflux = self._get_solar_flux(idx)
+                dataset = dataset / sflux * np.pi * 100
+                dataset.attrs["units"] = "%"
 
-        radiances.attrs["platform_name"] = self.platform_name
-        radiances.attrs["sensor"] = self.sensor
-        radiances.attrs.update(key.to_dict())
-        return radiances
+        dataset.attrs["platform_name"] = self.platform_name
+        dataset.attrs["sensor"] = self.sensor
+        dataset.attrs.update(key.to_dict())
+
+        return dataset
 
 
 class NCOLCI2(NCOLCIChannelBase):
@@ -228,6 +253,7 @@ def get_dataset(self, key, info):
         if self.channel is not None and self.channel != key["name"]:
             return
         logger.debug("Reading %s.", key["name"])
+
         if self.channel is not None and self.channel.startswith(self.reflectance_prefix):
             dataset = self.nc[self.channel + self.reflectance_suffix]
         else:
@@ -237,6 +263,7 @@ def get_dataset(self, key, info):
             dataset.attrs["_FillValue"] = 1
         elif key["name"] == "mask":
             dataset = self.getbitmask(dataset, self.mask_items)
+
         dataset.attrs["platform_name"] = self.platform_name
         dataset.attrs["sensor"] = self.sensor
         dataset.attrs.update(key.to_dict())
@@ -257,8 +284,8 @@ def delog(self, data_array):
     def getbitmask(self, wqsf, items=None):
         """Get the bitmask."""
         if items is None:
-            items = DEFAULT_MASK_ITEMS
-        bflags = BitFlags(wqsf)
+            items = DEFAULT_WQSF_MASK_ITEMS
+        bflags = BitFlags(wqsf, WQSF_FLAG_LIST)
         return reduce(np.logical_or, [bflags[item] for item in items])
 
 

satpy/tests/reader_tests/test_olci_nc.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-# Copyright (c) 2016-2018 Satpy developers
+# Copyright (c) 2016-2023 Satpy developers
 #
 # This file is part of satpy.
 #
@@ -94,7 +94,7 @@ def test_open_file_objects(self, mocked_open_dataset):
                 open_file.open.return_value == mocked_open_dataset.call_args[1].get("filename_or_obj"))
 
     @mock.patch("xarray.open_dataset")
-    def test_get_mask(self, mocked_dataset):
+    def test_get_l2_mask(self, mocked_dataset):
         """Test reading datasets."""
         import numpy as np
         import xarray as xr
@@ -118,7 +118,7 @@ def test_get_mask(self, mocked_dataset):
         np.testing.assert_array_equal(res.values, expected)
 
     @mock.patch("xarray.open_dataset")
-    def test_get_mask_with_alternative_items(self, mocked_dataset):
+    def test_get_l2_mask_with_alternative_items(self, mocked_dataset):
         """Test reading datasets."""
         import numpy as np
         import xarray as xr
@@ -137,6 +137,61 @@ def test_get_mask_with_alternative_items(self, mocked_dataset):
         expected = np.array([True] + [False] * 29).reshape(5, 6)
         np.testing.assert_array_equal(res.values, expected)
 
+
+    @mock.patch("xarray.open_dataset")
+    def test_get_l1b_default_mask(self, mocked_dataset):
+        """Test reading mask datasets from L1B products."""
+        import numpy as np
+        import xarray as xr
+
+        from satpy.readers.olci_nc import NCOLCI1B
+        from satpy.tests.utils import make_dataid
+        mocked_dataset.return_value = xr.Dataset({"quality_flags": (["rows", "columns"],
+                                                           np.array([1 << (x % 32) for x in range(35)]).reshape(5, 7))},
+                                                 coords={"rows": np.arange(5),
+                                                         "columns": np.arange(7)})
+        ds_id = make_dataid(name="mask")
+        filename_info = {"mission_id": "S3A", "dataset_name": "mask", "start_time": 0, "end_time": 0}
+        test = NCOLCI1B("somedir/somefile.nc", filename_info, "c")
+        res = test.get_dataset(ds_id, {"nc_key": "quality_flags"})
+        assert res.dtype == np.dtype("bool")
+
+        expected = np.array([[False, False, False, False, False, False, False],
+                             [False, False, False, False, False, False, False],
+                             [False, False, False, False, False, False, False],
+                             [True, True, True, True, True, True, True],
+                             [True, False, True, True, False, False, False]])
+
+        np.testing.assert_array_equal(res.values, expected)
+
+
+    @mock.patch("xarray.open_dataset")
+    def test_get_l1b_customized_mask(self, mocked_dataset):
+        """Test reading mask datasets from L1B products."""
+        import numpy as np
+        import xarray as xr
+
+        from satpy.readers.olci_nc import NCOLCI1B
+        from satpy.tests.utils import make_dataid
+        mocked_dataset.return_value = xr.Dataset({"quality_flags": (["rows", "columns"],
+                                                           np.array([1 << (x % 32) for x in range(35)]).reshape(5, 7))},
+                                                 coords={"rows": np.arange(5),
+                                                         "columns": np.arange(7)})
+        ds_id = make_dataid(name="mask")
+        filename_info = {"mission_id": "S3A", "dataset_name": "mask", "start_time": 0, "end_time": 0}
+        test = NCOLCI1B("somedir/somefile.nc", filename_info, "c", mask_items=["bright", "invalid"])
+        res = test.get_dataset(ds_id, {"nc_key": "quality_flags"})
+        assert res.dtype == np.dtype("bool")
+
+        expected = np.array([[False, False, False, False, False, False, False],
+                             [False, False, False, False, False, False, False],
+                             [False, False, False, False, False, False, False],
+                             [False, False, False, False, True, False, True],
+                             [False, False, False, False, False, False, False]])
+
+        np.testing.assert_array_equal(res.values, expected)
+
+
     @mock.patch("xarray.open_dataset")
     def test_olci_angles(self, mocked_dataset):
         """Test reading datasets."""
@@ -241,7 +296,7 @@ def test_chl_nn(self, mocked_dataset):
         assert res.values[-1, -1] == 1e29
 
 
-class TestBitFlags(unittest.TestCase):
+class TestL2BitFlags(unittest.TestCase):
     """Test the bitflag reading."""
 
     def test_bitflags(self):
@@ -373,3 +428,45 @@ def test_bitflags_with_custom_flag_list(self):
                              False, False, True, True, False, False, True,
                              False])
         assert all(mask == expected)
+
+
+class TestL1bBitFlags(unittest.TestCase):
+    """Test the bitflag reading."""
+
+    def test_bitflags(self):
+        """Test the BitFlags class."""
+        from functools import reduce
+
+        import numpy as np
+
+        from satpy.readers.olci_nc import BitFlags
+
+
+        L1B_QUALITY_FLAGS = ["saturated@Oa21", "saturated@Oa20", "saturated@Oa19", "saturated@Oa18",
+                             "saturated@Oa17", "saturated@Oa16", "saturated@Oa15", "saturated@Oa14",
+                             "saturated@Oa13", "saturated@Oa12", "saturated@Oa11", "saturated@Oa10",
+                             "saturated@Oa09", "saturated@Oa08", "saturated@Oa07", "saturated@Oa06",
+                             "saturated@Oa05", "saturated@Oa04", "saturated@Oa03", "saturated@Oa02",
+                             "saturated@Oa01", "dubious", "sun-glint_risk", "duplicated",
+                             "cosmetic", "invalid", "straylight_risk", "bright",
+                             "tidal_region", "fresh_inland_water", "coastline", "land"]
+
+        DEFAULT_L1B_MASK_ITEMS = ["dubious", "sun-glint_risk", "duplicated", "cosmetic", "invalid",
+                                  "straylight_risk", "bright", "tidal_region", "coastline", "land"]
+
+        bits = np.array([1 << x for x in range(len(L1B_QUALITY_FLAGS))])
+
+        bflags = BitFlags(bits, flag_list=L1B_QUALITY_FLAGS)
+
+        mask = reduce(np.logical_or, [bflags[item] for item in DEFAULT_L1B_MASK_ITEMS])
+
+        expected = np.array([False, False, False, False,
+                             False, False, False, False,
+                             False, False, False, False,
+                             False, False, False, False,
+                             False, False, False, False,
+                             False, True, True, True,
+                             True, True, True, True,
+                             True, False, True, True,
+                            ])
+        assert all(mask == expected)