Merge pull request pytroll#2905 from pdebuyl/mcd12q1_draft

Mcd12q1 draft

satpy/etc/readers/mcd12q1.yaml

@@ -0,0 +1,49 @@
+reader:
+  name: mcd12q1
+  short_name: MCD12Q1
+  long_name: MODIS Level 3 (mcd12Q1) data in HDF-EOS format
+  description: MODIS HDF-EOS MCD12Q1 L3 Reader
+  status: Beta
+  supports_fsspec: false
+  reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader
+  sensors: [modis]
+
+file_types:
+  modis_mcd12q1_hdf_eos:
+    file_patterns: ['MCD12Q1.A{start_time:%Y%j}.{tile_id}.{collection:03d}.{production_time:%Y%j%H%M%S}.hdf']
+    file_reader: !!python/name:satpy.readers.mcd12q1.MCD12Q1HDFFileHandler
+
+datasets:
+  LC_Type1:
+    name: LC_Type1
+    resolution: 500
+    file_type: modis_mcd12q1_hdf_eos
+
+  LC_Type2:
+    name: LC_Type2
+    resolution: 500
+    file_type: modis_mcd12q1_hdf_eos
+  LC_Type3:
+    name: LC_Type3
+    resolution: 500
+    file_type: modis_mcd12q1_hdf_eos
+  LC_Type4:
+    name: LC_Type4
+    resolution: 500
+    file_type: modis_mcd12q1_hdf_eos
+  LC_Type5:
+    name: LC_Type5
+    resolution: 500
+    file_type: modis_mcd12q1_hdf_eos
+  LC_Prop1:
+    name: LC_Prop1
+    resolution: 500
+    file_type: modis_mcd12q1_hdf_eos
+  LC_Prop2:
+    name: LC_Prop2
+    resolution: 500
+    file_type: modis_mcd12q1_hdf_eos
+  LC_Prop3:
+    name: LC_Prop3
+    resolution: 500
+    file_type: modis_mcd12q1_hdf_eos

satpy/readers/hdfeos_base.py

@@ -89,6 +89,12 @@ def _find_and_run_interpolation(interpolation_functions, src_resolution, dst_res
     logger.debug("Interpolating from {} to {}".format(src_resolution, dst_resolution))
     return interpolation_function(*args)
 
+def _modis_date(date):
+    """Transform a date and time string into a datetime object."""
+    if len(date) == 19:
+        return dt.datetime.strptime(date, "%Y-%m-%d %H:%M:%S")
+    else:
+        return dt.datetime.strptime(date, "%Y-%m-%d %H:%M:%S.%f")
 
 class HDFEOSBaseFileReader(BaseFileHandler):
     """Base file handler for HDF EOS data for both L1b and L2 products."""
@@ -183,7 +189,7 @@ def start_time(self):
         try:
             date = (self.metadata["INVENTORYMETADATA"]["RANGEDATETIME"]["RANGEBEGINNINGDATE"]["VALUE"] + " " +
                     self.metadata["INVENTORYMETADATA"]["RANGEDATETIME"]["RANGEBEGINNINGTIME"]["VALUE"])
-            return dt.datetime.strptime(date, "%Y-%m-%d %H:%M:%S.%f")
+            return _modis_date(date)
         except KeyError:
             return self._start_time_from_filename()
 
@@ -196,7 +202,7 @@ def end_time(self):
         try:
             date = (self.metadata["INVENTORYMETADATA"]["RANGEDATETIME"]["RANGEENDINGDATE"]["VALUE"] + " " +
                     self.metadata["INVENTORYMETADATA"]["RANGEDATETIME"]["RANGEENDINGTIME"]["VALUE"])
-            return dt.datetime.strptime(date, "%Y-%m-%d %H:%M:%S.%f")
+            return _modis_date(date)
         except KeyError:
             return self.start_time
 

satpy/readers/mcd12q1.py

@@ -0,0 +1,144 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2024 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/>.
+"""MCD12Q1 hdf-eos format reader.
+
+Introduction
+------------
+
+The ``mcd12q1`` reader reads MCD12Q1 products in HDF-EOS format.
+
+The 500m product is provided  on a sinusoidal grid.
+
+Reference documents and links:
+    - MODIS land products grid: https://modis-land.gsfc.nasa.gov/MODLAND_grid.html
+    - User guide: https://lpdaac.usgs.gov/documents/101/MCD12_User_Guide_V6.pdf
+    - MCD12Q1 v061: MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500 m SIN Grid
+
+The reader has been tested with:
+    - MCD12Q1: Land cover data.
+
+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 HDFEOSBaseFileReader
+
+logger = logging.getLogger(__name__)
+
+
+class MCD12Q1HDFFileHandler(HDFEOSBaseFileReader):
+    """File handler for MCD12Q1 HDF-EOS 500m granules."""
+    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()
+
+        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
+                continue
+            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]]:
+        res = self._get_res()
+        for var_name in self.sd.datasets().keys():
+            if var_name in handled_var_names:
+                # skip variables that YAML had configured
+                continue
+            common = {"file_type": "mcd12q1_500m_hdf",
+                      "resolution": res,
+                      "name": var_name}
+            yield True, common
+
+
+    def _get_res(self):
+        """Compute the resolution from the file metadata."""
+        gridname = self.metadata["GridStructure"]["GRID_1"]["GridName"]
+        if "MCD12Q1" not in gridname:
+            raise ValueError("Only MCD12Q1 grids are supported")
+
+        resolution_string = self.metadata["ARCHIVEDMETADATA"]["NADIRDATARESOLUTION"]["VALUE"]
+        if resolution_string[-1] == "m":
+            return int(resolution_string.removesuffix("m"))
+        else:
+            raise ValueError("Cannot parse resolution of MCD12Q1 grid")
+
+
+    def get_dataset(self, dataset_id, dataset_info):
+        """Get DataArray for specified dataset."""
+        dataset_name = dataset_id["name"]
+        # xxx
+        dataset = self.load_dataset(dataset_name, dataset_info.pop("category", False))
+
+        self._add_satpy_metadata(dataset_id, dataset)
+
+        return dataset
+
+    def _get_area_extent(self):
+        """Get the grid properties."""
+        # Now compute the data extent
+        upperleft = self.metadata["GridStructure"]["GRID_1"]["UpperLeftPointMtrs"]
+        lowerright = self.metadata["GridStructure"]["GRID_1"]["LowerRightMtrs"]
+
+        return upperleft[0], lowerright[1], lowerright[0], upperleft[1]
+
+    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.
+
+        The proj_param string comes from https://lpdaac.usgs.gov/documents/101/MCD12_User_Guide_V6.pdf
+        """
+        proj_param = "proj=sinu +a=6371007.181 +b=6371007.181 +units=m"
+
+        # Get the size of the dataset
+        nrows = self.metadata["GridStructure"]["GRID_1"]["YDim"]
+        ncols = self.metadata["GridStructure"]["GRID_1"]["XDim"]
+
+        # Construct the area definition
+        area = geometry.AreaDefinition("sinusoidal_modis",
+                                       "Tiled sinusoidal L3 MODIS area",
+                                       "sinusoidal",
+                                       proj_param,
+                                       ncols,
+                                       nrows,
+                                       self._get_area_extent())
+
+        return area

satpy/tests/reader_tests/modis_tests/_modis_fixtures.py

@@ -215,6 +215,30 @@ def _get_l1b_geo_variable_info(filename: str,
         variables_info.update(_get_angles_variable_info(geo_resolution))
     return variables_info
 
+def _get_l3_land_cover_info() -> dict:
+
+    lc_data = np.zeros((2400, 2400), dtype=np.uint8)
+
+    variables_info = \
+    {
+        "LC_Type1": {"data": lc_data,
+                     "type": SDC.UINT8,
+                     "fill_value": 255,
+                     "attrs": {
+                         "dim_labels": ["YDim:MCD12Q1", "XDim:MCD12Q1"],
+                     },
+                     },
+        "LC_Type2": {"data": lc_data,
+                     "type": SDC.UINT8,
+                     "fill_value": 255,
+                     "attrs": {
+                         "dim_labels": ["YDim:MCD12Q1", "XDim:MCD12Q1"],
+                     },
+                     },
+    }
+
+    return variables_info
+
 
 def generate_nasa_l1b_filename(prefix):
     """Generate a filename that follows NASA MODIS L1b convention."""
@@ -331,19 +355,30 @@ def _create_struct_metadata_cmg(ftype: str) -> str:
         gridline = 'GridName="MOD09CMG"\n'
         upleft = "UpperLeftPointMtrs=(-180000000.000000,90000000.000000)\n"
         upright = "LowerRightMtrs=(180000000.000000,-90000000.000000)\n"
+        XDim=7200
+        YDim=3600
+    # Case of a MCD12Q1 file
+    elif ftype == "MCD12Q1":
+        gridline = 'GridName="MCD12Q1"\n'
+        upleft = "UpperLeftPointMtrs=(-8895604.157333,-1111950.519667)\n"
+        upright = "LowerRightMtrs=(-7783653.637667,-2223901.039333)\n"
+        XDim=2400
+        YDim=2400
     # Case of a MCD43 file
     else:
         gridline = 'GridName="MCD_CMG_BRDF_0.05Deg"\n'
         upleft = "UpperLeftPointMtrs=(-180.000000,90.000000)\n"
         upright = "LowerRightMtrs=(180.000000,-90.000000)\n"
+        XDim=7200
+        YDim=3600
 
     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"XDim={XDim}\n"
+                              f"YDim={YDim}\n"
                               f"{upleft}\n"
                               f"{upright}\n"
                               "END_GROUP=GRID_1\n"
@@ -598,6 +633,10 @@ def generate_nasa_l2_filename(prefix: str) -> str:
     now = dt.datetime.now()
     return f"{prefix}_L2.A{now:%Y%j.%H%M}.061.{now:%Y%j%H%M%S}.hdf"
 
+def generate_nasa_l3_tile_filename(prefix: str) -> str:
+    """Generate a file name that follows MODIS sinusoidal grid tile pattern."""
+    now = dt.datetime.now()
+    return f"{prefix}.A{now:%Y}001.h34v07.061.{now:%Y%j%H%M%S}.hdf"
 
 @pytest.fixture(scope="session")
 def modis_l2_nasa_mod35_file(tmpdir_factory) -> list[str]:
@@ -613,6 +652,31 @@ def modis_l2_nasa_mod35_file(tmpdir_factory) -> list[str]:
                             _create_header_metadata())
     return [full_path]
 
+@pytest.fixture(scope="session")
+def modis_l3_nasa_mcd12q1_file(tmpdir_factory) -> list[str]:
+    """Create a single MOD35 L2 HDF4 file with headers."""
+    filename = generate_nasa_l3_tile_filename("MCD12Q1")
+    full_path = str(tmpdir_factory.mktemp("modis_l2").join(filename))
+    variable_infos = _get_l3_land_cover_info()
+    archive_header = \
+"""GROUP                  = ARCHIVEDMETADATA
+  GROUPTYPE            = MASTERGROUP
+
+  OBJECT                 = NADIRDATARESOLUTION
+    NUM_VAL              = 1
+    VALUE                = "500m"
+  END_OBJECT             = NADIRDATARESOLUTION
+
+END_GROUP              = ARCHIVEDMETADATA
+
+END
+"""
+    create_hdfeos_test_file(full_path,
+                            variable_infos,
+                            _create_struct_metadata_cmg("MCD12Q1"),
+                            _create_core_metadata("MCD12Q1"),
+                            archive_header)
+    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."""

satpy/tests/reader_tests/modis_tests/conftest.py

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