312 bug fix update geostationary readers to support multiple scan times

docs/source/releases/latest/312-bug-fix-update-geostationary-readers-to-support-multiple-scan-times.yaml

@@ -0,0 +1,22 @@
+bug fix:
+- description: |
+    *From GEOIPS#312: Bug Fix: Update geostationary readers to support multiple scan times*
+
+    Geostationary readers previously could only handle one scan time per run through.
+    This is a limiting factor to these readers and one or more scientists at CIRA needed
+    these readers to be able to handle multiple scan times per use. We've made the
+    corresponding updates to geostaionary readers that were requested could read
+    multiple scan times. If this PR goes well, it might be a good idea to support this
+    functionality on other readers in the future.
+  related-issue:
+    number: 312
+    repo_url: 'https://github.com/NRLMMD-GEOIPS/geoips/'
+  title: 'Update Geostationary readers to support multiple scan times'
+  files:
+    modified:
+      - geoips/interfaces/module_based/readers.py
+      - geoips/plugins/modules/readers/abi_netcdf.py
+      - geoips/plugins/modules/readers/ahi_hsd.py
+      - geoips/plugins/modules/readers/seviri_hrit.py
+      - tests/scripts/seviri.WV-Upper.unprojected_image.sh
+      - tests/outputs/seviri.WV-Upper.unprojected_image/20231211.080000.msg-2.seviri.WV-Upper.self_register.69p07.nesdisstar.10p0.png

geoips/interfaces/module_based/readers.py

@@ -3,6 +3,9 @@
 
 """Readers interface module."""
 
+from os.path import basename
+from xarray import Dataset
+
 from geoips.interfaces.base import BaseModuleInterface
 
 
@@ -19,5 +22,48 @@ class ReadersInterface(BaseModuleInterface):
         "standard": ["metadata_only", "chans", "area_def", "self_register"]
     }
 
+    def concatenate_metadata(self, all_metadata):
+        """Merge together metadata sourced from a list of files into one dictionary.
+
+        Where the structure of the merged metadata is a nested dictionary of metadata.
+
+        Ie. (xarray_obj has no data and is merely just a container for metadata):
+        {"METADATA": xobj.source_file_attributes: {fname: xobj, ..., "fnamex": xobj}}
+
+        Parameters
+        ----------
+        all_metadata: list of xarray.Datasets
+            - The incoming metadata from X number of files
+
+        Returns
+        -------
+        md: dict of xarray Datasets
+            - All metadata merged into a dictionary of xarray Datasets
+        """
+        md = {"METADATA": Dataset()}
+
+        for md_idx in range(len(all_metadata)):
+            # Set required attributes of the top-level metadata when this loop is
+            # started
+            if md_idx == 0:
+                md["METADATA"].attrs = all_metadata[md_idx].attrs
+                md["METADATA"].attrs["source_file_names"] = []
+                md["METADATA"].attrs["source_file_attributes"] = {}
+                md["METADATA"].attrs["source_file_datetimes"] = []
+                md["METADATA"].attrs["end_datetime"] = all_metadata[-1].end_datetime
+
+            # Add to optional attributes of the top-level metadata for each xobj
+            # provided
+            fname = all_metadata[md_idx].attrs["source_file_names"][0]
+            md["METADATA"].attrs["source_file_names"].append(basename(fname))
+            md["METADATA"].attrs["source_file_attributes"][basename(fname)] = (
+                all_metadata[md_idx]
+            )
+            md["METADATA"].attrs["source_file_datetimes"].append(
+                all_metadata[md_idx].start_datetime
+            )
+
+        return md
+
 
 readers = ReadersInterface()

geoips/plugins/modules/readers/abi_netcdf.py

@@ -12,6 +12,7 @@
 
 from scipy.ndimage import zoom
 
+from geoips.interfaces import readers
 from geoips.utils.context_managers import import_optional_dependencies
 from geoips.plugins.modules.readers.utils.geostationary_geolocation import (
     get_geolocation_cache_filename,
@@ -509,6 +510,100 @@ def call(fnames, metadata_only=False, chans=None, area_def=None, self_register=F
     """
     Read ABI NetCDF data from a list of filenames.
 
+    Parameters
+    ----------
+    fnames : list
+        * List of strings, full paths to files
+    metadata_only : bool, default=False
+        * Return before actually reading data if True
+    chans : list of str, default=None
+        * List of desired channels (skip unneeded variables as needed).
+        * Include all channels if None.
+    area_def : pyresample.AreaDefinition, default=None
+        * Specify region to read
+        * Read all data if None.
+    self_register : str or bool, default=False
+        * register all data to the specified dataset id (as specified in the
+          return dictionary keys).
+        * Read multiple resolutions of data if False.
+
+    Returns
+    -------
+    dict of xarray.Datasets
+        * dictionary of xarray.Dataset objects with required Variables and
+          Attributes.
+        * Dictionary keys can be any descriptive dataset ids.
+
+    See Also
+    --------
+    :ref:`xarray_standards`
+        Additional information regarding required attributes and variables
+        for GeoIPS-formatted xarray Datasets.
+    """
+    all_metadata = readers.concatenate_metadata(
+        [call_single_time([x], metadata_only=True)["METADATA"] for x in fnames]
+    )
+    if metadata_only:
+        # from IPython import embed as shell
+
+        # shell()
+        return all_metadata
+
+    start_times = [dt for dt in all_metadata["METADATA"].attrs["source_file_datetimes"]]
+    times = list(set(start_times))
+    import collections
+
+    ingested_xarrays = collections.defaultdict(list)
+    for time in times:
+        scan_time_files = [dt == time for dt in start_times]
+        data_dict = call_single_time(
+            np.array(fnames)[scan_time_files],
+            metadata_only=metadata_only,
+            chans=chans,
+            area_def=area_def,
+            self_register=self_register,
+        )
+        for (
+            dname,
+            dset,
+        ) in data_dict.items():
+            ingested_xarrays[dname].append(dset)
+
+    if len(times) == 1:
+        # No need to stack if we are only reading in one scan time
+        # This is likely temporary to maintain backwards compatibility
+
+        # This is not hit if we are provided multiple scan times.
+        return data_dict
+
+    import xarray
+
+    # merged_dset = xarray.Dataset()
+    # Now that we've ingested all scan times, stack along time dimension
+    dict_xarrays = {}
+    for dname, list_xarrays in ingested_xarrays.items():
+        if dname == "METADATA":
+            continue
+        merged_dset = xarray.concat(list_xarrays, dim="time_dim")
+        merged_dset.attrs["start_datetime"] = min(times)
+        merged_dset.attrs["end_datetime"] = max(times)
+        merged_dset = merged_dset.assign_coords({"time_dim": times})
+        dict_xarrays[dname] = merged_dset
+
+    metadata = all_metadata["METADATA"]
+    metadata.attrs["source_file_names"] = [os.path.basename(fname) for fname in fnames]
+    metadata.attrs["start_datetime"] = min(times)
+    metadata.attrs["end_datetime"] = max(times)
+    dict_xarrays["METADATA"] = metadata
+    return dict_xarrays
+
+
+def call_single_time(
+    fnames, metadata_only=False, chans=None, area_def=None, self_register=False
+):
+    """
+    Read ABI NetCDF data from a list of filenames.
+
     Parameters
     ----------
     fnames : list
@@ -629,6 +724,9 @@ def call(fnames, metadata_only=False, chans=None, area_def=None, self_register=F
     xarray_obj.attrs["end_datetime"] = edt
     xarray_obj.attrs["source_name"] = "abi"
     xarray_obj.attrs["data_provider"] = "noaa"
+    xarray_obj.attrs["source_file_names"] = [
+        os.path.basename(fname) for fname in fnames
+    ]
 
     # G16 -> goes-16
     xarray_obj.attrs["platform_name"] = highest_md["file_info"]["platform_ID"].replace(

geoips/plugins/modules/readers/ahi_hsd.py

@@ -15,6 +15,8 @@
 import xarray
 from scipy.ndimage import zoom
 
+# GeoIPS Libraries
+from geoips.interfaces import readers
 from geoips.utils.memusg import print_mem_usage
 from geoips.utils.context_managers import import_optional_dependencies
 from geoips.plugins.modules.readers.utils.geostationary_geolocation import (
@@ -946,6 +948,103 @@ def call(
     """
     Read AHI HSD data data from a list of filenames.
 
+    Parameters
+    ----------
+    fnames : list
+        * List of strings, full paths to files
+    metadata_only : bool, default=False
+        * Return before actually reading data if True
+    chans : list of str, default=None
+        * List of desired channels (skip unneeded variables as needed).
+        * Include all channels if None.
+    area_def : pyresample.AreaDefinition, default=None
+        * Specify region to read
+        * Read all data if None.
+    self_register : str or bool, default=False
+        * register all data to the specified dataset id (as specified in the
+          return dictionary keys).
+        * Read multiple resolutions of data if False.
+
+    Returns
+    -------
+    dict of xarray.Datasets
+        * dictionary of xarray.Dataset objects with required Variables and
+          Attributes.
+        * Dictionary keys can be any descriptive dataset ids.
+
+    See Also
+    --------
+    :ref:`xarray_standards`
+        Additional information regarding required attributes and variables
+        for GeoIPS-formatted xarray Datasets.
+    """
+    LOG.interactive("AHI reader test_arg: %s", test_arg)
+    all_metadata = readers.concatenate_metadata(
+        [call_single_time([x], metadata_only=True)["METADATA"] for x in fnames]
+    )
+    if metadata_only:
+        return all_metadata
+
+    start_times = [dt for dt in all_metadata["METADATA"].attrs["source_file_datetimes"]]
+    times = list(set(start_times))
+    import collections
+
+    ingested_xarrays = collections.defaultdict(list)
+    for time in times:
+        scan_time_files = [dt == time for dt in start_times]
+        data_dict = call_single_time(
+            np.array(fnames)[scan_time_files],
+            metadata_only=metadata_only,
+            chans=chans,
+            area_def=area_def,
+            self_register=self_register,
+        )
+        for (
+            dname,
+            dset,
+        ) in data_dict.items():
+            ingested_xarrays[dname].append(dset)
+
+    if len(times) == 1:
+        # No need to stack if we are only reading in one scan time
+        # This is likely temporary to maintain backwards compatibility
+
+        # This is not hit if we are provided multiple scan times.
+        return data_dict
+
+    import xarray
+
+    # merged_dset = xarray.Dataset()
+    # Now that we've ingested all scan times, stack along time dimension
+    dict_xarrays = {}
+    for dname, list_xarrays in ingested_xarrays.items():
+        if dname == "METADATA":
+            continue
+        merged_dset = xarray.concat(list_xarrays, dim="time_dim")
+        merged_dset.attrs["start_datetime"] = min(times)
+        merged_dset.attrs["end_datetime"] = max(times)
+        merged_dset = merged_dset.assign_coords({"time_dim": times})
+        dict_xarrays[dname] = merged_dset
+
+    metadata = all_metadata["METADATA"]
+    metadata.attrs["source_file_names"] = [os.path.basename(fname) for fname in fnames]
+    metadata.attrs["start_datetime"] = min(times)
+    metadata.attrs["end_datetime"] = max(times)
+    dict_xarrays["METADATA"] = metadata
+    return dict_xarrays
+
+
+def call_single_time(
+    fnames,
+    metadata_only=False,
+    chans=None,
+    area_def=None,
+    self_register=False,
+    test_arg="AHI Default Test Arg",
+):
+    """
+    Read AHI HSD data data from a list of filenames.
+
     Parameters
     ----------
     fnames : list
@@ -1097,6 +1196,9 @@ def call(
     xarray_obj.attrs["data_provider"] = "jma"
     xarray_obj.attrs["platform_name"] = highest_md["block_01"]["satellite_name"].lower()
     xarray_obj.attrs["area_definition"] = area_def
+    xarray_obj.attrs["source_file_names"] = [
+        os.path.basename(fname) for fname in fnames
+    ]
 
     # If metadata_only requested, return here.
     if metadata_only: