Merge pull request #3064 from mraspaud/feature-hrit-remote

Add remote reading for hrit seviri

Author: mraspaud

continuous_integration/environment.yaml

@@ -46,6 +46,7 @@ dependencies:
   - pytest
   - pytest-cov
   - fsspec
+  - universal_pathlib
   - botocore>=1.33
   - s3fs
   - python-geotiepoints

satpy/conftest.py

@@ -16,6 +16,7 @@
 # You should have received a copy of the GNU General Public License along with
 # satpy.  If not, see <http://www.gnu.org/licenses/>.
 """Pytest configuration and setup functions."""
+import pytest
 
 
 def pytest_configure(config):
@@ -28,3 +29,9 @@ def pytest_unconfigure(config):
     """Undo previous configurations."""
     from satpy import aux_download
     aux_download.RUNNING_TESTS = False
+
+
+@pytest.fixture(scope="session")
+def session_tmp_path(tmp_path_factory):
+    """Generate a single temp path to use for the entire session."""
+    return tmp_path_factory.mktemp("data")

satpy/readers/hrit_base.py

@@ -36,7 +36,6 @@
 from pyresample import geometry
 
 import satpy.readers.utils as utils
-from satpy.readers import FSFile
 from satpy.readers.eum_base import time_cds_short
 from satpy.readers.file_handlers import BaseFileHandler
 from satpy.readers.seviri_base import dec10216
@@ -88,14 +87,18 @@
                 }
 
 
-def decompress(infile):
+def decompress_file(infile) -> bytes:
     """Decompress an XRIT data file and return the decompressed buffer."""
-    from pyPublicDecompWT import xRITDecompress
-
     # decompress in-memory
     with open(infile, mode="rb") as fh:
-        xrit = xRITDecompress()
-        xrit.decompress(fh.read())
+        return decompress_buffer(fh.read())
+
+
+def decompress_buffer(buffer) -> bytes:
+    """Decompress buffer."""
+    from pyPublicDecompWT import xRITDecompress
+    xrit = xRITDecompress()
+    xrit.decompress(buffer)
 
     return xrit.data()
 
@@ -117,28 +120,31 @@ class HRITFileHandler(BaseFileHandler):
 
     def __init__(self, filename, filename_info, filetype_info, hdr_info):
         """Initialize the reader."""
-        super(HRITFileHandler, self).__init__(filename, filename_info,
-                                              filetype_info)
+        super().__init__(filename, filename_info, filetype_info)
 
         self.mda = {}
         self.hdr_info = hdr_info
         self._get_hd(self.hdr_info)
         self._start_time = filename_info["start_time"]
         self._end_time = self._start_time + dt.timedelta(minutes=15)
 
-    def _get_hd(self, hdr_info):
+    def _get_hd(self, hdr_info, verbose=False):
         """Open the file, read and get the basic file header info and set the mda dictionary."""
         hdr_map, variable_length_headers, text_headers = hdr_info
-
         with utils.generic_open(self.filename, mode="rb") as fp:
             total_header_length = 16
             while fp.tell() < total_header_length:
                 hdr_id = get_header_id(fp)
+                if verbose:
+                    print("hdr_id")  # noqa: T201
+                    print(f'np.void({hdr_id}, dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),')  # noqa: T201
                 the_type = hdr_map[hdr_id["hdr_id"]]
                 if the_type in variable_length_headers:
                     field_length = int((hdr_id["record_length"] - 3) /
                                        the_type.itemsize)
                     current_hdr = get_header_content(fp, the_type, field_length)
+                    if verbose:
+                        print(f"np.zeros(({field_length}, ), dtype={the_type}),")  # noqa: T201
                     key = variable_length_headers[the_type]
                     if key in self.mda:
                         if not isinstance(self.mda[key], list):
@@ -152,9 +158,13 @@ def _get_hd(self, hdr_info):
                     char = list(the_type.fields.values())[0][0].char
                     new_type = np.dtype(char + str(field_length))
                     current_hdr = get_header_content(fp, new_type)[0]
+                    if verbose:
+                        print(f'np.array({current_hdr}, dtype="{new_type}"),')  # noqa: T201
                     self.mda[text_headers[the_type]] = current_hdr
                 else:
                     current_hdr = get_header_content(fp, the_type)[0]
+                    if verbose:
+                        print(f"np.void({current_hdr}, dtype={the_type}),")  # noqa: T201
                     self.mda.update(
                         dict(zip(current_hdr.dtype.names, current_hdr)))
 
@@ -318,7 +328,7 @@ def _read_data_from_file(self):
         return self._read_data_from_disk()
 
     def _is_file_like(self):
-        return isinstance(self.filename, FSFile)
+        return not isinstance(self.filename, str)
 
     def _read_data_from_disk(self):
         # For reading the image data, unzip_context is faster than generic_open
@@ -327,7 +337,7 @@ def _read_data_from_disk(self):
 
             if self.compressed:
                 return np.frombuffer(
-                    decompress(fn),
+                    decompress_file(fn),
                     offset=self.offset,
                     dtype=dtype,
                     count=np.prod(shape)
@@ -344,12 +354,15 @@ def _read_file_like(self):
         # filename is likely to be a file-like object, already in memory
         dtype, shape = self._get_input_info()
         with utils.generic_open(self.filename, mode="rb") as fp:
+            decompressed_buffer = fp.read()
+            if self.compressed:
+                decompressed_buffer = decompress_buffer(decompressed_buffer)
             no_elements = np.prod(shape)
-            fp.seek(self.offset)
             return np.frombuffer(
-                fp.read(np.dtype(dtype).itemsize * no_elements),
+                decompressed_buffer,
                 dtype=np.dtype(dtype),
-                count=no_elements.item()
+                count=no_elements.item(),
+                offset=self.offset
             ).reshape(shape)
 
     def _get_input_info(self):

satpy/readers/seviri_l1b_hrit.py

@@ -312,7 +312,7 @@ class HRITMSGPrologueEpilogueBase(HRITFileHandler):
 
     def __init__(self, filename, filename_info, filetype_info, hdr_info):
         """Initialize the file handler for prologue and epilogue files."""
-        super(HRITMSGPrologueEpilogueBase, self).__init__(filename, filename_info, filetype_info, hdr_info)
+        super().__init__(filename, filename_info, filetype_info, hdr_info)
         self._reduced = None
 
     def _reduce(self, mda, max_size):
@@ -333,11 +333,11 @@ def __init__(self, filename, filename_info, filetype_info, calib_mode="nominal",
                  ext_calib_coefs=None, include_raw_metadata=False,
                  mda_max_array_size=None, fill_hrv=None, mask_bad_quality_scan_lines=None):
         """Initialize the reader."""
-        super(HRITMSGPrologueFileHandler, self).__init__(filename, filename_info,
-                                                         filetype_info,
-                                                         (msg_hdr_map,
-                                                          msg_variable_length_headers,
-                                                          msg_text_headers))
+        super().__init__(filename, filename_info,
+                         filetype_info,
+                         (msg_hdr_map,
+                          msg_variable_length_headers,
+                          msg_text_headers))
         self.prologue = {}
         self.read_prologue()
 

satpy/readers/utils.py

@@ -355,10 +355,10 @@ def generic_open(filename, *args, **kwargs):
             fp = filename.open(*args, **kwargs)
         except AttributeError:
             fp = open(filename, *args, **kwargs)
-
-    yield fp
-
-    fp.close()
+    try:
+        yield fp
+    finally:
+        fp.close()
 
 
 def fromfile(filename, dtype, count=1, offset=0):

satpy/tests/reader_tests/test_hrit_base.py

@@ -241,7 +241,7 @@ def test_read_band_filepath(self, stub_compressed_hrit_file):
         """Test reading a single band from a filepath."""
         filename = stub_compressed_hrit_file
 
-        with mock.patch("satpy.readers.hrit_base.decompress", side_effect=fake_decompress) as mock_decompress:
+        with mock.patch("satpy.readers.hrit_base.decompress_buffer", side_effect=fake_decompress) as mock_decompress:
             with mock.patch.object(HRITFileHandler, "_get_hd", side_effect=new_get_hd, autospec=True) as get_hd:
                 self.reader = HRITFileHandler(filename,
                                               {"platform_shortname": "MSG3",

satpy/tests/reader_tests/test_seviri_l1b_hrit.py

@@ -19,16 +19,21 @@
 """The HRIT msg reader tests package."""
 
 import datetime as dt
+import os
 import unittest
+import warnings
+import zipfile
 from unittest import mock
 
+import fsspec
 import numpy as np
 import pytest
 import xarray as xr
 from numpy import testing as npt
 from pyproj import CRS
 
 import satpy.tests.reader_tests.test_seviri_l1b_hrit_setup as setup
+from satpy.readers import FSFile
 from satpy.readers.seviri_l1b_hrit import HRITMSGEpilogueFileHandler, HRITMSGFileHandler, HRITMSGPrologueFileHandler
 from satpy.tests.reader_tests.test_seviri_base import ORBIT_POLYNOMIALS_INVALID
 from satpy.tests.reader_tests.test_seviri_l1b_calibration import TestFileHandlerCalibrationBase
@@ -502,3 +507,173 @@ def test_mask_bad_quality(self, file_handler):
         new_data[:, :] = np.nan
         expected = expected.copy(data=new_data)
         xr.testing.assert_equal(res, expected)
+
+
+@pytest.fixture(scope="session")
+def prologue_file(session_tmp_path, prologue_header_contents):
+    """Create a dummy prologue file."""
+    from satpy.readers.seviri_l1b_native_hdr import hrit_prologue
+    header = prologue_header_contents
+    contents = np.void(1, dtype=hrit_prologue)
+    contents["SatelliteStatus"]["SatelliteDefinition"]["SatelliteId"] = 324
+    return create_file(session_tmp_path / "prologue", header + [contents])
+
+
+@pytest.fixture(scope="session")
+def prologue_header_contents():
+    """Get the contents of the header."""
+    return [
+        # prime header
+        np.void((0, 16), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.void((128, 90, 3403688),
+                dtype=[("file_type", "u1"), ("total_header_length", ">u4"), ("data_field_length", ">u8")]),
+        # second header
+        np.void((4, 64), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.array(b"H-000-MSG4__-MSG4________-_________-PRO______-201802281500-__", dtype="|S61"),
+        # timestamp record
+        np.void((5, 10), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.void((64, (21973, 54911033)),
+                dtype=[("cds_p_field", "u1"), ("timestamp", [("Days", ">u2"), ("Milliseconds", ">u4")])])
+    ]
+
+
+@pytest.fixture(scope="session")
+def epilogue_file(session_tmp_path, epilogue_header_contents):
+    """Create a dummy epilogue file."""
+    from satpy.readers.seviri_l1b_native_hdr import hrit_epilogue
+    header = epilogue_header_contents
+    contents = np.void(1, dtype=hrit_epilogue)
+    return create_file(session_tmp_path / "epilogue", header + [contents])
+
+
+@pytest.fixture(scope="session")
+def epilogue_header_contents():
+    """Get the contents of the header."""
+    return [
+        np.void((0, 16), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.void((129, 90, 3042600),
+                dtype=[("file_type", "u1"), ("total_header_length", ">u4"), ("data_field_length", ">u8")]),
+        np.void((4, 64), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.array(b"H-000-MSG4__-MSG4________-_________-EPI______-201802281500-__", dtype="|S61"),
+        np.void((5, 10), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.void((64, (21973, 54911033)),
+                dtype=[("cds_p_field", "u1"), ("timestamp", [("Days", ">u2"), ("Milliseconds", ">u4")])]),
+    ]
+
+
+def create_file(filename, file_contents):
+    """Create an hrit file."""
+    with open(filename, "wb") as fh:
+        for array in file_contents:
+            array.tofile(fh)
+    return filename
+
+
+@pytest.fixture(scope="session")
+def segment_file(session_tmp_path):
+    """Create a segment_file."""
+    cols = 3712
+    lines = 464
+    bpp = 10
+    header = [
+        np.void((0, 16), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.void((0, 6198, 17223680), dtype=[("file_type", "u1"), ("total_header_length", ">u4"),
+                                            ("data_field_length", ">u8")]),
+        np.void((1, 9), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.void((bpp, cols, lines, 0), dtype=[("number_of_bits_per_pixel", "u1"), ("number_of_columns", ">u2"),
+                                           ("number_of_lines", ">u2"), ("compression_flag_for_data", "u1")]),
+        np.void((2, 51), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.void((b"GEOS(+000.0)                    ", -13642337, -13642337, 1856, 1856),
+                dtype=[("projection_name", "S32"),
+                       ("cfac", ">i4"), ("lfac", ">i4"),
+                       ("coff", ">i4"), ("loff", ">i4")]),
+        np.void((4, 64), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.array(b"H-000-MSG4__-MSG4________-VIS008___-000001___-201802281500-__", dtype="|S61"),
+        np.void((5, 10), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.void((64, (21973, 54911033)), dtype=[("cds_p_field", "u1"), ("timestamp", [("Days", ">u2"),
+                                                                                      ("Milliseconds", ">u4")])]),
+        np.void((128, 13), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.void((324, 2, 1, 1, 8, 0), dtype=[("GP_SC_ID", ">i2"), ("spectral_channel_id", "i1"),
+                                             ("segment_sequence_number", ">u2"),
+                                             ("planned_start_segment_number", ">u2"),
+                                             ("planned_end_segment_number", ">u2"),
+                                             ("data_field_representation", "i1")]),
+        np.void((129, 6035), dtype=[("hdr_id", "u1"), ("record_length", ">u2")]),
+        np.zeros((464, ), dtype=[("line_number_in_grid", ">i4"),
+                                 ("line_mean_acquisition", [("days", ">u2"), ("milliseconds", ">u4")]),
+                                 ("line_validity", "u1"), ("line_radiometric_quality", "u1"),
+                                 ("line_geometric_quality", "u1")]),
+        ]
+    contents = np.empty(cols * lines * bpp // 8, dtype="u1")
+
+    return create_file(session_tmp_path / "segment", header + [contents])
+
+
+def test_read_real_segment(prologue_file, epilogue_file, segment_file):
+    """Test reading an hrit segment."""
+    info = dict(start_time=dt.datetime(2018, 2, 28, 15, 0), service="")
+    prologue_fh = HRITMSGPrologueFileHandler(prologue_file, info, dict())
+    epilogue_fh = HRITMSGEpilogueFileHandler(epilogue_file, info, dict())
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", category=UserWarning, message="No orbit polynomial valid")
+        filehandler = HRITMSGFileHandler(segment_file, info, dict(), prologue_fh, epilogue_fh)
+    res = filehandler.get_dataset(dict(name="VIS008", calibration="counts"),
+                                  dict(units="", wavelength=0.8, standard_name="counts"))
+    res.compute()
+
+
+@pytest.fixture(scope="session")
+def compressed_seviri_hrit_files(session_tmp_path, prologue_file, epilogue_file, segment_file):
+    """Return the fsspec paths to the given seviri hrit files inside a zip file."""
+    zip_full_path = session_tmp_path / "test_seviri_hrit.zip"
+    with zipfile.ZipFile(zip_full_path, mode="w") as archive:
+        for filename in (prologue_file, epilogue_file, segment_file):
+            archive.write(filename, os.path.basename(filename))
+    return {hrit_file: f"zip://{hrit_file}::file://{zip_full_path.as_posix()}"
+            for hrit_file in ("prologue", "epilogue", "segment")}
+
+def test_read_real_segment_zipped(compressed_seviri_hrit_files):
+    """Test reading a remote hrit segment passed as FSFile."""
+    info = dict(start_time=dt.datetime(2018, 2, 28, 15, 0), service="")
+    prologue = FSFile(fsspec.open(compressed_seviri_hrit_files["prologue"]))
+    prologue_fh = HRITMSGPrologueFileHandler(prologue, info, dict())
+    epilogue = FSFile(fsspec.open(compressed_seviri_hrit_files["epilogue"]))
+    epilogue_fh = HRITMSGEpilogueFileHandler(epilogue, info, dict())
+    segment = FSFile(fsspec.open(compressed_seviri_hrit_files["segment"]))
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", category=UserWarning, message="No orbit polynomial valid")
+        filehandler = HRITMSGFileHandler(segment, info, dict(), prologue_fh, epilogue_fh)
+    res = filehandler.get_dataset(dict(name="VIS008", calibration="counts"),
+                                  dict(units="", wavelength=0.8, standard_name="counts"))
+    res.compute()
+
+
+def to_upath(fsfile):
+    """Convert FSFile instance to UPath."""
+    from upath import UPath
+    fsfile_fs = fsfile.fs.to_dict()
+    fsfile_fs.pop("cls")
+    path = UPath(os.fspath(fsfile), **fsfile_fs)
+    return path
+
+
+def test_read_real_segment_zipped_with_upath(compressed_seviri_hrit_files):
+    """Test reading a remote hrit segment passed as UPath."""
+    info = dict(start_time=dt.datetime(2018, 2, 28, 15, 0), service="")
+
+    prologue = FSFile(fsspec.open(compressed_seviri_hrit_files["prologue"]))
+    prologue = to_upath(prologue)
+    prologue_fh = HRITMSGPrologueFileHandler(prologue, info, dict())
+
+    epilogue = FSFile(fsspec.open(compressed_seviri_hrit_files["epilogue"]))
+    epilogue = to_upath(epilogue)
+    epilogue_fh = HRITMSGEpilogueFileHandler(epilogue, info, dict())
+
+    segment = FSFile(fsspec.open(compressed_seviri_hrit_files["segment"]))
+    segment = to_upath(segment)
+    with warnings.catch_warnings():
+        warnings.filterwarnings("ignore", category=UserWarning, message="No orbit polynomial valid")
+        filehandler = HRITMSGFileHandler(segment, info, dict(), prologue_fh, epilogue_fh)
+    res = filehandler.get_dataset(dict(name="VIS008", calibration="counts"),
+                                  dict(units="", wavelength=0.8, standard_name="counts"))
+    res.compute()