Merge pull request #2714 from djhoese/feat-viirs-aod-qc-filter

Add QC-based filtering to the VIIRS EDR AOD550 product

satpy/etc/readers/viirs_edr.yaml

@@ -25,7 +25,7 @@ file_types:
         file_patterns:
             - 'JRR-CloudHeight_{version}_{platform_shortname}_s{start_time:%Y%m%d%H%M%S%f}_e{end_time:%Y%m%d%H%M%S%f}_c{creation_time}.nc'
     jrr_aod:
-        file_reader: !!python/name:satpy.readers.viirs_edr.VIIRSJRRFileHandler
+        file_reader: !!python/name:satpy.readers.viirs_edr.VIIRSAODHandler
         file_patterns:
             - 'JRR-AOD_{version}_{platform_shortname}_s{start_time:%Y%m%d%H%M%S%f}_e{end_time:%Y%m%d%H%M%S%f}_c{creation_time}.nc'
     jrr_lst:

satpy/readers/viirs_edr.py

@@ -56,6 +56,21 @@
 
     scene = satpy.Scene(filenames, reader='viirs_edr', reader_kwargs={"filter_veg": False})
 
+AOD Filtering
+^^^^^^^^^^^^^
+
+The AOD (Aerosol Optical Depth) product can be optionally filtered based on
+Quality Control (QC) values in the file. By default no filtering is performed.
+By providing the ``aod_qc_filter`` keyword argument and specifying the maximum
+value of the ``QCAll`` variable to include (not mask). For example::
+
+    scene = satpy.Scene(filenames, reader='viirs_edr', reader_kwargs={"aod_qc_filter": 1})
+
+will only preserve AOD550 values where the quality is 0 ("high") or
+1 ("medium"). At the time of writing the ``QCAll`` variable has 1 ("medium"),
+2 ("low"), and 3 ("no retrieval").
+
+
 """
 from __future__ import annotations
 
@@ -76,7 +91,7 @@
 class VIIRSJRRFileHandler(BaseFileHandler):
     """NetCDF4 reader for VIIRS Active Fires."""
 
-    def __init__(self, filename, filename_info, filetype_info):
+    def __init__(self, filename, filename_info, filetype_info, **kwargs):
         """Initialize the geo filehandler."""
         super(VIIRSJRRFileHandler, self).__init__(filename, filename_info,
                                                   filetype_info)
@@ -343,3 +358,20 @@ def _scale_data(self):
             add_offset = self.nc[self._manual_scalings[var_name][1]]
             data_arr.data = data_arr.data * scale_factor.data + add_offset.data
             self.nc[var_name] = data_arr
+
+
+class VIIRSAODHandler(VIIRSJRRFileHandler):
+    """File handler for AOD data files."""
+
+    def __init__(self, *args, aod_qc_filter: int | None = None, **kwargs) -> None:
+        """Initialize file handler and keep track of QC filtering."""
+        super().__init__(*args, **kwargs)
+        self._aod_qc_filter = aod_qc_filter
+
+    def _mask_invalid(self, data_arr: xr.DataArray, ds_info: dict) -> xr.DataArray:
+        new_data_arr = super()._mask_invalid(data_arr, ds_info)
+        if self._aod_qc_filter is None or ds_info["name"] != "AOD550":
+            return new_data_arr
+        LOG.debug(f"Filtering AOD data to include quality <= {self._aod_qc_filter}")
+        qc_all = self.nc["QCAll"]
+        return new_data_arr.where(qc_all <= self._aod_qc_filter)

satpy/tests/reader_tests/test_viirs_edr.py

@@ -220,6 +220,14 @@ def aod_file(tmp_path_factory: TempPathFactory) -> Path:
     data_vars = _create_continuous_variables(
         ("AOD550",)
     )
+    qc_data = np.zeros(data_vars["AOD550"].shape, dtype=np.int8)
+    qc_data[-1, -1] = 2
+    data_vars["QCAll"] = xr.DataArray(
+        qc_data,
+        dims=data_vars["AOD550"].dims,
+        attrs={"valid_range": [0, 3]},
+    )
+    data_vars["QCAll"].encoding["_FillValue"] = -128
     return _create_fake_file(tmp_path_factory, fn, data_vars)
 
 
@@ -371,7 +379,6 @@ def test_get_dataset_surf_refl_with_veg_idx(
         ("var_names", "data_file"),
         [
             (("CldTopTemp", "CldTopHght", "CldTopPres"), lazy_fixture("cloud_height_file")),
-            (("AOD550",), lazy_fixture("aod_file")),
             (("VLST",), lazy_fixture("lst_file")),
         ]
     )
@@ -385,6 +392,31 @@ def test_get_dataset_generic(self, var_names, data_file):
         for var_name in var_names:
             _check_continuous_data_arr(scn[var_name])
 
+    @pytest.mark.parametrize(
+        ("aod_qc_filter", "exp_masked_pixel"),
+        [
+            (None, False),
+            (0, True),
+            (2, False)
+        ],
+    )
+    def test_get_aod_filtered(self, aod_file, aod_qc_filter, exp_masked_pixel):
+        """Test that the AOD product can be loaded and filtered."""
+        from satpy import Scene
+        bytes_in_m_row = 4 * 3200
+        with dask.config.set({"array.chunk-size": f"{bytes_in_m_row * 4}B"}):
+            scn = Scene(reader="viirs_edr", filenames=[aod_file], reader_kwargs={"aod_qc_filter": aod_qc_filter})
+            scn.load(["AOD550"])
+        _check_continuous_data_arr(scn["AOD550"])
+        data_np = scn["AOD550"].data.compute()
+        pixel_is_nan = np.isnan(data_np[-1, -1])
+        assert pixel_is_nan if exp_masked_pixel else not pixel_is_nan
+
+        # filtering should never affect geolocation
+        lons, lats = scn["AOD550"].attrs["area"].get_lonlats()
+        assert not np.isnan(lons[-1, -1].compute())
+        assert not np.isnan(lats[-1, -1].compute())
+
     @pytest.mark.parametrize(
         ("data_file", "exp_available"),
         [