Merge pull request #2209 from djhoese/feature-seadas-nc

Update seadas_l2 reader to handle alternative NetCDF file format

satpy/etc/readers/seadas_l2.yaml

@@ -15,41 +15,53 @@ file_types:
     - '{platform_indicator:1s}1.{start_time:%y%j.%H%M}.seadas.hdf'
     file_reader: !!python/name:satpy.readers.seadas_l2.SEADASL2HDFFileHandler
     geo_resolution: 1000
+  chlora_seadas_nc:
+    file_patterns:
+      # IMAPP-style filenames:
+      - '{platform_indicator:1s}1.{start_time:%y%j.%H%M}.seadas.nc'
+    file_reader: !!python/name:satpy.readers.seadas_l2.SEADASL2NetCDFFileHandler
+    geo_resolution: 1000
   chlora_seadas_viirs:
     # SEADAS_npp_d20211118_t1728125_e1739327.hdf
     file_patterns:
       - 'SEADAS_{platform_indicator:s}_d{start_time:%Y%m%d_t%H%M%S%f}_e{end_time:%H%M%S%f}.hdf'
     file_reader: !!python/name:satpy.readers.seadas_l2.SEADASL2HDFFileHandler
     geo_resolution: 750
+  chlora_seadas_viirs_nc:
+    # SEADAS_npp_d20211118_t1728125_e1739327.nc
+    file_patterns:
+      - 'SEADAS_{platform_indicator:s}_d{start_time:%Y%m%d_t%H%M%S%f}_e{end_time:%H%M%S%f}.nc'
+    file_reader: !!python/name:satpy.readers.seadas_l2.SEADASL2NetCDFFileHandler
+    geo_resolution: 750
 
 datasets:
   longitude:
     name: longitude
-    file_type: [chlora_seadas, seadas_hdf_viirs]
-    file_key: longitude
+    file_type: [chlora_seadas_viirs_nc, chlora_seadas_nc, seadas_hdf_viirs, chlora_sedas]
+    file_key: ["navigation_data/longitude", "longitude"]
     resolution:
       1000:
-        file_type: chlora_seadas
+        file_type: [chlora_seadas_nc, chlora_seadas]
       750:
-        file_type: chlora_seadas_viirs
+        file_type: [chlora_seadas_viirs_nc, chlora_seadas_viirs]
 
   latitude:
     name: latitude
-    file_type: [chlora_seadas, seadas_hdf_viirs]
-    file_key: latitude
+    file_type: [chlora_seadas_viirs_nc, chlora_seadas_nc, seadas_hdf_viirs, chlora_sedas]
+    file_key: ["navigation_data/latitude", "latitude"]
     resolution:
       1000:
-        file_type: chlora_seadas
+        file_type: [chlora_seadas_nc, chlora_seadas]
       750:
-        file_type: chlora_seadas_viirs
+        file_type: [chlora_seadas_viirs_nc, chlora_seadas_viirs]
 
   chlor_a:
     name: chlor_a
-    file_type: [chlora_seadas, seadas_hdf_viirs]
-    file_key: chlor_a
+    file_type: [chlora_seadas_viirs_nc, chlora_seadas_nc, seadas_hdf_viirs, chlora_sedas]
+    file_key: ["geophysical_data/chlor_a", "chlor_a"]
     resolution:
       1000:
-        file_type: chlora_seadas
+        file_type: [chlora_seadas_nc, chlora_seadas]
       750:
-        file_type: chlora_seadas_viirs
+        file_type: [chlora_seadas_viirs_nc, chlora_seadas_viirs]
     coordinates: [longitude, latitude]

satpy/readers/seadas_l2.py

@@ -31,17 +31,16 @@
 from datetime import datetime
 
 from .hdf4_utils import HDF4FileHandler
+from .netcdf_utils import NetCDF4FileHandler
 
-TIME_FORMAT = "%Y%j%H%M%S"
 
-
-class SEADASL2HDFFileHandler(HDF4FileHandler):
+class _SEADASL2Base:
     """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
+        self.apply_quality_flags = apply_quality_flags and self.l2_flags_var_name in self
 
     def _add_satpy_metadata(self, data):
         data.attrs["sensor"] = self.sensor_names
@@ -57,7 +56,7 @@ def _rows_per_scan(self):
         raise ValueError(f"Don't know how to read data for sensors: {self.sensor_names}")
 
     def _platform_name(self):
-        platform = self["/attr/Mission"]
+        platform = self[self.platform_attr_name]
         platform_dict = {'NPP': 'Suomi-NPP',
                          'JPSS-1': 'NOAA-20',
                          'JPSS-2': 'NOAA-21'}
@@ -66,38 +65,93 @@ def _platform_name(self):
     @property
     def start_time(self):
         """Get the starting observation time of this file's data."""
-        start_time = self["/attr/Start Time"]
-        return datetime.strptime(start_time[:-3], TIME_FORMAT)
+        start_time = self[self.start_time_attr_name]
+        return datetime.strptime(start_time[:-3], self.time_format)
 
     @property
     def end_time(self):
         """Get the ending observation time of this file's data."""
-        end_time = self["/attr/End Time"]
-        return datetime.strptime(end_time[:-3], TIME_FORMAT)
+        end_time = self[self.end_time_attr_name]
+        return datetime.strptime(end_time[:-3], self.time_format)
 
     @property
     def sensor_names(self):
         """Get sensor for the current file's data."""
         # Example: MODISA or VIIRSN or VIIRSJ1
-        sensor_name = self["/attr/Sensor Name"].lower()
+        sensor_name = self[self.sensor_attr_name].lower()
         if sensor_name.startswith("modis"):
             return {"modis"}
         return {"viirs"}
 
     def get_dataset(self, data_id, dataset_info):
         """Get DataArray for the specified DataID."""
-        file_key = dataset_info.get("file_key", data_id["name"])
-        data = self[file_key]
-        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)
+        file_key, data = self._get_file_key_and_variable(data_id, dataset_info)
+        data = self._filter_by_valid_min_max(data)
+        data = self._rename_2d_dims_if_necessary(data)
+        data = self._mask_based_on_l2_flags(data)
         for attr_name in ("standard_name", "long_name", "units"):
             val = data.attrs[attr_name]
             if val[-1] == "\x00":
                 data.attrs[attr_name] = data.attrs[attr_name][:-1]
         data = self._add_satpy_metadata(data)
         return data
+
+    def _get_file_key_and_variable(self, data_id, dataset_info):
+        file_keys = dataset_info.get("file_key", data_id["name"])
+        if not isinstance(file_keys, list):
+            file_keys = [file_keys]
+        for file_key in file_keys:
+            try:
+                data = self[file_key]
+                return file_key, data
+            except KeyError:
+                continue
+        raise KeyError(f"Unable to find any of the possible keys for {data_id}: {file_keys}")
+
+    def _rename_2d_dims_if_necessary(self, data_arr):
+        if data_arr.ndim != 2 or data_arr.dims == ("y", "x"):
+            return data_arr
+        return data_arr.rename(dict(zip(data_arr.dims, ("y", "x"))))
+
+    def _filter_by_valid_min_max(self, data_arr):
+        valid_range = self._valid_min_max(data_arr)
+        data_arr = data_arr.where(valid_range[0] <= data_arr)
+        data_arr = data_arr.where(data_arr <= valid_range[1])
+        return data_arr
+
+    def _valid_min_max(self, data_arr):
+        try:
+            return data_arr.attrs["valid_range"]
+        except KeyError:
+            return data_arr.attrs["valid_min"], data_arr.attrs["valid_max"]
+
+    def _mask_based_on_l2_flags(self, data_arr):
+        standard_name = data_arr.attrs.get("standard_name", "")
+        if self.apply_quality_flags and not ("lon" in standard_name or "lat" in standard_name):
+            l2_flags = self[self.l2_flags_var_name]
+            l2_flags = self._rename_2d_dims_if_necessary(l2_flags)
+            mask = (l2_flags & 0b00000000010000000000000000000000) != 0
+            data_arr = data_arr.where(~mask)
+        return data_arr
+
+
+class SEADASL2NetCDFFileHandler(_SEADASL2Base, NetCDF4FileHandler):
+    """Simple handler of SEADAS L2 NetCDF4 files."""
+
+    start_time_attr_name = "/attr/time_coverage_start"
+    end_time_attr_name = "/attr/time_coverage_end"
+    time_format = "%Y-%m-%dT%H:%M:%S.%f"
+    platform_attr_name = "/attr/platform"
+    sensor_attr_name = "/attr/instrument"
+    l2_flags_var_name = "geophysical_data/l2_flags"
+
+
+class SEADASL2HDFFileHandler(_SEADASL2Base, HDF4FileHandler):
+    """Simple handler of SEADAS L2 HDF4 files."""
+
+    start_time_attr_name = "/attr/Start Time"
+    end_time_attr_name = "/attr/End Time"
+    time_format = "%Y%j%H%M%S"
+    platform_attr_name = "/attr/Mission"
+    sensor_attr_name = "/attr/Sensor Name"
+    l2_flags_var_name = "l2_flags"

satpy/tests/reader_tests/test_seadas_l2.py

@@ -19,7 +19,6 @@
 
 import numpy as np
 import pytest
-from pyhdf.SD import SD, SDC
 from pyresample.geometry import SwathDefinition
 from pytest_lazyfixture import lazy_fixture
 
@@ -31,26 +30,27 @@ def seadas_l2_modis_chlor_a(tmp_path_factory):
     """Create MODIS SEADAS file."""
     filename = "a1.21322.1758.seadas.hdf"
     full_path = str(tmp_path_factory.mktemp("seadas_l2") / filename)
-    return _create_seadas_chlor_a_file(full_path, "Aqua", "MODISA")
+    return _create_seadas_chlor_a_hdf4_file(full_path, "Aqua", "MODISA")
 
 
 @pytest.fixture(scope="module")
 def seadas_l2_viirs_npp_chlor_a(tmp_path_factory):
     """Create VIIRS NPP SEADAS file."""
     filename = "SEADAS_npp_d20211118_t1728125_e1739327.hdf"
     full_path = str(tmp_path_factory.mktemp("seadas") / filename)
-    return _create_seadas_chlor_a_file(full_path, "NPP", "VIIRSN")
+    return _create_seadas_chlor_a_hdf4_file(full_path, "NPP", "VIIRSN")
 
 
 @pytest.fixture(scope="module")
 def seadas_l2_viirs_j01_chlor_a(tmp_path_factory):
     """Create VIIRS JPSS-01 SEADAS file."""
     filename = "SEADAS_j01_d20211118_t1728125_e1739327.hdf"
     full_path = str(tmp_path_factory.mktemp("seadas") / filename)
-    return _create_seadas_chlor_a_file(full_path, "JPSS-1", "VIIRSJ1")
+    return _create_seadas_chlor_a_hdf4_file(full_path, "JPSS-1", "VIIRSJ1")
 
 
-def _create_seadas_chlor_a_file(full_path, mission, sensor):
+def _create_seadas_chlor_a_hdf4_file(full_path, mission, sensor):
+    from pyhdf.SD import SD, SDC
     h = SD(full_path, SDC.WRITE | SDC.CREATE)
     setattr(h, "Sensor Name", sensor)
     h.Mission = mission
@@ -79,8 +79,8 @@ def _create_seadas_chlor_a_file(full_path, mission, sensor):
             "valid_range": (-90.0, 90.0),
         }
     }
-    _add_variable_to_file(h, "longitude", lon_info)
-    _add_variable_to_file(h, "latitude", lat_info)
+    _add_variable_to_hdf4_file(h, "longitude", lon_info)
+    _add_variable_to_hdf4_file(h, "latitude", lat_info)
 
     chlor_a_info = {
         "type": SDC.FLOAT32,
@@ -93,7 +93,7 @@ def _create_seadas_chlor_a_file(full_path, mission, sensor):
             "valid_range": (0.001, 100.0),
         }
     }
-    _add_variable_to_file(h, "chlor_a", chlor_a_info)
+    _add_variable_to_hdf4_file(h, "chlor_a", chlor_a_info)
 
     l2_flags = np.zeros((5, 5), dtype=np.int32)
     l2_flags[2, 2] = -1
@@ -103,11 +103,11 @@ def _create_seadas_chlor_a_file(full_path, mission, sensor):
         "dim_labels": ["Number of Scan Lines", "Number of Pixel Control Points"],
         "attrs": {},
     }
-    _add_variable_to_file(h, "l2_flags", l2_flags_info)
+    _add_variable_to_hdf4_file(h, "l2_flags", l2_flags_info)
     return [full_path]
 
 
-def _add_variable_to_file(h, var_name, var_info):
+def _add_variable_to_hdf4_file(h, var_name, var_info):
     v = h.create(var_name, var_info['type'], var_info['data'].shape)
     v[:] = var_info['data']
     for dim_count, dimension_name in enumerate(var_info['dim_labels']):
@@ -118,6 +118,85 @@ def _add_variable_to_file(h, var_name, var_info):
         setattr(v, attr_key, attr_val)
 
 
+@pytest.fixture(scope="module")
+def seadas_l2_modis_chlor_a_netcdf(tmp_path_factory):
+    """Create MODIS SEADAS NetCDF file."""
+    filename = "t1.21332.1758.seadas.nc"
+    full_path = str(tmp_path_factory.mktemp("seadas_l2") / filename)
+    return _create_seadas_chlor_a_netcdf_file(full_path, "Terra", "MODIS")
+
+
+def _create_seadas_chlor_a_netcdf_file(full_path, mission, sensor):
+    from netCDF4 import Dataset
+    nc = Dataset(full_path, "w")
+    nc.createDimension("number_of_lines", 5)
+    nc.createDimension("pixels_per_line", 5)
+    nc.instrument = sensor
+    nc.platform = mission
+    nc.time_coverage_start = "2021-11-18T17:58:53.191Z"
+    nc.time_coverage_end = "2021-11-18T18:05:51.214Z"
+
+    lon_info = {
+        "data": np.zeros((5, 5), dtype=np.float32),
+        "dim_labels": ("number_of_lines", "pixels_per_line"),
+        "attrs": {
+            "long_name": "Longitude",
+            "standard_name": "longitude",
+            "units": "degrees_east",
+            "valid_min": -180.0,
+            "valid_max": 180.0,
+        }
+    }
+    lat_info = {
+        "data": np.zeros((5, 5), np.float32),
+        "dim_labels": ("number_of_lines", "pixels_per_line"),
+        "attrs": {
+            "long_name": "Latitude",
+            "standard_name": "latitude",
+            "units": "degrees_north",
+            "valid_min": -90.0,
+            "valid_max": 90.0,
+        }
+    }
+    nav_group = nc.createGroup("navigation_data")
+    _add_variable_to_netcdf_file(nav_group, "longitude", lon_info)
+    _add_variable_to_netcdf_file(nav_group, "latitude", lat_info)
+
+    chlor_a_info = {
+        "data": np.ones((5, 5), np.float32),
+        "dim_labels": ("number_of_lines", "pixels_per_line"),
+        "attrs": {
+            "long_name": "Chlorophyll Concentration, OCI Algorithm",
+            "units": "mg m^-3",
+            "standard_name": "mass_concentration_of_chlorophyll_in_sea_water",
+            "valid_min": 0.001,
+            "valid_max": 100.0,
+        }
+    }
+    l2_flags = np.zeros((5, 5), dtype=np.int32)
+    l2_flags[2, 2] = -1
+    l2_flags_info = {
+        "data": l2_flags,
+        "dim_labels": ("number_of_lines", "pixels_per_line"),
+        "attrs": {
+            "valid_min": -2147483648,
+            "valid_max": 2147483647,
+        },
+    }
+    geophys_group = nc.createGroup("geophysical_data")
+    _add_variable_to_netcdf_file(geophys_group, "chlor_a", chlor_a_info)
+    _add_variable_to_netcdf_file(geophys_group, "l2_flags", l2_flags_info)
+    return [full_path]
+
+
+def _add_variable_to_netcdf_file(nc, var_name, var_info):
+    v = nc.createVariable(var_name, var_info["data"].dtype.str[1:], dimensions=var_info["dim_labels"],
+                          fill_value=var_info.get("fill_value"))
+    v[:] = var_info['data']
+    for attr_key, attr_val in var_info['attrs'].items():
+        setattr(v, attr_key, attr_val)
+
+
 class TestSEADAS:
     """Test the SEADAS L2 file reader."""
 
@@ -145,6 +224,7 @@ def test_scene_available_datasets(self, input_files):
             (lazy_fixture("seadas_l2_modis_chlor_a"), "Aqua", {"modis"}, 10),
             (lazy_fixture("seadas_l2_viirs_npp_chlor_a"), "Suomi-NPP", {"viirs"}, 16),
             (lazy_fixture("seadas_l2_viirs_j01_chlor_a"), "NOAA-20", {"viirs"}, 16),
+            (lazy_fixture("seadas_l2_modis_chlor_a_netcdf"), "Terra", {"modis"}, 10),
         ])
     @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):
@@ -153,6 +233,7 @@ def test_load_chlor_a(self, input_files, exp_plat, exp_sensor, exp_rps, apply_qu
         scene = Scene(reader='seadas_l2', filenames=input_files, reader_kwargs=reader_kwargs)
         scene.load(['chlor_a'])
         data_arr = scene['chlor_a']
+        assert data_arr.dims == ("y", "x")
         assert data_arr.attrs['platform_name'] == exp_plat
         assert data_arr.attrs['sensor'] == exp_sensor
         assert data_arr.attrs['units'] == 'mg m^-3'