From 9f9e34ddad87fb69ebfdf9b4b89cd1f5872938fb Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 10:20:18 +0000 Subject: [PATCH 01/30] Add a reader for MODIS Level 3 files in CMG format. --- satpy/readers/hdfeos_base.py | 5 +- satpy/readers/modis_l3.py | 111 +++++++++++++++++++++++++++++++++++ 2 files changed, 115 insertions(+), 1 deletion(-) create mode 100644 satpy/readers/modis_l3.py diff --git a/satpy/readers/hdfeos_base.py b/satpy/readers/hdfeos_base.py index f60040a46f..9964eeb2e1 100644 --- a/satpy/readers/hdfeos_base.py +++ b/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: diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py new file mode 100644 index 0000000000..3c60efba87 --- /dev/null +++ b/satpy/readers/modis_l3.py @@ -0,0 +1,111 @@ +#!/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 . +"""Modis level 3 hdf-eos format reader. + +Introduction +------------ + +The ``modis_l3`` reader reads Modis L3 products in hdf-eos format. +Since there are a multitude of different level 3 datasets not all of theses are implemented (yet). + + +Currently the reader supports: + - mcd43c1: BRDF/Albedo Model Parameters dataset + - mcd43c3: BRDF/Albedo Albedo dataset + +To get a list of the available datasets for a given file refer to the "Load data" section in :doc:`../reading`. + +""" +import logging + +from pyresample import geometry + +from satpy.readers.hdfeos_base import HDFEOSGeoReader +from satpy.utils import get_legacy_chunk_size + +logger = logging.getLogger(__name__) +CHUNK_SIZE = get_legacy_chunk_size() + + +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 + gridname = self.metadata["GridStructure"]["GRID_1"]["GridName"] + if "CMG" not in gridname: + raise ValueError("Only CMG grids are supported") + + # Get the grid resolution + pos = gridname.rfind("_") + 1 + pos2 = gridname.rfind("Deg") + self.resolution = float(gridname[pos:pos2]) + + upperleft = self.metadata["GridStructure"]["GRID_1"]["UpperLeftPointMtrs"] + lowerright = self.metadata["GridStructure"]["GRID_1"]["LowerRightMtrs"] + + self.area_extent = (upperleft[0], lowerright[1], lowerright[0], upperleft[1]) + + + def available_datasets(self, configured_datasets=None): + """Automatically determine datasets provided by this file.""" + logger.debug("Available_datasets begin...") + + ds_dict = self.sd.datasets() + + yield from super().available_datasets(configured_datasets) + common = {"file_type": "mcd43_cmg_hdf", "resolution": self.resolution} + for key in ds_dict.keys(): + if "/" in key: # not a dataset + continue + yield True, {"name": key} | 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) + self.area = area + + return self.area From 7c7d5fc5670e76574f815396c239d552688b9dfd Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 10:40:04 +0000 Subject: [PATCH 02/30] Update reader for MODIS Level 3 files in CMG format. --- satpy/etc/readers/modis_l3.yaml | 16 ++++++++++++++++ satpy/readers/modis_l3.py | 13 ++++++++++++- 2 files changed, 28 insertions(+), 1 deletion(-) create mode 100644 satpy/etc/readers/modis_l3.yaml diff --git a/satpy/etc/readers/modis_l3.yaml b/satpy/etc/readers/modis_l3.yaml new file mode 100644 index 0000000000..5ad2f32e04 --- /dev/null +++ b/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: + mcd43_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 diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 3c60efba87..2b9387ed58 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -61,11 +61,22 @@ def __init__(self, filename, filename_info, filetype_info, **kwargs): # Get the grid resolution pos = gridname.rfind("_") + 1 pos2 = gridname.rfind("Deg") - self.resolution = float(gridname[pos:pos2]) + + # 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]) 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: + upperleft = tuple(val / 1e6 for val in upperleft) + lowerright = tuple(val / 1e6 for val in lowerright) + self.area_extent = (upperleft[0], lowerright[1], lowerright[0], upperleft[1]) From 0de217a4c7f7db8e343d6e1173ca0103cc37c57a Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 10:44:15 +0000 Subject: [PATCH 03/30] Update MODIS L3 docstring. --- satpy/readers/modis_l3.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 2b9387ed58..86560624a8 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -21,12 +21,13 @@ ------------ The ``modis_l3`` reader reads Modis L3 products in hdf-eos format. -Since there are a multitude of different level 3 datasets not all of theses are implemented (yet). +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. -Currently the reader supports: - - mcd43c1: BRDF/Albedo Model Parameters dataset - - mcd43c3: BRDF/Albedo Albedo dataset +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`. From b4cf2c43447a6c87ceb9b1a767b0ffc5b036e3d0 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 10:47:12 +0000 Subject: [PATCH 04/30] Restructure the L3 MODIS reader. --- satpy/readers/modis_l3.py | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 86560624a8..bff33b190c 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -46,15 +46,8 @@ 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 + def _sort_grid(self): + """Get the grid properties.""" gridname = self.metadata["GridStructure"]["GRID_1"]["GridName"] if "CMG" not in gridname: raise ValueError("Only CMG grids are supported") @@ -65,9 +58,9 @@ def __init__(self, filename, filename_info, filetype_info, **kwargs): # Some products don't have resolution listed. if pos < 0 or pos2 < 0: - self.resolution = 360. / self.ncols + resolution = 360. / self.ncols else: - self.resolution = float(gridname[pos:pos2]) + resolution = float(gridname[pos:pos2]) upperleft = self.metadata["GridStructure"]["GRID_1"]["UpperLeftPointMtrs"] lowerright = self.metadata["GridStructure"]["GRID_1"]["LowerRightMtrs"] @@ -78,7 +71,21 @@ def __init__(self, filename, filename_info, filetype_info, **kwargs): upperleft = tuple(val / 1e6 for val in upperleft) lowerright = tuple(val / 1e6 for val in lowerright) - self.area_extent = (upperleft[0], lowerright[1], lowerright[0], upperleft[1]) + area_extent = (upperleft[0], lowerright[1], lowerright[0], upperleft[1]) + + return resolution, area_extent + + + 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.resolution, self.area_extent = self._sort_grid() def available_datasets(self, configured_datasets=None): From 1457ac67e7457ab1c05d2472a73c0b632864278a Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 11:44:39 +0000 Subject: [PATCH 05/30] Add tests for MODIS L3 reader. --- .../modis_tests/_modis_fixtures.py | 80 ++++++++++++++++++- .../reader_tests/modis_tests/conftest.py | 2 + .../reader_tests/modis_tests/test_modis_l3.py | 78 ++++++++++++++++++ 3 files changed, 159 insertions(+), 1 deletion(-) create mode 100644 satpy/tests/reader_tests/modis_tests/test_modis_l3.py diff --git a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py index 66221c613e..aff84de7be 100644 --- a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py +++ b/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: @@ -252,7 +255,10 @@ def create_hdfeos_test_file(filename: str, if geo_resolution is None or file_shortname is None: raise ValueError("'geo_resolution' and 'file_shortname' are required when including metadata.") setattr(h, 'CoreMetadata.0', _create_core_metadata(file_shortname)) # noqa - setattr(h, 'StructMetadata.0', _create_struct_metadata(geo_resolution)) # noqa + if geo_resolution == -999 or geo_resolution == -9999: + setattr(h, 'StructMetadata.0', _create_struct_metadata_cmg(geo_resolution)) # noqa + else: + 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,6 +332,31 @@ def _create_struct_metadata(geo_resolution: int) -> str: 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" return archive_metadata_header @@ -471,6 +502,28 @@ def _get_cloud_mask_variable_info(var_name: str, resolution: int) -> dict: } +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 +590,31 @@ def modis_l2_nasa_mod35_file(tmpdir_factory) -> list[str]: 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" + + +@pytest.fixture(scope="session") +def modis_l3_nasa_mod09_file(tmpdir_factory) -> list[str]: + """Create a single MOD09 L3 HDF4 file with headers.""" + filename = generate_nasa_l3_filename("MOD09CMG") + full_path = str(tmpdir_factory.mktemp("modis_l3").join(filename)) + variable_infos = _get_l3_refl_variable_info("Coarse_Resolution_Surface_Reflectance_Band_2") + create_hdfeos_test_file(full_path, variable_infos, geo_resolution=-999, file_shortname="MOD09") + return [full_path] + +@pytest.fixture(scope="session") +def modis_l3_nasa_mod43_file(tmpdir_factory) -> list[str]: + """Create a single MOD09 L3 HDF4 file with headers.""" + filename = generate_nasa_l3_filename("MCD43C1") + full_path = str(tmpdir_factory.mktemp("modis_l3").join(filename)) + variable_infos = _get_l3_refl_variable_info("BRDF_Albedo_Parameter1_Band2") + create_hdfeos_test_file(full_path, variable_infos, geo_resolution=-9999, file_shortname="MCD43C1") + return [full_path] + + @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.""" diff --git a/satpy/tests/reader_tests/modis_tests/conftest.py b/satpy/tests/reader_tests/modis_tests/conftest.py index e6a8432653..309b16321f 100644 --- a/satpy/tests/reader_tests/modis_tests/conftest.py +++ b/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, ) diff --git a/satpy/tests/reader_tests/modis_tests/test_modis_l3.py b/satpy/tests/reader_tests/modis_tests/test_modis_l3.py new file mode 100644 index 0000000000..1203ecf205 --- /dev/null +++ b/satpy/tests/reader_tests/modis_tests/test_modis_l3.py @@ -0,0 +1,78 @@ +#!/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 . +"""Unit tests for MODIS L3 HDF reader.""" + +from __future__ import annotations + +import dask.array as da +import pytest +from pyresample import geometry +from pytest_lazyfixture import lazy_fixture + +from satpy import Scene, available_readers + +from ._modis_fixtures import _shape_for_resolution + + +def _expected_area(): + proj_param = "EPSG:4326" + + return geometry.AreaDefinition("gridded_modis", + "A gridded L3 MODIS area", + "longlat", + proj_param, + 7200, + 3600, + (-180, -90, 180, 90)) + + +class TestModisL3: + """Test MODIS L3 reader.""" + + def test_available_reader(self): + """Test that MODIS L3 reader is available.""" + assert "modis_l3" in available_readers() + + @pytest.mark.parametrize( + ("loadable", "filename"), + [ + ("Coarse_Resolution_Surface_Reflectance_Band_2", lazy_fixture("modis_l3_nasa_mod09_file")), + ("BRDF_Albedo_Parameter1_Band2",lazy_fixture("modis_l3_nasa_mod43_file")), + ] + ) + def test_scene_available_datasets(self, loadable, filename): + """Test that datasets are available.""" + scene = Scene(reader="modis_l3", filenames=filename) + available_datasets = scene.all_dataset_names() + assert len(available_datasets) > 0 + assert loadable in available_datasets + + def test_load_l3_dataset(self, modis_l3_nasa_mod09_file): + """Load and check an L2 variable.""" + scene = Scene(reader="modis_l3", filenames=modis_l3_nasa_mod09_file) + + ds_name = "Coarse_Resolution_Surface_Reflectance_Band_2" + scene.load([ds_name]) + + data_arr = scene[ds_name] + assert isinstance(data_arr.data, da.Array) + data_arr = data_arr.compute() + + assert data_arr.shape == _shape_for_resolution(-999) + assert data_arr.attrs.get("resolution") == 0.05 + assert data_arr.attrs.get("area") == _expected_area() From 1f501168138ab74c143975fc160ab5fc1129119d Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 11:51:39 +0000 Subject: [PATCH 06/30] Simplify the MODIS L3 code + tests somewhat. --- satpy/readers/modis_l3.py | 25 ++++++++++++------- .../modis_tests/_modis_fixtures.py | 4 +-- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index bff33b190c..dc4600790d 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -46,22 +46,31 @@ class ModisL3GriddedHDFFileHandler(HDFEOSGeoReader): """File handler for MODIS HDF-EOS Level 3 CMG gridded files.""" - def _sort_grid(self): - """Get the grid properties.""" + 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 + # 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: - resolution = 360. / self.ncols + self.resolution = 360. / self.ncols else: - resolution = float(gridname[pos:pos2]) + self.resolution = float(gridname[pos:pos2]) + + + 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"] @@ -71,9 +80,7 @@ def _sort_grid(self): upperleft = tuple(val / 1e6 for val in upperleft) lowerright = tuple(val / 1e6 for val in lowerright) - area_extent = (upperleft[0], lowerright[1], lowerright[0], upperleft[1]) - - return resolution, area_extent + self.area_extent = (upperleft[0], lowerright[1], lowerright[0], upperleft[1]) def __init__(self, filename, filename_info, filetype_info, **kwargs): @@ -85,7 +92,7 @@ def __init__(self, filename, filename_info, filetype_info, **kwargs): self.ncols = self.metadata["GridStructure"]["GRID_1"]["XDim"] # Get the grid name and other projection info - self.resolution, self.area_extent = self._sort_grid() + self._sort_grid() def available_datasets(self, configured_datasets=None): diff --git a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py index aff84de7be..f076e100ba 100644 --- a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py @@ -254,11 +254,11 @@ def create_hdfeos_test_file(filename: str, 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.") - setattr(h, 'CoreMetadata.0', _create_core_metadata(file_shortname)) # noqa - if geo_resolution == -999 or geo_resolution == -9999: + elif geo_resolution == -999 or geo_resolution == -9999: setattr(h, 'StructMetadata.0', _create_struct_metadata_cmg(geo_resolution)) # noqa else: setattr(h, 'StructMetadata.0', _create_struct_metadata(geo_resolution)) # noqa + setattr(h, 'CoreMetadata.0', _create_core_metadata(file_shortname)) # noqa setattr(h, 'ArchiveMetadata.0', _create_header_metadata()) # noqa for var_name, var_info in variable_infos.items(): From 7afa5773f4e03f478e4725a8bb091e5602fbdd8c Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 11:57:40 +0000 Subject: [PATCH 07/30] Further simplify the MODIS L3 tests. --- .../modis_tests/_modis_fixtures.py | 32 ++++++++++++------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py index f076e100ba..40e448e067 100644 --- a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py @@ -596,23 +596,33 @@ def generate_nasa_l3_filename(prefix: str) -> str: 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.""" - filename = generate_nasa_l3_filename("MOD09CMG") - full_path = str(tmpdir_factory.mktemp("modis_l3").join(filename)) - variable_infos = _get_l3_refl_variable_info("Coarse_Resolution_Surface_Reflectance_Band_2") - create_hdfeos_test_file(full_path, variable_infos, geo_resolution=-999, file_shortname="MOD09") - return [full_path] + 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 MOD09 L3 HDF4 file with headers.""" - filename = generate_nasa_l3_filename("MCD43C1") - full_path = str(tmpdir_factory.mktemp("modis_l3").join(filename)) - variable_infos = _get_l3_refl_variable_info("BRDF_Albedo_Parameter1_Band2") - create_hdfeos_test_file(full_path, variable_infos, geo_resolution=-9999, file_shortname="MCD43C1") - return [full_path] + """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") From 35434d3ae172d795c0716e1125077d285fcbdf3a Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 15:40:01 +0000 Subject: [PATCH 08/30] Update satpy/readers/modis_l3.py Co-authored-by: David Hoese --- satpy/readers/modis_l3.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index dc4600790d..b5313d89cc 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -61,9 +61,6 @@ def _get_res(self): self.resolution = 360. / self.ncols else: self.resolution = float(gridname[pos:pos2]) - - - def _sort_grid(self): """Get the grid properties.""" From fb59a77e279f2cf8d9447ab8429ca1075ebc6bcf Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 15:41:39 +0000 Subject: [PATCH 09/30] Update MODIS L3 reader for review comments. --- satpy/readers/modis_l3.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index b5313d89cc..f8f94372cc 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -37,10 +37,8 @@ from pyresample import geometry from satpy.readers.hdfeos_base import HDFEOSGeoReader -from satpy.utils import get_legacy_chunk_size logger = logging.getLogger(__name__) -CHUNK_SIZE = get_legacy_chunk_size() class ModisL3GriddedHDFFileHandler(HDFEOSGeoReader): @@ -61,6 +59,7 @@ def _get_res(self): self.resolution = 360. / self.ncols else: self.resolution = float(gridname[pos:pos2]) + def _sort_grid(self): """Get the grid properties.""" @@ -77,7 +76,7 @@ def _sort_grid(self): upperleft = tuple(val / 1e6 for val in upperleft) lowerright = tuple(val / 1e6 for val in lowerright) - self.area_extent = (upperleft[0], lowerright[1], lowerright[0], upperleft[1]) + return (upperleft[0], lowerright[1], lowerright[0], upperleft[1]) def __init__(self, filename, filename_info, filetype_info, **kwargs): @@ -89,7 +88,7 @@ def __init__(self, filename, filename_info, filetype_info, **kwargs): self.ncols = self.metadata["GridStructure"]["GRID_1"]["XDim"] # Get the grid name and other projection info - self._sort_grid() + self.area_extent = self._sort_grid() def available_datasets(self, configured_datasets=None): From 261fca6e38cba8a3a0fe145f8a5078e90bb25723 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 15:46:13 +0000 Subject: [PATCH 10/30] Update MODIS test fixtures to simplify. --- .../modis_tests/_modis_fixtures.py | 35 +++++++++++-------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py index 40e448e067..e4272373b3 100644 --- a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py @@ -229,6 +229,16 @@ def generate_imapp_filename(suffix): 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, @@ -254,10 +264,7 @@ def create_hdfeos_test_file(filename: str, 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.") - elif geo_resolution == -999 or geo_resolution == -9999: - setattr(h, 'StructMetadata.0', _create_struct_metadata_cmg(geo_resolution)) # noqa - else: - setattr(h, 'StructMetadata.0', _create_struct_metadata(geo_resolution)) # noqa + h = _add_geo_metadata(h, geo_resolution) setattr(h, 'CoreMetadata.0', _create_core_metadata(file_shortname)) # noqa setattr(h, 'ArchiveMetadata.0', _create_header_metadata()) # noqa @@ -344,16 +351,16 @@ def _create_struct_metadata_cmg(res) -> str: 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") + "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 From bdcddab6ca033d26d2363513ddb24fc07ddb2679 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 15:46:40 +0000 Subject: [PATCH 11/30] Remove rogue blank line. --- satpy/readers/modis_l3.py | 1 - 1 file changed, 1 deletion(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index f8f94372cc..dfddc0732b 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -62,7 +62,6 @@ def _get_res(self): def _sort_grid(self): """Get the grid properties.""" - # First, get the grid resolution self._get_res() From 9fd639475113feaf7fdb884bcd8087f409b7c0e1 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 8 Nov 2023 16:47:23 +0000 Subject: [PATCH 12/30] Update MODIS L3 tests to check data types. --- .../tests/reader_tests/modis_tests/test_modis_l3.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/satpy/tests/reader_tests/modis_tests/test_modis_l3.py b/satpy/tests/reader_tests/modis_tests/test_modis_l3.py index 1203ecf205..23c1af6fc1 100644 --- a/satpy/tests/reader_tests/modis_tests/test_modis_l3.py +++ b/satpy/tests/reader_tests/modis_tests/test_modis_l3.py @@ -20,6 +20,7 @@ from __future__ import annotations import dask.array as da +import numpy as np import pytest from pyresample import geometry from pytest_lazyfixture import lazy_fixture @@ -71,8 +72,12 @@ def test_load_l3_dataset(self, modis_l3_nasa_mod09_file): data_arr = scene[ds_name] assert isinstance(data_arr.data, da.Array) - data_arr = data_arr.compute() + data_arr_comp = data_arr.compute() - assert data_arr.shape == _shape_for_resolution(-999) - assert data_arr.attrs.get("resolution") == 0.05 - assert data_arr.attrs.get("area") == _expected_area() + # Check types + assert data_arr_comp.dtype == data_arr.dtype + assert data_arr_comp.dtype == np.float32 + + assert data_arr_comp.shape == _shape_for_resolution(-999) + assert data_arr_comp.attrs.get("resolution") == 0.05 + assert data_arr_comp.attrs.get("area") == _expected_area() From d8c594b53255d4f82b756de807c4d5f434280997 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 14:45:15 +0000 Subject: [PATCH 13/30] Add line to HDF-EOS tests to ensure that line splitting in attrs is handled correctly. --- satpy/tests/reader_tests/modis_tests/_modis_fixtures.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py index e4272373b3..0b79e00854 100644 --- a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py @@ -365,7 +365,9 @@ def _create_struct_metadata_cmg(res) -> str: 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 From 19d1b9969c9464d65735f5ce867deadac7e1bfc6 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 14:47:16 +0000 Subject: [PATCH 14/30] Fix docstring in MODIS L3 reader. --- satpy/readers/modis_l3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index dfddc0732b..68c7b435ed 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -20,7 +20,7 @@ Introduction ------------ -The ``modis_l3`` reader reads Modis L3 products in hdf-eos format. +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. From 59e22405c03b408e67f067272b6369aa9a66fa24 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 14:49:46 +0000 Subject: [PATCH 15/30] Re-order functions in the MODIS L3 reader. --- satpy/readers/modis_l3.py | 45 ++++++++++++++++++++------------------- 1 file changed, 23 insertions(+), 22 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 68c7b435ed..2055d041d9 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -44,21 +44,17 @@ class ModisL3GriddedHDFFileHandler(HDFEOSGeoReader): """File handler for MODIS HDF-EOS Level 3 CMG gridded files.""" - 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") + def __init__(self, filename, filename_info, filetype_info, **kwargs): + """Init the file handler.""" + super().__init__(filename, filename_info, filetype_info, **kwargs) - # Get the grid resolution from the grid name - pos = gridname.rfind("_") + 1 - pos2 = gridname.rfind("Deg") + # 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() - # 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 _sort_grid(self): """Get the grid properties.""" @@ -75,19 +71,24 @@ def _sort_grid(self): 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]) + return upperleft[0], lowerright[1], lowerright[0], upperleft[1] - def __init__(self, filename, filename_info, filetype_info, **kwargs): - """Init the file handler.""" - super().__init__(filename, filename_info, filetype_info, **kwargs) + 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") - # 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 resolution from the grid name + pos = gridname.rfind("_") + 1 + pos2 = gridname.rfind("Deg") - # Get the grid name and other projection info - self.area_extent = self._sort_grid() + # 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): From 876d683d758986acc9f277b6266de774eaeca752 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 14:50:55 +0000 Subject: [PATCH 16/30] Re-order code in the `available_datasets` function for MODIS L3. --- satpy/readers/modis_l3.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 2055d041d9..ef0a55975e 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -93,12 +93,11 @@ def _get_res(self): def available_datasets(self, configured_datasets=None): """Automatically determine datasets provided by this file.""" - logger.debug("Available_datasets begin...") - - ds_dict = self.sd.datasets() yield from super().available_datasets(configured_datasets) common = {"file_type": "mcd43_cmg_hdf", "resolution": self.resolution} + ds_dict = self.sd.datasets() + for key in ds_dict.keys(): if "/" in key: # not a dataset continue From 60a4a297427d37252a0a7dab766e7437549ace53 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 14:53:01 +0000 Subject: [PATCH 17/30] Remove `self.area` from the MODIS L3 code. --- satpy/readers/modis_l3.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index ef0a55975e..316fa5b9ae 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -127,6 +127,5 @@ def get_area_def(self, dsid): self.ncols, self.nrows, self.area_extent) - self.area = area - return self.area + return area From 01f09aa10ba029d7c05b6703b07f0963fc738a40 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 14:53:41 +0000 Subject: [PATCH 18/30] Add missing space in MODIS L3 tests. --- satpy/tests/reader_tests/modis_tests/test_modis_l3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/tests/reader_tests/modis_tests/test_modis_l3.py b/satpy/tests/reader_tests/modis_tests/test_modis_l3.py index 23c1af6fc1..c7a0d3e9cf 100644 --- a/satpy/tests/reader_tests/modis_tests/test_modis_l3.py +++ b/satpy/tests/reader_tests/modis_tests/test_modis_l3.py @@ -53,7 +53,7 @@ def test_available_reader(self): ("loadable", "filename"), [ ("Coarse_Resolution_Surface_Reflectance_Band_2", lazy_fixture("modis_l3_nasa_mod09_file")), - ("BRDF_Albedo_Parameter1_Band2",lazy_fixture("modis_l3_nasa_mod43_file")), + ("BRDF_Albedo_Parameter1_Band2", lazy_fixture("modis_l3_nasa_mod43_file")), ] ) def test_scene_available_datasets(self, loadable, filename): From ef16a1a46105c23a2188586511d159dcd5472f24 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 15:07:54 +0000 Subject: [PATCH 19/30] Rename MODIS L3 file type to be more generic and refactor the dynamic dataset handler. --- satpy/etc/readers/modis_l3.yaml | 2 +- satpy/readers/modis_l3.py | 39 +++++++++++++++++++++++++++++---- 2 files changed, 36 insertions(+), 5 deletions(-) diff --git a/satpy/etc/readers/modis_l3.yaml b/satpy/etc/readers/modis_l3.yaml index 5ad2f32e04..608d15601a 100644 --- a/satpy/etc/readers/modis_l3.yaml +++ b/satpy/etc/readers/modis_l3.yaml @@ -9,7 +9,7 @@ reader: sensors: [modis] file_types: - mcd43_cmg_hdf: + 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' diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 316fa5b9ae..2c2f331c29 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -33,6 +33,7 @@ """ import logging +from typing import Iterable from pyresample import geometry @@ -94,14 +95,44 @@ def _get_res(self): 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) - common = {"file_type": "mcd43_cmg_hdf", "resolution": self.resolution} + ds_dict = self.sd.datasets() - for key in ds_dict.keys(): - if "/" in key: # not a dataset + 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 - yield True, {"name": key} | common + 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.""" From 9794f3df5a49546323d3ac7644cff308d8ac8204 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 15:18:38 +0000 Subject: [PATCH 20/30] Add additional check on MODIS L3 bounds to ensure values are scaled correctly. --- satpy/readers/modis_l3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 2c2f331c29..4ff0bd1dc6 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -68,7 +68,7 @@ def _sort_grid(self): # For some reason, a few of the CMG products multiply their # decimal degree extents by one million. This fixes it. - if lowerright[0] > 1e6: + if lowerright[0] > 1e6 or upperleft[0] > 1e6: upperleft = tuple(val / 1e6 for val in upperleft) lowerright = tuple(val / 1e6 for val in lowerright) From 23f2e0f97bc477adfdfb5eae63bc94f1f1ac1b49 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 15:45:59 +0000 Subject: [PATCH 21/30] Remove rogue `yield from` in MODIS L3 reader. --- satpy/readers/modis_l3.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 4ff0bd1dc6..1a205545ac 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -98,8 +98,6 @@ def available_datasets(self, configured_datasets=None): # 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 []): From a0e8386b6e332e770bba2edb43ad5a192dc90aca Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 15:56:22 +0000 Subject: [PATCH 22/30] Properly handle incorrect datasets for the VIIRS EDR and MODIS --- satpy/readers/modis_l3.py | 1 + satpy/readers/viirs_edr.py | 1 + 2 files changed, 2 insertions(+) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 1a205545ac..2370609e7a 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -116,6 +116,7 @@ def available_datasets(self, configured_datasets=None): 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) diff --git a/satpy/readers/viirs_edr.py b/satpy/readers/viirs_edr.py index 646d7e0d17..b0eaf7b7ba 100644 --- a/satpy/readers/viirs_edr.py +++ b/satpy/readers/viirs_edr.py @@ -215,6 +215,7 @@ def available_datasets(self, configured_datasets=None): 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 self.nc, ds_info yield from self._dynamic_variables_from_file(handled_var_names) From b748479cf9d197b48d2bf830c6edc9cfbcfd030c Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 16:16:18 +0000 Subject: [PATCH 23/30] Simplify the MODIS L3 code and add tests for the dynamic dataset availability. --- satpy/readers/modis_l3.py | 62 +++++++++---------- .../reader_tests/modis_tests/test_modis_l3.py | 23 +++++++ 2 files changed, 52 insertions(+), 33 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 2370609e7a..485fc1031f 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -49,30 +49,8 @@ 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] + # Get the grid resolution, name and other projection info + self.resolution = self._get_res() def _get_res(self): @@ -85,11 +63,12 @@ def _get_res(self): pos = gridname.rfind("_") + 1 pos2 = gridname.rfind("Deg") + # Initialise number of rows and columns # Some products don't have resolution listed. if pos < 0 or pos2 < 0: - self.resolution = 360. / self.ncols + return 360. / self.metadata["GridStructure"]["GRID_1"]["XDim"] else: - self.resolution = float(gridname[pos:pos2]) + return float(gridname[pos:pos2]) def available_datasets(self, configured_datasets=None): @@ -104,7 +83,7 @@ def available_datasets(self, configured_datasets=None): 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 + # 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: @@ -122,11 +101,9 @@ def available_datasets(self, configured_datasets=None): 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 + # skip variables that YAML had configured continue common = {"file_type": "modis_l3_cmg_hdf", "resolution": self.resolution, @@ -141,6 +118,20 @@ def get_dataset(self, dataset_id, dataset_info): 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"] + + # 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_area_def(self, dsid): """Get the area definition. @@ -150,12 +141,17 @@ def get_area_def(self, dsid): """ proj_param = "EPSG:4326" + # 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("gridded_modis", "A gridded L3 MODIS area", "longlat", proj_param, - self.ncols, - self.nrows, - self.area_extent) + ncols, + nrows, + self._get_area_extent()) return area diff --git a/satpy/tests/reader_tests/modis_tests/test_modis_l3.py b/satpy/tests/reader_tests/modis_tests/test_modis_l3.py index c7a0d3e9cf..3f6a9e8250 100644 --- a/satpy/tests/reader_tests/modis_tests/test_modis_l3.py +++ b/satpy/tests/reader_tests/modis_tests/test_modis_l3.py @@ -63,6 +63,29 @@ def test_scene_available_datasets(self, loadable, filename): assert len(available_datasets) > 0 assert loadable in available_datasets + from satpy.readers.modis_l3 import ModisL3GriddedHDFFileHandler + fh = ModisL3GriddedHDFFileHandler(filename[0], {}, {"file_type": "modis_l3_cmg_hdf"}) + configured_datasets = [[None, {"name": "none_ds", "file_type": "modis_l3_cmg_hdf"}], + [True, {"name": "true_ds", "file_type": "modis_l3_cmg_hdf"}], + [False, {"name": "false_ds", "file_type": "modis_l3_cmg_hdf"}], + [None, {"name": "other_ds", "file_type": "modis_l2_random"}]] + for status, mda in fh.available_datasets(configured_datasets): + if mda["name"] == "none_ds": + assert mda["file_type"] == "modis_l3_cmg_hdf" + assert status is False + elif mda["name"] == "true_ds": + assert mda["file_type"] == "modis_l3_cmg_hdf" + assert status + elif mda["name"] == "false_ds": + assert mda["file_type"] == "modis_l3_cmg_hdf" + assert status is False + elif mda["name"] == "other_ds": + assert mda["file_type"] == "modis_l2_random" + assert status is None + elif mda["name"] == loadable: + assert mda["file_type"] == "modis_l3_cmg_hdf" + assert status + def test_load_l3_dataset(self, modis_l3_nasa_mod09_file): """Load and check an L2 variable.""" scene = Scene(reader="modis_l3", filenames=modis_l3_nasa_mod09_file) From 3dab2c6587b92f98214ce4d204958074051f2a7b Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 17:20:46 +0000 Subject: [PATCH 24/30] Refactor MODIS tests. --- .../modis_tests/_modis_fixtures.py | 80 +++++++++---------- .../reader_tests/modis_tests/test_modis_l3.py | 4 +- 2 files changed, 38 insertions(+), 46 deletions(-) diff --git a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py index 0b79e00854..c702752b28 100644 --- a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py @@ -51,9 +51,6 @@ 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: @@ -229,21 +226,9 @@ def generate_imapp_filename(suffix): 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, - file_shortname: Optional[str] = None, - include_metadata: bool = True): + metadata_dict: Optional[dict] = {}): """Create a fake MODIS L1b HDF4 file with headers. Args: @@ -256,17 +241,27 @@ def create_hdfeos_test_file(filename: str, file_shortname: Short name of the file to be stored in global metadata attributes. Only used if ``include_metadata`` is ``True`` (default). - include_metadata: Include global metadata attributes (default: True). + metadata_dict: A dictionary of metadata to be added to the file. """ h = SD(filename, SDC.WRITE | SDC.CREATE) - 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, 'ArchiveMetadata.0', _create_header_metadata()) # noqa + if metadata_dict is not None and metadata_dict != {}: + # Check if we're dealing with an L3 file + if "l3_type" not in metadata_dict.keys(): + if "file_shortname" not in metadata_dict["file_shortname"].keys(): + raise ValueError("'file_shortname' is required when including metadata.") + # For L1 and L2 files we need to know the resolution + if "geo_resolution" not in metadata_dict.keys(): + raise ValueError("'geo_resolution' is required when including L1/L2 metadata.") + setattr(h, "StructMetadata.0", _create_struct_metadata(metadata_dict["geo_resolution"])) + setattr(h, "CoreMetadata.0", _create_core_metadata(metadata_dict["file_shortname"])) # noqa + else: + # For an L3 file, we just call the relevant metadata creator + setattr(h, "StructMetadata.0", _create_struct_metadata_cmg(metadata_dict["l3_type"])) + setattr(h, "CoreMetadata.0", _create_core_metadata(metadata_dict["l3_type"])) # noqa + + setattr(h, "ArchiveMetadata.0", _create_header_metadata()) # noqa for var_name, var_info in variable_infos.items(): _add_variable_to_file(h, var_name, var_info) @@ -339,13 +334,14 @@ def _create_struct_metadata(geo_resolution: int) -> str: return struct_metadata_header -def _create_struct_metadata_cmg(res) -> str: +def _create_struct_metadata_cmg(ftype: str) -> 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 + if ftype == "MOD09": + gridline = 'GridName="MOD09CMG"\n' + upleft = "UpperLeftPointMtrs=(-180000000.000000,90000000.000000)\n" + upright = "LowerRightMtrs=(180000000.000000,-90000000.000000)\n" + # Case of a MCD43 file + elif ftype == "MCD43C1": gridline = 'GridName="MCD_CMG_BRDF_0.05Deg"\n' upleft = "UpperLeftPointMtrs=(-180.000000,90.000000)\n" upright = "LowerRightMtrs=(180.000000,-90.000000)\n" @@ -381,7 +377,7 @@ def modis_l1b_nasa_mod021km_file(tmpdir_factory) -> list[str]: variable_infos.update(_get_visible_variable_info("EV_500_Aggr1km_RefSB", 1000, AVAILABLE_HKM_PRODUCT_NAMES)) variable_infos.update(_get_visible_variable_info("EV_250_Aggr1km_RefSB", 1000, AVAILABLE_QKM_PRODUCT_NAMES)) variable_infos.update(_get_emissive_variable_info("EV_1KM_Emissive", 1000, AVAILABLE_1KM_IR_PRODUCT_NAMES)) - create_hdfeos_test_file(full_path, variable_infos, geo_resolution=5000, file_shortname="MOD021KM") + create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 5000, "file_shortname": "MOD021KM"}) return [full_path] @@ -395,7 +391,7 @@ def modis_l1b_imapp_1000m_file(tmpdir_factory) -> list[str]: variable_infos.update(_get_visible_variable_info("EV_500_Aggr1km_RefSB", 1000, AVAILABLE_HKM_PRODUCT_NAMES)) variable_infos.update(_get_visible_variable_info("EV_250_Aggr1km_RefSB", 1000, AVAILABLE_QKM_PRODUCT_NAMES)) variable_infos.update(_get_emissive_variable_info("EV_1KM_Emissive", 1000, AVAILABLE_1KM_IR_PRODUCT_NAMES)) - create_hdfeos_test_file(full_path, variable_infos, geo_resolution=5000, file_shortname="MOD021KM") + create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 5000, "file_shortname": "MOD021KM"}) return [full_path] @@ -406,7 +402,7 @@ def modis_l1b_nasa_mod02hkm_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l1b").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 1000, include_angles=False) variable_infos.update(_get_visible_variable_info("EV_500_RefSB", 250, AVAILABLE_QKM_PRODUCT_NAMES)) - create_hdfeos_test_file(full_path, variable_infos, geo_resolution=1000, file_shortname="MOD02HKM") + create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 1000, "file_shortname": "MOD02HKM"}) return [full_path] @@ -417,7 +413,7 @@ def modis_l1b_nasa_mod02qkm_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l1b").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 1000, include_angles=False) variable_infos.update(_get_visible_variable_info("EV_250_RefSB", 250, AVAILABLE_QKM_PRODUCT_NAMES)) - create_hdfeos_test_file(full_path, variable_infos, geo_resolution=1000, file_shortname="MOD02QKM") + create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 1000, "file_shortname": "MOD02QKM"}) return [full_path] @@ -427,7 +423,7 @@ def modis_l1b_nasa_mod03_file(tmpdir_factory) -> list[str]: filename = generate_nasa_l1b_filename("MOD03") full_path = str(tmpdir_factory.mktemp("modis_l1b").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 1000, include_angles=True) - create_hdfeos_test_file(full_path, variable_infos, geo_resolution=1000, file_shortname="MOD03") + create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 1000, "file_shortname": "MOD03"}) return [full_path] @@ -437,7 +433,7 @@ def modis_l1b_imapp_geo_file(tmpdir_factory) -> list[str]: filename = generate_imapp_filename("geo") full_path = str(tmpdir_factory.mktemp("modis_l1b").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 1000, include_angles=True) - create_hdfeos_test_file(full_path, variable_infos, geo_resolution=1000, file_shortname="MOD03") + create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 1000, "file_shortname": "MOD03"}) return [full_path] @@ -595,7 +591,7 @@ def modis_l2_nasa_mod35_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l2").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 5000, include_angles=True) variable_infos.update(_get_cloud_mask_variable_info("Cloud_Mask", 1000)) - create_hdfeos_test_file(full_path, variable_infos, geo_resolution=5000, file_shortname="MOD35") + create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 5000, "file_shortname": "MOD35"}) return [full_path] @@ -605,12 +601,12 @@ def generate_nasa_l3_filename(prefix: str) -> str: 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): +def modis_l3_file(tmpdir_factory, f_prefix, var_name, 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) + create_hdfeos_test_file(full_path, variable_infos, {"l3_type": f_short}) return [full_path] @@ -620,7 +616,6 @@ def modis_l3_nasa_mod09_file(tmpdir_factory) -> list[str]: return modis_l3_file(tmpdir_factory, "MOD09CMG", "Coarse_Resolution_Surface_Reflectance_Band_2", - -999, "MOD09") @@ -630,7 +625,6 @@ def modis_l3_nasa_mod43_file(tmpdir_factory) -> list[str]: return modis_l3_file(tmpdir_factory, "MCD43C1", "BRDF_Albedo_Parameter1_Band2", - -9999, "MCD43C1") @@ -647,7 +641,7 @@ def modis_l2_nasa_mod06_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l2").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 5000, include_angles=True) variable_infos.update(_get_basic_variable_info("Surface_Pressure", 5000)) - create_hdfeos_test_file(full_path, variable_infos, geo_resolution=5000, file_shortname="MOD06") + create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 5000, "file_shortname": "MOD06"}) return [full_path] @@ -658,7 +652,7 @@ def modis_l2_imapp_snowmask_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l2").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 5000, include_angles=False) variable_infos.update(_get_basic_variable_info("Snow_Mask", 1000)) - create_hdfeos_test_file(full_path, variable_infos, include_metadata=False) + create_hdfeos_test_file(full_path, variable_infos, {}) return [full_path] @@ -675,7 +669,7 @@ def modis_l2_imapp_mask_byte1_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l2").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 5000, include_angles=False) variable_infos.update(_get_mask_byte1_variable_info()) - create_hdfeos_test_file(full_path, variable_infos, include_metadata=False) + create_hdfeos_test_file(full_path, variable_infos, {}) return [full_path] diff --git a/satpy/tests/reader_tests/modis_tests/test_modis_l3.py b/satpy/tests/reader_tests/modis_tests/test_modis_l3.py index 3f6a9e8250..de8ff682a1 100644 --- a/satpy/tests/reader_tests/modis_tests/test_modis_l3.py +++ b/satpy/tests/reader_tests/modis_tests/test_modis_l3.py @@ -27,8 +27,6 @@ from satpy import Scene, available_readers -from ._modis_fixtures import _shape_for_resolution - def _expected_area(): proj_param = "EPSG:4326" @@ -101,6 +99,6 @@ def test_load_l3_dataset(self, modis_l3_nasa_mod09_file): assert data_arr_comp.dtype == data_arr.dtype assert data_arr_comp.dtype == np.float32 - assert data_arr_comp.shape == _shape_for_resolution(-999) + assert data_arr_comp.shape == (3600, 7200) assert data_arr_comp.attrs.get("resolution") == 0.05 assert data_arr_comp.attrs.get("area") == _expected_area() From 1b136b9c19f15f121c4468b97cfcdf9ff4380662 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 17:21:58 +0000 Subject: [PATCH 25/30] Fix bug in MODIS tests. --- satpy/tests/reader_tests/modis_tests/_modis_fixtures.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py index c702752b28..7c20091d31 100644 --- a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py @@ -249,7 +249,7 @@ def create_hdfeos_test_file(filename: str, if metadata_dict is not None and metadata_dict != {}: # Check if we're dealing with an L3 file if "l3_type" not in metadata_dict.keys(): - if "file_shortname" not in metadata_dict["file_shortname"].keys(): + if "file_shortname" not in metadata_dict.keys(): raise ValueError("'file_shortname' is required when including metadata.") # For L1 and L2 files we need to know the resolution if "geo_resolution" not in metadata_dict.keys(): From d0774021726ab09b17c1fbdd4bde8e4650e9a034 Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 17:25:08 +0000 Subject: [PATCH 26/30] Remove mutable argument from modis tests. Co-authored-by: David Hoese --- satpy/tests/reader_tests/modis_tests/_modis_fixtures.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py index 7c20091d31..84ac7fc5ae 100644 --- a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py @@ -228,7 +228,7 @@ def generate_imapp_filename(suffix): def create_hdfeos_test_file(filename: str, variable_infos: dict, - metadata_dict: Optional[dict] = {}): + metadata_dict: Optional[dict] = None): """Create a fake MODIS L1b HDF4 file with headers. Args: From fe554dd1763de12226d53a5ec643623ac380c1ae Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 17:59:46 +0000 Subject: [PATCH 27/30] Additional refactoring of MODIS tests. --- .../modis_tests/_modis_fixtures.py | 98 ++++++++++++------- 1 file changed, 61 insertions(+), 37 deletions(-) diff --git a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py index 84ac7fc5ae..e792b70d89 100644 --- a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py @@ -19,7 +19,6 @@ from __future__ import annotations from datetime import datetime, timedelta -from typing import Optional import numpy as np import pytest @@ -228,40 +227,29 @@ def generate_imapp_filename(suffix): def create_hdfeos_test_file(filename: str, variable_infos: dict, - metadata_dict: Optional[dict] = None): + struct_meta: str = "", + core_meta: str = "", + archive_meta: str = "", + ) -> None: """Create a fake MODIS L1b HDF4 file with headers. Args: filename: Full path of filename to be created. variable_infos: Dictionary mapping HDF4 variable names to dictionary of variable information (see ``_add_variable_to_file``). - geo_resolution: Resolution of geolocation datasets to be stored in the - metadata strings stored in the global metadata attributes. Only - used if ``include_metadata`` is ``True`` (default). - file_shortname: Short name of the file to be stored in global metadata - attributes. Only used if ``include_metadata`` is ``True`` - (default). - metadata_dict: A dictionary of metadata to be added to the file. + struct_meta: Contents of the 'StructMetadata.0' header. + core_meta: Contents of the 'CoreMetadata.0' header. + archive_meta:Contents of the 'ArchiveMetadata.0' header. """ h = SD(filename, SDC.WRITE | SDC.CREATE) - if metadata_dict is not None and metadata_dict != {}: - # Check if we're dealing with an L3 file - if "l3_type" not in metadata_dict.keys(): - if "file_shortname" not in metadata_dict.keys(): - raise ValueError("'file_shortname' is required when including metadata.") - # For L1 and L2 files we need to know the resolution - if "geo_resolution" not in metadata_dict.keys(): - raise ValueError("'geo_resolution' is required when including L1/L2 metadata.") - setattr(h, "StructMetadata.0", _create_struct_metadata(metadata_dict["geo_resolution"])) - setattr(h, "CoreMetadata.0", _create_core_metadata(metadata_dict["file_shortname"])) # noqa - else: - # For an L3 file, we just call the relevant metadata creator - setattr(h, "StructMetadata.0", _create_struct_metadata_cmg(metadata_dict["l3_type"])) - setattr(h, "CoreMetadata.0", _create_core_metadata(metadata_dict["l3_type"])) # noqa - - setattr(h, "ArchiveMetadata.0", _create_header_metadata()) # noqa + if struct_meta != "": + setattr(h, "StructMetadata.0", struct_meta) + if core_meta != "": + setattr(h, "CoreMetadata.0", core_meta) + if archive_meta != "": + setattr(h, "ArchiveMetadata.0", archive_meta) for var_name, var_info in variable_infos.items(): _add_variable_to_file(h, var_name, var_info) @@ -341,7 +329,7 @@ def _create_struct_metadata_cmg(ftype: str) -> str: upleft = "UpperLeftPointMtrs=(-180000000.000000,90000000.000000)\n" upright = "LowerRightMtrs=(180000000.000000,-90000000.000000)\n" # Case of a MCD43 file - elif ftype == "MCD43C1": + else: gridline = 'GridName="MCD_CMG_BRDF_0.05Deg"\n' upleft = "UpperLeftPointMtrs=(-180.000000,90.000000)\n" upright = "LowerRightMtrs=(180.000000,-90.000000)\n" @@ -377,7 +365,11 @@ def modis_l1b_nasa_mod021km_file(tmpdir_factory) -> list[str]: variable_infos.update(_get_visible_variable_info("EV_500_Aggr1km_RefSB", 1000, AVAILABLE_HKM_PRODUCT_NAMES)) variable_infos.update(_get_visible_variable_info("EV_250_Aggr1km_RefSB", 1000, AVAILABLE_QKM_PRODUCT_NAMES)) variable_infos.update(_get_emissive_variable_info("EV_1KM_Emissive", 1000, AVAILABLE_1KM_IR_PRODUCT_NAMES)) - create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 5000, "file_shortname": "MOD021KM"}) + create_hdfeos_test_file(full_path, + variable_infos, + _create_struct_metadata(5000), + _create_core_metadata("MOD021KM"), + _create_header_metadata()) return [full_path] @@ -391,7 +383,11 @@ def modis_l1b_imapp_1000m_file(tmpdir_factory) -> list[str]: variable_infos.update(_get_visible_variable_info("EV_500_Aggr1km_RefSB", 1000, AVAILABLE_HKM_PRODUCT_NAMES)) variable_infos.update(_get_visible_variable_info("EV_250_Aggr1km_RefSB", 1000, AVAILABLE_QKM_PRODUCT_NAMES)) variable_infos.update(_get_emissive_variable_info("EV_1KM_Emissive", 1000, AVAILABLE_1KM_IR_PRODUCT_NAMES)) - create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 5000, "file_shortname": "MOD021KM"}) + create_hdfeos_test_file(full_path, + variable_infos, + _create_struct_metadata(5000), + _create_core_metadata("MOD021KM"), + _create_header_metadata()) return [full_path] @@ -402,7 +398,11 @@ def modis_l1b_nasa_mod02hkm_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l1b").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 1000, include_angles=False) variable_infos.update(_get_visible_variable_info("EV_500_RefSB", 250, AVAILABLE_QKM_PRODUCT_NAMES)) - create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 1000, "file_shortname": "MOD02HKM"}) + create_hdfeos_test_file(full_path, + variable_infos, + _create_struct_metadata(1000), + _create_core_metadata("MOD02HKM"), + _create_header_metadata()) return [full_path] @@ -413,7 +413,11 @@ def modis_l1b_nasa_mod02qkm_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l1b").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 1000, include_angles=False) variable_infos.update(_get_visible_variable_info("EV_250_RefSB", 250, AVAILABLE_QKM_PRODUCT_NAMES)) - create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 1000, "file_shortname": "MOD02QKM"}) + create_hdfeos_test_file(full_path, + variable_infos, + _create_struct_metadata(1000), + _create_core_metadata("MOD02QKM"), + _create_header_metadata()) return [full_path] @@ -423,7 +427,11 @@ def modis_l1b_nasa_mod03_file(tmpdir_factory) -> list[str]: filename = generate_nasa_l1b_filename("MOD03") full_path = str(tmpdir_factory.mktemp("modis_l1b").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 1000, include_angles=True) - create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 1000, "file_shortname": "MOD03"}) + create_hdfeos_test_file(full_path, + variable_infos, + _create_struct_metadata(1000), + _create_core_metadata("MOD03"), + _create_header_metadata()) return [full_path] @@ -433,7 +441,11 @@ def modis_l1b_imapp_geo_file(tmpdir_factory) -> list[str]: filename = generate_imapp_filename("geo") full_path = str(tmpdir_factory.mktemp("modis_l1b").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 1000, include_angles=True) - create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 1000, "file_shortname": "MOD03"}) + create_hdfeos_test_file(full_path, + variable_infos, + _create_struct_metadata(1000), + _create_core_metadata("MOD03"), + _create_header_metadata()) return [full_path] @@ -591,7 +603,11 @@ def modis_l2_nasa_mod35_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l2").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 5000, include_angles=True) variable_infos.update(_get_cloud_mask_variable_info("Cloud_Mask", 1000)) - create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 5000, "file_shortname": "MOD35"}) + create_hdfeos_test_file(full_path, + variable_infos, + _create_struct_metadata(5000), + _create_core_metadata("MOD35"), + _create_header_metadata()) return [full_path] @@ -606,7 +622,11 @@ def modis_l3_file(tmpdir_factory, f_prefix, var_name, f_short): 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, {"l3_type": f_short}) + create_hdfeos_test_file(full_path, + variable_infos, + _create_struct_metadata_cmg(f_short), + _create_core_metadata(f_short), + _create_header_metadata()) return [full_path] @@ -641,7 +661,11 @@ def modis_l2_nasa_mod06_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l2").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 5000, include_angles=True) variable_infos.update(_get_basic_variable_info("Surface_Pressure", 5000)) - create_hdfeos_test_file(full_path, variable_infos, {"geo_resolution": 5000, "file_shortname": "MOD06"}) + create_hdfeos_test_file(full_path, + variable_infos, + _create_struct_metadata(5000), + _create_core_metadata("MOD06"), + _create_header_metadata()) return [full_path] @@ -652,7 +676,7 @@ def modis_l2_imapp_snowmask_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l2").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 5000, include_angles=False) variable_infos.update(_get_basic_variable_info("Snow_Mask", 1000)) - create_hdfeos_test_file(full_path, variable_infos, {}) + create_hdfeos_test_file(full_path, variable_infos) return [full_path] @@ -669,7 +693,7 @@ def modis_l2_imapp_mask_byte1_file(tmpdir_factory) -> list[str]: full_path = str(tmpdir_factory.mktemp("modis_l2").join(filename)) variable_infos = _get_l1b_geo_variable_info(filename, 5000, include_angles=False) variable_infos.update(_get_mask_byte1_variable_info()) - create_hdfeos_test_file(full_path, variable_infos, {}) + create_hdfeos_test_file(full_path, variable_infos) return [full_path] From a211edd9ba77c169b8f5989521a1be42d879b02f Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 19:34:15 +0000 Subject: [PATCH 28/30] Refactor MODIS L3 reader and improve tests. --- satpy/readers/modis_l3.py | 48 ++++++++----------- .../modis_tests/_modis_fixtures.py | 13 ++--- 2 files changed, 27 insertions(+), 34 deletions(-) diff --git a/satpy/readers/modis_l3.py b/satpy/readers/modis_l3.py index 485fc1031f..29e0247fdc 100644 --- a/satpy/readers/modis_l3.py +++ b/satpy/readers/modis_l3.py @@ -44,33 +44,6 @@ 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) - - # Get the grid resolution, name and other projection info - self.resolution = self._get_res() - - - 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") - - # Initialise number of rows and columns - # Some products don't have resolution listed. - if pos < 0 or pos2 < 0: - return 360. / self.metadata["GridStructure"]["GRID_1"]["XDim"] - else: - return float(gridname[pos:pos2]) - - def available_datasets(self, configured_datasets=None): """Automatically determine datasets provided by this file.""" @@ -101,15 +74,34 @@ def available_datasets(self, configured_datasets=None): 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": "modis_l3_cmg_hdf", - "resolution": self.resolution, + "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 "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") + + # Initialise number of rows and columns + # Some products don't have resolution listed. + if pos < 0 or pos2 < 0: + return 360. / self.metadata["GridStructure"]["GRID_1"]["XDim"] + else: + return float(gridname[pos:pos2]) + def get_dataset(self, dataset_id, dataset_info): """Get DataArray for specified dataset.""" dataset_name = dataset_id["name"] diff --git a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py index e792b70d89..6dc4bf2d05 100644 --- a/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py +++ b/satpy/tests/reader_tests/modis_tests/_modis_fixtures.py @@ -19,6 +19,7 @@ from __future__ import annotations from datetime import datetime, timedelta +from typing import Optional import numpy as np import pytest @@ -227,9 +228,9 @@ def generate_imapp_filename(suffix): def create_hdfeos_test_file(filename: str, variable_infos: dict, - struct_meta: str = "", - core_meta: str = "", - archive_meta: str = "", + struct_meta: Optional[str] = None, + core_meta: Optional[str] = None, + archive_meta: Optional[str] = None, ) -> None: """Create a fake MODIS L1b HDF4 file with headers. @@ -244,11 +245,11 @@ def create_hdfeos_test_file(filename: str, """ h = SD(filename, SDC.WRITE | SDC.CREATE) - if struct_meta != "": + if struct_meta: setattr(h, "StructMetadata.0", struct_meta) - if core_meta != "": + if core_meta: setattr(h, "CoreMetadata.0", core_meta) - if archive_meta != "": + if archive_meta: setattr(h, "ArchiveMetadata.0", archive_meta) for var_name, var_info in variable_infos.items(): From 328816ec1dac80f1d3773ced4a0a7a2ab3b801cf Mon Sep 17 00:00:00 2001 From: Simon Proud Date: Wed, 15 Nov 2023 20:18:06 +0000 Subject: [PATCH 29/30] Update HDFEOS code to always use `maxsplit=1` when splitting attrs. --- satpy/readers/hdfeos_base.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/satpy/readers/hdfeos_base.py b/satpy/readers/hdfeos_base.py index 6e25fd40a8..37fe714435 100644 --- a/satpy/readers/hdfeos_base.py +++ b/satpy/readers/hdfeos_base.py @@ -148,10 +148,7 @@ def _read_mda(cls, lines, element=None): @classmethod def _split_line(cls, line, lines): - try: - key, val = line.split("=") - except ValueError: - key, val = line.split("=", maxsplit=1) + key, val = line.split("=", maxsplit=1) key = key.strip() val = val.strip() try: From 0a18d2da91893472fd802175161405157c3854e6 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 16 Nov 2023 09:29:53 -0600 Subject: [PATCH 30/30] Add test for VIIRS EDR available datasets fix --- satpy/tests/reader_tests/test_viirs_edr.py | 27 ++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/satpy/tests/reader_tests/test_viirs_edr.py b/satpy/tests/reader_tests/test_viirs_edr.py index 9b13f384e2..a6932520c0 100644 --- a/satpy/tests/reader_tests/test_viirs_edr.py +++ b/satpy/tests/reader_tests/test_viirs_edr.py @@ -285,6 +285,33 @@ def _create_fake_dataset(vars_dict: dict[str, xr.DataArray]) -> xr.Dataset: return ds +def test_available_datasets(aod_file): + """Test that available datasets doesn't claim non-filetype datasets. + + For example, if a YAML-configured dataset's file type is not loaded + then the available status is `None` and should remain `None`. This + means no file type knows what to do with this dataset. If it is + `False` then that means that a file type knows of the dataset, but + that the variable is not available in the file. In the below test + this isn't the case so the YAML-configured dataset should be + provided once and have a `None` availability. + + """ + from satpy.readers.viirs_edr import VIIRSJRRFileHandler + file_handler = VIIRSJRRFileHandler( + aod_file, + {"platform_shortname": "npp"}, + {"file_type": "jrr_aod"}, + ) + fake_yaml_datasets = [ + (None, {"file_key": "fake", "file_type": "fake_file", "name": "fake"}), + ] + available_datasets = list(file_handler.available_datasets(configured_datasets=fake_yaml_datasets)) + fake_availables = [avail_tuple for avail_tuple in available_datasets if avail_tuple[1]["name"] == "fake"] + assert len(fake_availables) == 1 + assert fake_availables[0][0] is None + + class TestVIIRSJRRReader: """Test the VIIRS JRR L2 reader."""