Merge pull request #1993 from djhoese/feature-seadas-l2-flags

Add 'l2_flags' quality filtering to 'seadas_l2' reader

satpy/readers/seadas_l2.py

@@ -19,6 +19,13 @@
 
 This reader currently only supports MODIS and VIIRS Chlorophyll A from SEADAS.
 
+The reader includes an additional keyword argument ``apply_quality_flags``
+which can be used to mask out low-quality pixels based on quality flags
+contained in the file (``l2_flags``). This option defaults to ``False``, but
+when set to ``True`` the "CHLWARN" pixels of the ``l2_flags`` variable
+are masked out. These pixels represent data where the chlorophyll algorithm
+warned about the quality of the result.
+
 """
 
 from datetime import datetime
@@ -31,6 +38,11 @@
 class SEADASL2HDFFileHandler(HDF4FileHandler):
     """Simple handler of SEADAS L2 files."""
 
+    def __init__(self, filename, filename_info, filetype_info, apply_quality_flags=False):
+        """Initialize file handler and determine if data quality flags should be applied."""
+        super().__init__(filename, filename_info, filetype_info)
+        self.apply_quality_flags = apply_quality_flags and "l2_flags" in self
+
     def _add_satpy_metadata(self, data):
         data.attrs["sensor"] = self.sensor_names
         data.attrs["platform_name"] = self._platform_name()
@@ -79,6 +91,10 @@ def get_dataset(self, data_id, dataset_info):
         valid_range = data.attrs["valid_range"]
         data = data.where(valid_range[0] <= data)
         data = data.where(data <= valid_range[1])
+        if self.apply_quality_flags and not ("lon" in file_key or "lat" in file_key):
+            l2_flags = self["l2_flags"]
+            mask = (l2_flags & 0b00000000010000000000000000000000) != 0
+            data = data.where(~mask)
         for attr_name in ("standard_name", "long_name", "units"):
             val = data.attrs[attr_name]
             if val[-1] == "\x00":

satpy/tests/reader_tests/test_seadas_l2.py

@@ -84,7 +84,7 @@ def _create_seadas_chlor_a_file(full_path, mission, sensor):
 
     chlor_a_info = {
         "type": SDC.FLOAT32,
-        "data": np.zeros((5, 5), np.float32),
+        "data": np.ones((5, 5), np.float32),
         "dim_labels": ["Number of Scan Lines", "Number of Pixel Control Points"],
         "attrs": {
             "long_name": "Chlorophyll Concentration, OCI Algorithm\x00",
@@ -94,6 +94,16 @@ def _create_seadas_chlor_a_file(full_path, mission, sensor):
         }
     }
     _add_variable_to_file(h, "chlor_a", chlor_a_info)
+
+    l2_flags = np.zeros((5, 5), dtype=np.int32)
+    l2_flags[2, 2] = -1
+    l2_flags_info = {
+        "type": SDC.INT32,
+        "data": l2_flags,
+        "dim_labels": ["Number of Scan Lines", "Number of Pixel Control Points"],
+        "attrs": {},
+    }
+    _add_variable_to_file(h, "l2_flags", l2_flags_info)
     return [full_path]
 
 
@@ -136,9 +146,11 @@ def test_scene_available_datasets(self, input_files):
             (lazy_fixture("seadas_l2_viirs_npp_chlor_a"), "Suomi-NPP", {"viirs"}, 16),
             (lazy_fixture("seadas_l2_viirs_j01_chlor_a"), "NOAA-20", {"viirs"}, 16),
         ])
-    def test_load_chlor_a(self, input_files, exp_plat, exp_sensor, exp_rps):
+    @pytest.mark.parametrize("apply_quality_flags", [False, True])
+    def test_load_chlor_a(self, input_files, exp_plat, exp_sensor, exp_rps, apply_quality_flags):
         """Test that we can load 'chlor_a'."""
-        scene = Scene(reader='seadas_l2', filenames=input_files)
+        reader_kwargs = {"apply_quality_flags": apply_quality_flags}
+        scene = Scene(reader='seadas_l2', filenames=input_files, reader_kwargs=reader_kwargs)
         scene.load(['chlor_a'])
         data_arr = scene['chlor_a']
         assert data_arr.attrs['platform_name'] == exp_plat
@@ -147,3 +159,9 @@ def test_load_chlor_a(self, input_files, exp_plat, exp_sensor, exp_rps):
         assert data_arr.dtype.type == np.float32
         assert isinstance(data_arr.attrs["area"], SwathDefinition)
         assert data_arr.attrs["rows_per_scan"] == exp_rps
+        data = data_arr.data.compute()
+        if apply_quality_flags:
+            assert np.isnan(data[2, 2])
+            assert np.count_nonzero(np.isnan(data)) == 1
+        else:
+            assert np.count_nonzero(np.isnan(data)) == 0