Add a reader for MODIS Level 3 files in CMG format.

satpy/etc/readers/modis_l3.yaml

@@ -0,0 +1,16 @@
+reader:
+  name: modis_l3
+  short_name: MODIS l3
+  long_name: MODIS Level 3 (mcd43) data in HDF-EOS format
+  description: MODIS HDF-EOS L3 Reader
+  status: Beta
+  supports_fsspec: false
+  reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader
+  sensors: [modis]
+
+file_types:
+  modis_l3_cmg_hdf:
+    file_patterns:
+    - 'MCD43C{prod_type}.A{start_time:%Y%j}.{collection:03d}.{production_time:%Y%j%H%M%S}.hdf'
+    - 'M{platform_indicator:1s}D09CMG.A{start_time:%Y%j}.{collection:03d}.{production_time:%Y%j%H%M%S}.hdf'
+    file_reader: !!python/name:satpy.readers.modis_l3.ModisL3GriddedHDFFileHandler

satpy/readers/hdfeos_base.py

@@ -148,7 +148,10 @@ def _read_mda(cls, lines, element=None):
 
     @classmethod
     def _split_line(cls, line, lines):
-        key, val = line.split("=")
+        try:
+            key, val = line.split("=")
+        except ValueError:
+            key, val = line.split("=", maxsplit=1)
         key = key.strip()
         val = val.strip()
         try:

satpy/readers/modis_l3.py

@@ -0,0 +1,162 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2023 Satpy developers
+#
+# This file is part of satpy.
+#
+# satpy is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# satpy is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# satpy.  If not, see <http://www.gnu.org/licenses/>.
+"""Modis level 3 hdf-eos format reader.
+
+Introduction
+------------
+
+The ``modis_l3`` reader reads MODIS L3 products in HDF-EOS format.
+
+There are multiple level 3 products, including some on sinusoidal grids and some on the climate modeling grid (CMG).
+This reader supports the CMG products at present, and the sinusoidal products will be added if there is demand.
+
+The reader has been tested with:
+    - MCD43c*: BRDF/Albedo data, such as parameters, albedo and nbar
+    - MOD09CMG: Surface Reflectance on climate monitoring grid.
+
+To get a list of the available datasets for a given file refer to the "Load data" section in :doc:`../reading`.
+
+"""
+import logging
+from typing import Iterable
+
+from pyresample import geometry
+
+from satpy.readers.hdfeos_base import HDFEOSGeoReader
+
+logger = logging.getLogger(__name__)
+
+
+class ModisL3GriddedHDFFileHandler(HDFEOSGeoReader):
+    """File handler for MODIS HDF-EOS Level 3 CMG gridded files."""
+
+    def __init__(self, filename, filename_info, filetype_info, **kwargs):
+        """Init the file handler."""
+        super().__init__(filename, filename_info, filetype_info, **kwargs)
+
+        # Initialise number of rows and columns
+        self.nrows = self.metadata["GridStructure"]["GRID_1"]["YDim"]
+        self.ncols = self.metadata["GridStructure"]["GRID_1"]["XDim"]
+
+        # Get the grid name and other projection info
+        self.area_extent = self._sort_grid()
+
+
+    def _sort_grid(self):
+        """Get the grid properties."""
+        # First, get the grid resolution
+        self._get_res()
+
+        # Now compute the data extent
+        upperleft = self.metadata["GridStructure"]["GRID_1"]["UpperLeftPointMtrs"]
+        lowerright = self.metadata["GridStructure"]["GRID_1"]["LowerRightMtrs"]
+
+        # For some reason, a few of the CMG products multiply their
+        # decimal degree extents by one million. This fixes it.
+        if lowerright[0] > 1e6 or upperleft[0] > 1e6:
+            upperleft = tuple(val / 1e6 for val in upperleft)
+            lowerright = tuple(val / 1e6 for val in lowerright)
+
+        return upperleft[0], lowerright[1], lowerright[0], upperleft[1]
+
+
+    def _get_res(self):
+        """Compute the resolution from the file metadata."""
+        gridname = self.metadata["GridStructure"]["GRID_1"]["GridName"]
+        if "CMG" not in gridname:
+            raise ValueError("Only CMG grids are supported")
+
+        # Get the grid resolution from the grid name
+        pos = gridname.rfind("_") + 1
+        pos2 = gridname.rfind("Deg")
+
+        # Some products don't have resolution listed.
+        if pos < 0 or pos2 < 0:
+            self.resolution = 360. / self.ncols
+        else:
+            self.resolution = float(gridname[pos:pos2])
+
+
+    def available_datasets(self, configured_datasets=None):
+        """Automatically determine datasets provided by this file."""
+
+        # Initialise set of variable names to carry through code
+        handled_var_names = set()
+
+        yield from super().available_datasets(configured_datasets)
+
+        ds_dict = self.sd.datasets()
+
+        for is_avail, ds_info in (configured_datasets or []):
+            file_key = ds_info.get("file_key", ds_info["name"])
+            # we must add all variables here even if another file handler has
+            # claimed the variable. It could be another instance of this file
+            # type and we don't want to add that variable dynamically if the
+            # other file handler defined it by the YAML definition.
+            handled_var_names.add(file_key)
+            if is_avail is not None:
+                # some other file handler said it has this dataset
+                # we don't know any more information than the previous
+                # file handler so let's yield early
+                yield is_avail, ds_info
+                continue
+            if self.file_type_matches(ds_info["file_type"]) is None:
+                # this is not the file type for this dataset
+                yield None, ds_info
+            yield file_key in ds_dict.keys(), ds_info
+
+        yield from self._dynamic_variables_from_file(handled_var_names)
+
+    def _dynamic_variables_from_file(self, handled_var_names: set) -> Iterable[tuple[bool, dict]]:
+
+        for var_name in self.sd.datasets().keys():
+            if var_name in handled_var_names:
+                # skip variables that YAML had configured, but allow lon/lats
+                # to be reprocessed due to our dynamic coordinate naming
+                continue
+            common = {"file_type": "modis_l3_cmg_hdf",
+                      "resolution": self.resolution,
+                      "name": var_name}
+            yield True, common
+
+    def get_dataset(self, dataset_id, dataset_info):
+        """Get DataArray for specified dataset."""
+        dataset_name = dataset_id["name"]
+        dataset = self.load_dataset(dataset_name, dataset_info.pop("category", False))
+        self._add_satpy_metadata(dataset_id, dataset)
+
+        return dataset
+
+
+    def get_area_def(self, dsid):
+        """Get the area definition.
+
+        This is fixed, but not defined in the file. So we must
+        generate it ourselves with some assumptions.
+        """
+        proj_param = "EPSG:4326"
+
+        area = geometry.AreaDefinition("gridded_modis",
+                                       "A gridded L3 MODIS area",
+                                       "longlat",
+                                       proj_param,
+                                       self.ncols,
+                                       self.nrows,
+                                       self.area_extent)
+
+        return area

satpy/tests/reader_tests/modis_tests/_modis_fixtures.py

@@ -51,6 +51,9 @@
 
 
 def _shape_for_resolution(resolution: int) -> tuple[int, int]:
+    # Case of a CMG 0.05 degree file, for L3 tests
+    if resolution == -999:
+        return 3600, 7200
     assert resolution in RES_TO_REPEAT_FACTOR
     factor = RES_TO_REPEAT_FACTOR[resolution]
     if factor == 1:
@@ -226,6 +229,16 @@
     return f"t1.{now:%y%j.%H%M}.{suffix}.hdf"
 
 
+def _add_geo_metadata(h, geo_res):
+    """Add the geoinfo metadata to the fake file."""
+    if geo_res == -999 or geo_res == -9999:
+        setattr(h, 'StructMetadata.0', _create_struct_metadata_cmg(geo_res))  # noqa
+    else:
+        setattr(h, 'StructMetadata.0', _create_struct_metadata(geo_res))  # noqa
+
+    return h
+
+
 def create_hdfeos_test_file(filename: str,
                             variable_infos: dict,
                             geo_resolution: Optional[int] = None,
@@ -251,8 +264,8 @@
     if include_metadata:
         if geo_resolution is None or file_shortname is None:
             raise ValueError("'geo_resolution' and 'file_shortname' are required when including metadata.")
+        h = _add_geo_metadata(h, geo_resolution)
         setattr(h, 'CoreMetadata.0', _create_core_metadata(file_shortname))  # noqa
-        setattr(h, 'StructMetadata.0', _create_struct_metadata(geo_resolution))  # noqa
         setattr(h, 'ArchiveMetadata.0', _create_header_metadata())  # noqa
 
     for var_name, var_info in variable_infos.items():
@@ -326,8 +339,35 @@
     return struct_metadata_header
 
 
+def _create_struct_metadata_cmg(res) -> str:
+    # Case of a MOD09 file
+    gridline = 'GridName="MOD09CMG"\n'
+    upleft = "UpperLeftPointMtrs=(-180000000.000000,90000000.000000)\n"
+    upright = "LowerRightMtrs=(180000000.000000,-90000000.000000)\n"
+    if res == -9999:
+        # Case of a MCD43 file
+        gridline = 'GridName="MCD_CMG_BRDF_0.05Deg"\n'
+        upleft = "UpperLeftPointMtrs=(-180.000000,90.000000)\n"
+        upright = "LowerRightMtrs=(180.000000,-90.000000)\n"
+
+    struct_metadata_header = ("GROUP=SwathStructure\n"
+                              "END_GROUP=SwathStructure\n"
+                              "GROUP=GridStructure\n"
+                              "GROUP=GRID_1\n"
+                              f"{gridline}\n"
+                              "XDim=7200\n"
+                              "YDim=3600\n"
+                              f"{upleft}\n"
+                              f"{upright}\n"
+                              "END_GROUP=GRID_1\n"
+                              "END_GROUP=GridStructure\nEND")
+    return struct_metadata_header
+
+
 def _create_header_metadata() -> str:
-    archive_metadata_header = "GROUP = ARCHIVEDMETADATA\nEND_GROUP = ARCHIVEDMETADATA\nEND"
+    archive_metadata_header = ("GROUP = ARCHIVEDMETADATA\n"
+                               'TEST_URL = "http://modis.gsfc.nasa.gov/?some_val=100"\n'
+                               "END_GROUP = ARCHIVEDMETADATA\nEND")
     return archive_metadata_header
 
 
@@ -471,6 +511,28 @@
     }
 
 
+def _get_l3_refl_variable_info(var_name: str) -> dict:
+    shape = (3600, 7200)
+    data = np.zeros((shape[0], shape[1]), dtype=np.int16)
+    row_dim_name = "XDim"
+    col_dim_name = "YDim"
+    return {
+        var_name: {
+            "data": data,
+            "type": SDC.INT16,
+            "fill_value": -28672,
+            "attrs": {
+                # dim_labels are just unique dimension names, may not match exactly with real world files
+                "dim_labels": [row_dim_name,
+                               col_dim_name],
+                "valid_range": (-100, 16000),
+                "scale_factor": 1e-4,
+                "add_offset": 0.,
+            },
+        },
+    }
+
+
 def _get_mask_byte1_variable_info() -> dict:
     shape = _shape_for_resolution(1000)
     data = np.zeros((shape[0], shape[1]), dtype=np.uint16)
@@ -537,6 +599,41 @@
     return [full_path]
 
 
+def generate_nasa_l3_filename(prefix: str) -> str:
+    """Generate a file name that follows MODIS 09 L3 convention in a temporary directory."""
+    now = datetime.now()
+    return f"{prefix}.A{now:%Y%j}.061.{now:%Y%j%H%M%S}.hdf"
+
+
+def modis_l3_file(tmpdir_factory, f_prefix, var_name, geo_res, f_short):
+    """Create a MODIS L3 file of the desired type."""
+    filename = generate_nasa_l3_filename(f_prefix)
+    full_path = str(tmpdir_factory.mktemp("modis_l3").join(filename))
+    variable_infos = _get_l3_refl_variable_info(var_name)
+    create_hdfeos_test_file(full_path, variable_infos, geo_resolution=geo_res, file_shortname=f_short)
+    return [full_path]
+
+
+@pytest.fixture(scope="session")
+def modis_l3_nasa_mod09_file(tmpdir_factory) -> list[str]:
+    """Create a single MOD09 L3 HDF4 file with headers."""
+    return modis_l3_file(tmpdir_factory,
+                         "MOD09CMG",
+                         "Coarse_Resolution_Surface_Reflectance_Band_2",
+                         -999,
+                         "MOD09")
+
+
+@pytest.fixture(scope="session")
+def modis_l3_nasa_mod43_file(tmpdir_factory) -> list[str]:
+    """Create a single MVCD43 L3 HDF4 file with headers."""
+    return modis_l3_file(tmpdir_factory,
+                         "MCD43C1",
+                         "BRDF_Albedo_Parameter1_Band2",
+                         -9999,
+                         "MCD43C1")
+
+
 @pytest.fixture(scope="session")
 def modis_l2_nasa_mod35_mod03_files(modis_l2_nasa_mod35_file, modis_l1b_nasa_mod03_file) -> list[str]:
     """Create a MOD35 L2 HDF4 file and MOD03 L1b geolocation file."""

satpy/tests/reader_tests/modis_tests/conftest.py

@@ -32,4 +32,6 @@
     modis_l2_nasa_mod06_file,
     modis_l2_nasa_mod35_file,
     modis_l2_nasa_mod35_mod03_files,
+    modis_l3_nasa_mod09_file,
+    modis_l3_nasa_mod43_file,
 )