From 8e904fbc71d2ddf4231500edcd0df89bbd8c3daf Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Thu, 27 Apr 2023 16:43:08 +0200 Subject: [PATCH 01/23] Added file defining IASI L1C format Source: https://gitlab.eumetsat.int/open-source/data-tailor/-/blob/d1115baf74d6ce9eb41189264d73fee175826383/epct/epct/resources/eps_iasil1c_9a_custom.xml Apache License 2.0 --- satpy/etc/eps_iasil1c_9a_custom.xml | 687 ++++++++++++++++++++++++++++ 1 file changed, 687 insertions(+) create mode 100644 satpy/etc/eps_iasil1c_9a_custom.xml diff --git a/satpy/etc/eps_iasil1c_9a_custom.xml b/satpy/etc/eps_iasil1c_9a_custom.xml new file mode 100644 index 0000000000..769a5ebf2e --- /dev/null +++ b/satpy/etc/eps_iasil1c_9a_custom.xml @@ -0,0 +1,687 @@ + + + + + + + + + + 110 + current + + + 9a + + + + EPS IASI Level 1C Format + + + + + + IASI_xxx_1C_*Z* + + + + + + + + + + + + + + + + + + + + + + + + Geolocation + IASI-1C Location of pixel centre in geodetic coordinates for sounder pixels + mdr-1c[].GGeoSondLoc[][][1] + mdr-1c[].GGeoSondLoc[][][0] + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Detailed quality flag + Detailed quality flag for the system + mdr-1c[].GQisFlagQualDetailed[][] + Geolocation + + + + GEUM AVHRR 1B cloud fraction + Cloud fraction in IASI FOV from AVHRR 1B in IASI FOV + mdr-1c[].GEUMAvhrr1BCldFrac[][] + Geolocation + + + + GEUM AVHRR 1B land fraction + Land and Coast fraction in IASI FOV from AVHRR 1B + mdr-1c[].GEUMAvhrr1BLandFrac[][] + Geolocation + + + + GEUM AVHRR 1B quality + Quality indicator for cloud and Land and Coast fraction + mdr-1c[].GEUMAvhrr1BQual[][] + Geolocation + + + + Spectra measurements + Level 1C spectra + mdr-1c[].GS1cSpect[][][] + Geolocation + + + + Azimuth angles + Azimuth angle for each sounder pixel + mdr-1c[].GGeoSondAnglesMETOP[][][1] + Geolocation + + + + Zenith angles + Zenith angle for each sounder pixel + mdr-1c[].GGeoSondAnglesMETOP[][][0] + Geolocation + + + + Solar azimuth angles + Solar azimuth angle at the surface for each sounder pixel + mdr-1c[].GGeoSondAnglesSUN[][][1] + Geolocation + + + + Solar zenith angles + Solar zenith angle at the surface for each sounder pixel + mdr-1c[].GGeoSondAnglesSUN[][][0] + Geolocation + + + eps-product + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From cc5e746cce5446d99ba12ba15f96e91302aeb922 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Fri, 28 Apr 2023 17:48:31 +0200 Subject: [PATCH 02/23] Refactor EPS reader, fix IASI EPS definition file MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Started refactoring the AVHRR L1b EPS reader in preparation for the IASI L2 EPS reader. Replaced the wrong XML format definition by the correct one. Source for the EPS IASI L2 format definition is eugene. After installing eugene with ``mamba install -c eumetsat eugene``, I found the definition file at ``envs/py311/share/eugene/formats/eps_iasil2_9.0.xml``. --- satpy/etc/eps_iasil1c_9a_custom.xml | 687 ----------------- satpy/etc/eps_iasil2_9.0.xml | 909 +++++++++++++++++++++++ satpy/etc/readers/iasi_l2_eps.yaml | 22 + satpy/readers/eps_base.py | 132 ++++ satpy/readers/eps_l1b.py | 108 +-- satpy/readers/iasi_l2_eps.py | 35 + satpy/tests/reader_tests/test_eps_l1b.py | 6 +- 7 files changed, 1118 insertions(+), 781 deletions(-) delete mode 100644 satpy/etc/eps_iasil1c_9a_custom.xml create mode 100644 satpy/etc/eps_iasil2_9.0.xml create mode 100644 satpy/etc/readers/iasi_l2_eps.yaml create mode 100644 satpy/readers/eps_base.py create mode 100644 satpy/readers/iasi_l2_eps.py diff --git a/satpy/etc/eps_iasil1c_9a_custom.xml b/satpy/etc/eps_iasil1c_9a_custom.xml deleted file mode 100644 index 769a5ebf2e..0000000000 --- a/satpy/etc/eps_iasil1c_9a_custom.xml +++ /dev/null @@ -1,687 +0,0 @@ - - - - - - - - - - 110 - current - - - 9a - - - - EPS IASI Level 1C Format - - - - - - IASI_xxx_1C_*Z* - - - - - - - - - - - - - - - - - - - - - - - - Geolocation - IASI-1C Location of pixel centre in geodetic coordinates for sounder pixels - mdr-1c[].GGeoSondLoc[][][1] - mdr-1c[].GGeoSondLoc[][][0] - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Detailed quality flag - Detailed quality flag for the system - mdr-1c[].GQisFlagQualDetailed[][] - Geolocation - - - - GEUM AVHRR 1B cloud fraction - Cloud fraction in IASI FOV from AVHRR 1B in IASI FOV - mdr-1c[].GEUMAvhrr1BCldFrac[][] - Geolocation - - - - GEUM AVHRR 1B land fraction - Land and Coast fraction in IASI FOV from AVHRR 1B - mdr-1c[].GEUMAvhrr1BLandFrac[][] - Geolocation - - - - GEUM AVHRR 1B quality - Quality indicator for cloud and Land and Coast fraction - mdr-1c[].GEUMAvhrr1BQual[][] - Geolocation - - - - Spectra measurements - Level 1C spectra - mdr-1c[].GS1cSpect[][][] - Geolocation - - - - Azimuth angles - Azimuth angle for each sounder pixel - mdr-1c[].GGeoSondAnglesMETOP[][][1] - Geolocation - - - - Zenith angles - Zenith angle for each sounder pixel - mdr-1c[].GGeoSondAnglesMETOP[][][0] - Geolocation - - - - Solar azimuth angles - Solar azimuth angle at the surface for each sounder pixel - mdr-1c[].GGeoSondAnglesSUN[][][1] - Geolocation - - - - Solar zenith angles - Solar zenith angle at the surface for each sounder pixel - mdr-1c[].GGeoSondAnglesSUN[][][0] - Geolocation - - - eps-product - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/satpy/etc/eps_iasil2_9.0.xml b/satpy/etc/eps_iasil2_9.0.xml new file mode 100644 index 0000000000..e45415785e --- /dev/null +++ b/satpy/etc/eps_iasil2_9.0.xml @@ -0,0 +1,909 @@ + + + + + + + + 11 + 0 + current + + + 9 + 0 + + + EPS IASI Level 2 Format + + This IASI 2 description was generated using the IASI PFS Excel document Issue 6 Revision 3 + + + IASI_???_02_*Z* + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Geolocation + Geolocation Coverage (Latitude, Longitude) + mdr[].EARTH_LOCATION[][0] + mdr[].EARTH_LOCATION[][1] + + + + PRESSURE + + + Atmospheric Temperature + Atmospheric Temperature + mdr[].ATMOSPHERIC_TEMPERATURE[][$PRESSURE] + Geolocation + + + + PRESSURE + + + Water Vapour + Atmospheric Water Vapour + mdr[].ATMOSPHERIC_WATER_VAPOUR[][$PRESSURE] + Geolocation + + + + PRESSURE + + + Ozone + Atmospheric Ozone + mdr[].ATMOSPHERIC_OZONE[][$PRESSURE] + Geolocation + + + Integrated Ozone + Integrated Ozone + mdr[].INTEGRATED_OZONE[] + Geolocation + + + Surface Temperature 1 + Surface Temperature 1 + mdr[].SURFACE_TEMPERATURE[][0] + Geolocation + + + Surface Temperature 2 + Surface Temperature 2 + mdr[].SURFACE_TEMPERATURE[][1] + Geolocation + + + Integrated N20 + Integrated N20 + mdr[].INEGRATED_N2O[] + Geolocation + + + Integrated CO + Integrated CO + mdr[].INTEGRATED_CO[] + Geolocation + + + Integrated CH4 + Integrated CH4 + mdr[].INTEGRATED_CH4[] + Geolocation + + + Integrated CO2 + Integrated CO2 + mdr[].INTEGRATED_CO2[] + Geolocation + + + + WAVELENGTH + + + Surface Emissivity + Surface Emissivity + mdr[].SURFACE_EMISSIVITY[][$WAVELENGTH] + Geolocation + + + Number of cloud formations + Number of cloud formations + mdr[].NUMBER_CLOUD_FORMATIONS[] + Geolocation + + + + CLOUD + + + Fractional cloud cover + Fractional cloud cover + mdr[].FRACTIONAL_CLOUD_COVER[][$CLOUD] + Geolocation + + + + CLOUD + + + Cloud top temperature + Cloud top temperature + mdr[].CLOUD_TOP_TEMPERATURE[][$CLOUD] + Geolocation + + + + CLOUD + + + Cloud top pressure + Cloud top pressure + mdr[].CLOUD_TOP_PRESSURE[][$CLOUD] + Geolocation + + + + CLOUD + + + Cloud phase + Cloud phase (0=no cloud, 1=liquid, 2=ice, 3=mixed) + mdr[].CLOUD_PHASE[][$CLOUD] + Geolocation + + + FLG_AVHAVL + Completness flag for AVHRR products + mdr[].FLG_AVHAVL[] + Geolocation + + + FLG_CHNSEL + Channel selection flag + mdr[].FLG_CHNSEL[] + Geolocation + + + FLG_CLDAVH + Cloud formations anaylsed in AVHRR flag + mdr[].FLG_CLDAVH[] + Geolocation + + + FLG_CLDPHA + Cloud phase flag + mdr[].FLG_CLDPHA[] + Geolocation + + + FLG_DAYNIT + Flag discriminating between day and night + mdr[].FLG_DAYNIT[] + Geolocation + + + FLG_IASICLR + IASI-l2 IFOV clear, partly cloudy or cloudy flag + mdr[].FLG_IASICLR[] + Geolocation + + + FLG_ITCONV + Convergence of the iterative retrieval flag + mdr[].FLG_ITCONV[] + Geolocation + + + FLG_ITRBOU + Validation of iterated state vector flag + mdr[].FLG_ITRBOU[] + Geolocation + + + FLG_LANSEA + Flag specifying surface type + mdr[].FLG_LANSEA[] + Geolocation + + + FLG_NUMIT + Number of iterations used for retrieval + mdr[].FLG_NUMIT[] + Geolocation + + + FLG_NWPBAD + Validation flag of NWP forecast + mdr[].FLG_NWPBAD[] + Geolocation + + + FLG_QUAL + Quality and completness of the retrieval + mdr[].FLG_QUAL[] + Geolocation + + + FLG_RESID + Retrieval residual acceptance flag + mdr[].FLG_RESID[] + Geolocation + + + FLG_SATMAN + Flag indicating satellite manoeuvre + mdr[].FLG_SATMAN[] + Geolocation + + + FLG_SELBAC + Selection of background state flag + mdr[].FLG_SELBAC[] + Geolocation + + + FLG_SUNGLNT + Flag identifying of sun glint + mdr[].FLG_SUNGLNT[] + Geolocation + + + FLG_SUPADI + Flag indicating super-adiabatic conditions in final retrieval + mdr[].FLG_SUPADI[] + Geolocation + + + FLG_SUPSAT + Flag indicating super-saturation in final retrieval + mdr[].FLG_SUPSAT[] + Geolocation + + + FLG_THICIR + Thin cirrus cloud test + mdr[].FLG_THICIR[] + Geolocation + + + FLG_THICOR + Thin cirrus has been corrected for + mdr[].FLG_THICOR[] + Geolocation + + + FLG_VARCLR + Cloud clearing by variational analysis + mdr[].FLG_VARCLR[] + Geolocation + + eps-product + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/satpy/etc/readers/iasi_l2_eps.yaml b/satpy/etc/readers/iasi_l2_eps.yaml new file mode 100644 index 0000000000..98e419d52a --- /dev/null +++ b/satpy/etc/readers/iasi_l2_eps.yaml @@ -0,0 +1,22 @@ +reader: + name: iasi_l2_eps + short_name: IASI l2 EPS + long_name: IASI Combined Sounding Products - Metop + description: > + Reader for IASI Combined Sounding Products - Metop files. + Data and documentation are available + from https://data.eumetsat.int/product/EO:EUM:DAT:METOP:IASSND02 and + from the EUMETSAT Data Store under ID EO:EUM:DAT:METOP:IASSND02. + status: Alpha + supports_fsspec: False + reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader + sensors: [iasi] + default_datasets: + +file_types: + iasi_l2_eps: + file_reader: !!python/name:satpy.readers.iasi_l2_eps.EPSIASIFile + # see https://www.eumetsat.int/media/40048 §7 for file pattern + file_patterns: + - IASI_SND_02_{platform_short_name}_{start_time:%Y%m%d%H%M%SZ}_{end_time:%Y%m%d%H%M%SZ}_{processing_mode}_{disposition_mode}_{creation_time:%Y%m%d%H%M%SZ} + - IASI_SND_02_{platform_short_name}_{start_time:%Y%m%d%H%M%SZ}_{end_time:%Y%m%d%H%M%SZ}_{processing_mode}_{disposition_mode}_{creation_time:%Y%m%d%H%M%SZ}.nat diff --git a/satpy/readers/eps_base.py b/satpy/readers/eps_base.py new file mode 100644 index 0000000000..b00e66a272 --- /dev/null +++ b/satpy/readers/eps_base.py @@ -0,0 +1,132 @@ +# 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 . + +"""Reader for eps level 1b data. Uses xml files as a format description.""" + +import logging + +import dask.array as da +import numpy as np +import xarray as xr + +from satpy._config import get_config_path +from satpy.readers.file_handlers import BaseFileHandler +from satpy.readers.xmlformat import XMLFormat +from satpy.utils import get_legacy_chunk_size + +logger = logging.getLogger(__name__) + +CHUNK_SIZE = get_legacy_chunk_size() + +C1 = 1.191062e-05 # mW/(m2*sr*cm-4) +C2 = 1.4387863 # K/cm-1 + + +record_class = ["Reserved", "mphr", "sphr", + "ipr", "geadr", "giadr", + "veadr", "viadr", "mdr"] + + +def read_records(filename, format_definition): + """Read ``filename`` without scaling it afterwards. + + Args: + filename (str): Filename to read. + format_definition (str): XML file containing the format definition. + Must be findable in the Satpy config path. + """ + format_fn = get_config_path(format_definition) + form = XMLFormat(format_fn) + + grh_dtype = np.dtype([("record_class", "|i1"), + ("INSTRUMENT_GROUP", "|i1"), + ("RECORD_SUBCLASS", "|i1"), + ("RECORD_SUBCLASS_VERSION", "|i1"), + ("RECORD_SIZE", ">u4"), + ("RECORD_START_TIME", "S6"), + ("RECORD_STOP_TIME", "S6")]) + + max_lines = np.floor((CHUNK_SIZE ** 2) / 2048) + + dtypes = [] + cnt = 0 + counts = [] + classes = [] + prev = None + with open(filename, "rb") as fdes: + while True: + grh = np.fromfile(fdes, grh_dtype, 1) + if grh.size == 0: + break + rec_class = record_class[int(grh["record_class"])] + sub_class = grh["RECORD_SUBCLASS"][0] + + expected_size = int(grh["RECORD_SIZE"]) + bare_size = expected_size - grh_dtype.itemsize + try: + the_type = form.dtype((rec_class, sub_class)) + # the_descr = grh_dtype.descr + the_type.descr + except KeyError: + the_type = np.dtype([('unknown', 'V%d' % bare_size)]) + the_descr = grh_dtype.descr + the_type.descr + the_type = np.dtype(the_descr) + if the_type.itemsize < expected_size: + padding = [('unknown%d' % cnt, 'V%d' % (expected_size - the_type.itemsize))] + cnt += 1 + the_descr += padding + new_dtype = np.dtype(the_descr) + key = (rec_class, sub_class) + if key == prev: + counts[-1] += 1 + else: + dtypes.append(new_dtype) + counts.append(1) + classes.append(key) + prev = key + fdes.seek(expected_size - grh_dtype.itemsize, 1) + + sections = {} + offset = 0 + for dtype, count, rec_class in zip(dtypes, counts, classes): + fdes.seek(offset) + if rec_class == ('mdr', 2): + record = da.from_array(np.memmap(fdes, mode='r', dtype=dtype, shape=count, offset=offset), + chunks=(max_lines,)) + else: + record = np.fromfile(fdes, dtype=dtype, count=count) + offset += dtype.itemsize * count + if rec_class in sections: + logger.debug('Multiple records for ', str(rec_class)) + sections[rec_class] = np.hstack((sections[rec_class], record)) + else: + sections[rec_class] = record + + return sections, form + + +def create_xarray(arr): + """Create xarray with correct dimensions.""" + res = arr + res = xr.DataArray(res, dims=['y', 'x']) + return res + + +class EPSFile(BaseFileHandler): + """Base class for EPS level 1b readers.""" + + spacecrafts = {"M01": "Metop-B", + "M02": "Metop-A", + "M03": "Metop-C", } diff --git a/satpy/readers/eps_l1b.py b/satpy/readers/eps_l1b.py index 1cc098a612..d58f29771e 100644 --- a/satpy/readers/eps_l1b.py +++ b/satpy/readers/eps_l1b.py @@ -1,6 +1,4 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# Copyright (c) 2017-2020 Satpy developers +# Copyright (c) 2017-2023 Satpy developers # # This file is part of satpy. # @@ -16,7 +14,7 @@ # You should have received a copy of the GNU General Public License along with # satpy. If not, see . -"""Reader for eps level 1b data. Uses xml files as a format description.""" +"""Reader for EPS level 1b data. Uses xml files as a format description.""" import functools import logging @@ -28,9 +26,7 @@ from pyresample.geometry import SwathDefinition from satpy._compat import cached_property -from satpy._config import get_config_path -from satpy.readers.file_handlers import BaseFileHandler -from satpy.readers.xmlformat import XMLFormat +from satpy.readers.eps_base import EPSFile, XMLFormat, read_records, record_class # noqa from satpy.utils import get_legacy_chunk_size logger = logging.getLogger(__name__) @@ -51,82 +47,6 @@ def radiance_to_refl(arr, solar_flux): return arr * np.pi * 100.0 / solar_flux -record_class = ["Reserved", "mphr", "sphr", - "ipr", "geadr", "giadr", - "veadr", "viadr", "mdr"] - - -def read_records(filename): - """Read *filename* without scaling it afterwards.""" - format_fn = get_config_path("eps_avhrrl1b_6.5.xml") - form = XMLFormat(format_fn) - - grh_dtype = np.dtype([("record_class", "|i1"), - ("INSTRUMENT_GROUP", "|i1"), - ("RECORD_SUBCLASS", "|i1"), - ("RECORD_SUBCLASS_VERSION", "|i1"), - ("RECORD_SIZE", ">u4"), - ("RECORD_START_TIME", "S6"), - ("RECORD_STOP_TIME", "S6")]) - - max_lines = np.floor((CHUNK_SIZE ** 2) / 2048) - - dtypes = [] - cnt = 0 - counts = [] - classes = [] - prev = None - with open(filename, "rb") as fdes: - while True: - grh = np.fromfile(fdes, grh_dtype, 1) - if grh.size == 0: - break - rec_class = record_class[int(grh["record_class"])] - sub_class = grh["RECORD_SUBCLASS"][0] - - expected_size = int(grh["RECORD_SIZE"]) - bare_size = expected_size - grh_dtype.itemsize - try: - the_type = form.dtype((rec_class, sub_class)) - # the_descr = grh_dtype.descr + the_type.descr - except KeyError: - the_type = np.dtype([('unknown', 'V%d' % bare_size)]) - the_descr = grh_dtype.descr + the_type.descr - the_type = np.dtype(the_descr) - if the_type.itemsize < expected_size: - padding = [('unknown%d' % cnt, 'V%d' % (expected_size - the_type.itemsize))] - cnt += 1 - the_descr += padding - new_dtype = np.dtype(the_descr) - key = (rec_class, sub_class) - if key == prev: - counts[-1] += 1 - else: - dtypes.append(new_dtype) - counts.append(1) - classes.append(key) - prev = key - fdes.seek(expected_size - grh_dtype.itemsize, 1) - - sections = {} - offset = 0 - for dtype, count, rec_class in zip(dtypes, counts, classes): - fdes.seek(offset) - if rec_class == ('mdr', 2): - record = da.from_array(np.memmap(fdes, mode='r', dtype=dtype, shape=count, offset=offset), - chunks=(max_lines,)) - else: - record = np.fromfile(fdes, dtype=dtype, count=count) - offset += dtype.itemsize * count - if rec_class in sections: - logger.debug('Multiple records for ', str(rec_class)) - sections[rec_class] = np.hstack((sections[rec_class], record)) - else: - sections[rec_class] = record - - return sections, form - - def create_xarray(arr): """Create xarray with correct dimensions.""" res = arr @@ -134,12 +54,8 @@ def create_xarray(arr): return res -class EPSAVHRRFile(BaseFileHandler): - """Eps level 1b reader for AVHRR data.""" - - spacecrafts = {"M01": "Metop-B", - "M02": "Metop-A", - "M03": "Metop-C", } +class EPSAVHRRFile(EPSFile): + """EPS level 1b reader for AVHRR data.""" sensors = {"AVHR": "avhrr-3"} @@ -167,7 +83,7 @@ def __init__(self, filename, filename_info, filetype_info): def _read_all(self): logger.debug("Reading %s", self.filename) - self.sections, self.form = read_records(self.filename) + self.sections, self.form = read_records(self.filename, "eps_avhrrl1b_6.5.xml") self.scanlines = self['TOTAL_MDR'] if self.scanlines != len(self.sections[('mdr', 2)]): logger.warning("Number of declared records doesn't match number of scanlines in the file.") @@ -398,3 +314,15 @@ def end_time(self): """Get end time.""" # return datetime.strptime(self["SENSING_END"], "%Y%m%d%H%M%SZ") return self._end_time + + +class EPSIASIFile(EPSFile): + """EPS level 2 reader for IASI data. + + Reader for the IASI Level 2 combined sounding products in native format. + + Overview of the data including links to the product user guide, product format + specification, validation reports, and other documents, can be found at the + EUMETSAT Data Services at https://data.eumetsat.int/product/EO:EUM:DAT:METOP:IASSND02 + + """ diff --git a/satpy/readers/iasi_l2_eps.py b/satpy/readers/iasi_l2_eps.py new file mode 100644 index 0000000000..466f37a79e --- /dev/null +++ b/satpy/readers/iasi_l2_eps.py @@ -0,0 +1,35 @@ +# 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 . + +"""Reader for EPS level 2 IASI data. Uses xml files as a format description.""" + +from .eps_base import EPSFile + + +class EPSIASIFile(EPSFile): + """EPS level 2 reader for IASI data. + + Reader for the IASI Level 2 combined sounding products in native format. + + Overview of the data including links to the product user guide, product format + specification, validation reports, and other documents, can be found at the + EUMETSAT Data Services at https://data.eumetsat.int/product/EO:EUM:DAT:METOP:IASSND02 + + """ + + def __init__(self, filename, filename_info, filetype_info): + """Initialise Filehandler.""" + raise NotImplementedError() diff --git a/satpy/tests/reader_tests/test_eps_l1b.py b/satpy/tests/reader_tests/test_eps_l1b.py index 3d085d5106..4e3435ad23 100644 --- a/satpy/tests/reader_tests/test_eps_l1b.py +++ b/satpy/tests/reader_tests/test_eps_l1b.py @@ -1,6 +1,4 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# Copyright (c) 2019, 2022 Satpy developers +# Copyright (c) 2019, 2022-2023 Satpy developers # # This file is part of satpy. # @@ -15,7 +13,7 @@ # # You should have received a copy of the GNU General Public License along with # satpy. If not, see . -"""Test the eps l1b format.""" +"""Test the EPS l1b format reader.""" import os from contextlib import suppress From ffc5f50b76bea4855f23003373354d9acadc9be6 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Tue, 2 May 2023 09:33:59 +0200 Subject: [PATCH 03/23] Add docu link --- satpy/readers/iasi_l2_eps.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/satpy/readers/iasi_l2_eps.py b/satpy/readers/iasi_l2_eps.py index 466f37a79e..dccaf9d0c4 100644 --- a/satpy/readers/iasi_l2_eps.py +++ b/satpy/readers/iasi_l2_eps.py @@ -28,6 +28,9 @@ class EPSIASIFile(EPSFile): specification, validation reports, and other documents, can be found at the EUMETSAT Data Services at https://data.eumetsat.int/product/EO:EUM:DAT:METOP:IASSND02 + Product format specification for IASI L2 data: https://www.eumetsat.int/media/41105 + + Generic EPS product format specification: https://www.eumetsat.int/media/40048 """ def __init__(self, filename, filename_info, filetype_info): From 7261b8a61c927d0ea265b1835a3849b40b0d3cfb Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Tue, 4 Jul 2023 09:03:00 +0200 Subject: [PATCH 04/23] Added test for parsing XML configuration --- satpy/tests/reader_tests/test_xmlformat.py | 27 ++++++++++++++++++++++ 1 file changed, 27 insertions(+) create mode 100644 satpy/tests/reader_tests/test_xmlformat.py diff --git a/satpy/tests/reader_tests/test_xmlformat.py b/satpy/tests/reader_tests/test_xmlformat.py new file mode 100644 index 0000000000..7b7f5bce91 --- /dev/null +++ b/satpy/tests/reader_tests/test_xmlformat.py @@ -0,0 +1,27 @@ +# Copyright (c) 2018-2020 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 . +"""Tests for reading a format from an XML file.""" + +import pytest + + +@pytest.mark.parametrize("form", ["eps_avhrrl1b_6.5.xml", "eps_iasil2_9.0.xml"]) +def test_parse_format(form): + """Test parsing the XML format.""" + from satpy._config import get_config_path + from satpy.readers.xmlformat import parse_format + filename = get_config_path(form) + parse_format(filename) From 2ab420d39c73e32e15541d24c4de45ee2996e0f5 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Tue, 4 Jul 2023 09:57:02 +0200 Subject: [PATCH 05/23] Add missing definitions for reading IASI L2 For the XML format specification, add missing type configurations to xmlformat.py. Also uncomment parameter definitions from the XML definition, which are needed to interpret the rest of the file. I don't know why those were commented out. --- satpy/etc/eps_iasil2_9.0.xml | 18 +++++++++--------- satpy/readers/xmlformat.py | 8 ++++++-- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/satpy/etc/eps_iasil2_9.0.xml b/satpy/etc/eps_iasil2_9.0.xml index e45415785e..f5219776b5 100644 --- a/satpy/etc/eps_iasil2_9.0.xml +++ b/satpy/etc/eps_iasil2_9.0.xml @@ -48,25 +48,25 @@ - - - - + + + + - + - + - + - - + + diff --git a/satpy/readers/xmlformat.py b/satpy/readers/xmlformat.py index 0c46a3595e..1327a2cae9 100644 --- a/satpy/readers/xmlformat.py +++ b/satpy/readers/xmlformat.py @@ -19,7 +19,7 @@ from __future__ import annotations -from xml.etree.ElementTree import ElementTree +from xml.etree.ElementTree import ElementTree # nosec import numpy as np @@ -28,8 +28,12 @@ TYPEC = {"boolean": ">i1", "integer2": ">i2", "integer4": ">i4", + "vinteger4": ">i4", + "uinteger1": ">u1", "uinteger2": ">u2", - "uinteger4": ">u4", } + "vuinteger2": ">u2", + "uinteger4": ">u4", + "vuinteger4": ">u4"} def process_delimiter(elt, ascii=False): From 741b31de0ea3ecac0b05787070480a356c7229f2 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Wed, 5 Jul 2023 11:53:15 +0200 Subject: [PATCH 06/23] Mroe refactoring and improve tests Change the distribution of responsibility between EPSBaseFileHandler and EPSAVHRRFile. Added a tests for xmlformat.process_array. --- satpy/etc/readers/iasi_l2_eps.yaml | 2 +- satpy/readers/eps_base.py | 77 ++++++++++++- satpy/readers/eps_l1b.py | 120 ++++----------------- satpy/tests/reader_tests/test_xmlformat.py | 31 ++++++ 4 files changed, 130 insertions(+), 100 deletions(-) diff --git a/satpy/etc/readers/iasi_l2_eps.yaml b/satpy/etc/readers/iasi_l2_eps.yaml index 98e419d52a..315770a36f 100644 --- a/satpy/etc/readers/iasi_l2_eps.yaml +++ b/satpy/etc/readers/iasi_l2_eps.yaml @@ -15,7 +15,7 @@ reader: file_types: iasi_l2_eps: - file_reader: !!python/name:satpy.readers.iasi_l2_eps.EPSIASIFile + file_reader: !!python/name:satpy.readers.eps_l2.EPSIASIL2FileHandler # see https://www.eumetsat.int/media/40048 §7 for file pattern file_patterns: - IASI_SND_02_{platform_short_name}_{start_time:%Y%m%d%H%M%SZ}_{end_time:%Y%m%d%H%M%SZ}_{processing_mode}_{disposition_mode}_{creation_time:%Y%m%d%H%M%SZ} diff --git a/satpy/readers/eps_base.py b/satpy/readers/eps_base.py index b00e66a272..d33e13ebeb 100644 --- a/satpy/readers/eps_base.py +++ b/satpy/readers/eps_base.py @@ -21,6 +21,7 @@ import dask.array as da import numpy as np import xarray as xr +from pyresample.geometry import SwathDefinition from satpy._config import get_config_path from satpy.readers.file_handlers import BaseFileHandler @@ -124,9 +125,83 @@ def create_xarray(arr): return res -class EPSFile(BaseFileHandler): +class EPSBaseFileHandler(BaseFileHandler): """Base class for EPS level 1b readers.""" spacecrafts = {"M01": "Metop-B", "M02": "Metop-A", "M03": "Metop-C", } + + xml_conf: str + mdr_subclass: int + + def __init__(self, filename, filename_info, filetype_info): + """Initialize FileHandler.""" + super().__init__(filename, filename_info, filetype_info) + + self.area = None + self._start_time = filename_info['start_time'] + self._end_time = filename_info['end_time'] + self.form = None + self.scanlines = None + self.sections = None + + def _read_all(self): + logger.debug("Reading %s", self.filename) + self.sections, self.form = read_records(self.filename, self.xml_conf) + self.scanlines = self['TOTAL_MDR'] + if self.scanlines != len(self.sections[('mdr', self.mdr_subclass)]): + logger.warning("Number of declared records doesn't match number of scanlines in the file.") + self.scanlines = len(self.sections[('mdr', self.mdr_subclass)]) + + def __getitem__(self, key): + """Get value for given key.""" + for altkey in self.form.scales: + try: + try: + return self.sections[altkey][key] * self.form.scales[altkey][key] + except TypeError: + val = self.sections[altkey][key].item().decode().split("=")[1] + try: + return float(val) * self.form.scales[altkey][key].item() + except ValueError: + return val.strip() + except (KeyError, ValueError): + continue + raise KeyError("No matching value for " + str(key)) + + def keys(self): + """List of reader's keys.""" + keys = [] + for val in self.form.scales.values(): + keys += val.dtype.fields.keys() + return keys + + def get_lonlats(self): + """Get lonlats.""" + if self.area is None: + lons, lats = self.get_full_lonlats() + self.area = SwathDefinition(lons, lats) + self.area.name = '_'.join([self.platform_name, str(self.start_time), + str(self.end_time)]) + return self.area + + @property + def platform_name(self): + """Get platform name.""" + return self.spacecrafts[self["SPACECRAFT_ID"]] + + @property + def sensor_name(self): + """Get sensor name.""" + return self.sensors[self["INSTRUMENT_ID"]] + + @property + def start_time(self): + """Get start time.""" + return self._start_time + + @property + def end_time(self): + """Get end time.""" + return self._end_time diff --git a/satpy/readers/eps_l1b.py b/satpy/readers/eps_l1b.py index d58f29771e..23a6071d7a 100644 --- a/satpy/readers/eps_l1b.py +++ b/satpy/readers/eps_l1b.py @@ -23,10 +23,9 @@ import numpy as np import xarray as xr from dask.delayed import delayed -from pyresample.geometry import SwathDefinition from satpy._compat import cached_property -from satpy.readers.eps_base import EPSFile, XMLFormat, read_records, record_class # noqa +from satpy.readers.eps_base import EPSBaseFileHandler, XMLFormat, read_records, record_class # noqa from satpy.utils import get_legacy_chunk_size logger = logging.getLogger(__name__) @@ -54,7 +53,7 @@ def create_xarray(arr): return res -class EPSAVHRRFile(EPSFile): +class EPSAVHRRFile(EPSBaseFileHandler): """EPS level 1b reader for AVHRR data.""" sensors = {"AVHR": "avhrr-3"} @@ -62,57 +61,25 @@ class EPSAVHRRFile(EPSFile): units = {"reflectance": "%", "brightness_temperature": "K"} + xml_conf = "eps_avhrrl1b_6.5.xml" + mdr_subclass = 2 + def __init__(self, filename, filename_info, filetype_info): """Initialize FileHandler.""" - super(EPSAVHRRFile, self).__init__( - filename, filename_info, filetype_info) - - self.area = None - self._start_time = filename_info['start_time'] - self._end_time = filename_info['end_time'] - self.form = None - self.scanlines = None - self.pixels = None - self.sections = None + super().__init__(filename, filename_info, filetype_info) + self.get_full_angles = functools.lru_cache(maxsize=1)( self._get_full_angles_uncached ) self.get_full_lonlats = functools.lru_cache(maxsize=1)( self._get_full_lonlats_uncached ) + self.pixels = None def _read_all(self): - logger.debug("Reading %s", self.filename) - self.sections, self.form = read_records(self.filename, "eps_avhrrl1b_6.5.xml") - self.scanlines = self['TOTAL_MDR'] - if self.scanlines != len(self.sections[('mdr', 2)]): - logger.warning("Number of declared records doesn't match number of scanlines in the file.") - self.scanlines = len(self.sections[('mdr', 2)]) + super()._read_all() self.pixels = self["EARTH_VIEWS_PER_SCANLINE"] - def __getitem__(self, key): - """Get value for given key.""" - for altkey in self.form.scales: - try: - try: - return self.sections[altkey][key] * self.form.scales[altkey][key] - except TypeError: - val = self.sections[altkey][key].item().decode().split("=")[1] - try: - return float(val) * self.form.scales[altkey][key].item() - except ValueError: - return val.strip() - except (KeyError, ValueError): - continue - raise KeyError("No matching value for " + str(key)) - - def keys(self): - """List of reader's keys.""" - keys = [] - for val in self.form.scales.values(): - keys += val.dtype.fields.keys() - return keys - def _get_full_lonlats_uncached(self): """Get the interpolated longitudes and latitudes.""" raw_lats = np.hstack((self["EARTH_LOCATION_FIRST"][:, [0]], @@ -198,6 +165,19 @@ def get_bounding_box(self): self["EARTH_LOCATION_FIRST"][-1, [1]]]) return lons.ravel(), lats.ravel() + def _get_angle_dataarray(self, key): + """Get an angle dataarray.""" + sun_azi, sun_zen, sat_azi, sat_zen = self.get_full_angles() + if key['name'] == 'solar_zenith_angle': + return create_xarray(sun_zen) + if key['name'] == 'solar_azimuth_angle': + return create_xarray(sun_azi) + if key['name'] == 'satellite_zenith_angle': + return create_xarray(sat_zen) + if key['name'] == 'satellite_azimuth_angle': + return create_xarray(sat_azi) + raise ValueError(f"Unknown angle data-array: {key['name']:s}") + def get_dataset(self, key, info): """Get calibrated channel data.""" if self.sections is None: @@ -227,19 +207,6 @@ def get_dataset(self, key, info): dataset.attrs.update(key.to_dict()) return dataset - def _get_angle_dataarray(self, key): - """Get an angle dataarray.""" - sun_azi, sun_zen, sat_azi, sat_zen = self.get_full_angles() - if key['name'] == 'solar_zenith_angle': - dataset = create_xarray(sun_zen) - elif key['name'] == 'solar_azimuth_angle': - dataset = create_xarray(sun_azi) - if key['name'] == 'satellite_zenith_angle': - dataset = create_xarray(sat_zen) - elif key['name'] == 'satellite_azimuth_angle': - dataset = create_xarray(sat_azi) - return dataset - @cached_property def three_a_mask(self): """Mask for 3A.""" @@ -283,46 +250,3 @@ def _get_calibrated_dataarray(self, key): if mask is not None: dataset = dataset.where(~mask) return dataset - - def get_lonlats(self): - """Get lonlats.""" - if self.area is None: - lons, lats = self.get_full_lonlats() - self.area = SwathDefinition(lons, lats) - self.area.name = '_'.join([self.platform_name, str(self.start_time), - str(self.end_time)]) - return self.area - - @property - def platform_name(self): - """Get platform name.""" - return self.spacecrafts[self["SPACECRAFT_ID"]] - - @property - def sensor_name(self): - """Get sensor name.""" - return self.sensors[self["INSTRUMENT_ID"]] - - @property - def start_time(self): - """Get start time.""" - # return datetime.strptime(self["SENSING_START"], "%Y%m%d%H%M%SZ") - return self._start_time - - @property - def end_time(self): - """Get end time.""" - # return datetime.strptime(self["SENSING_END"], "%Y%m%d%H%M%SZ") - return self._end_time - - -class EPSIASIFile(EPSFile): - """EPS level 2 reader for IASI data. - - Reader for the IASI Level 2 combined sounding products in native format. - - Overview of the data including links to the product user guide, product format - specification, validation reports, and other documents, can be found at the - EUMETSAT Data Services at https://data.eumetsat.int/product/EO:EUM:DAT:METOP:IASSND02 - - """ diff --git a/satpy/tests/reader_tests/test_xmlformat.py b/satpy/tests/reader_tests/test_xmlformat.py index 7b7f5bce91..0340de1a8a 100644 --- a/satpy/tests/reader_tests/test_xmlformat.py +++ b/satpy/tests/reader_tests/test_xmlformat.py @@ -15,6 +15,9 @@ # satpy. If not, see . """Tests for reading a format from an XML file.""" +from xml.etree.ElementTree import Element + +import numpy as np import pytest @@ -25,3 +28,31 @@ def test_parse_format(form): from satpy.readers.xmlformat import parse_format filename = get_config_path(form) parse_format(filename) + + +def test_process_array(monkeypatch): + """Test processing an array tag.""" + from satpy.readers import xmlformat + from satpy.readers.xmlformat import process_array + elt = Element( + "array", + {"name": "SCENE_RADIANCES", + "length": "5", + "labels": "channel1,channel2,channel3,channel4,channel5"}) + elt2 = Element( + "array", + {"length": "$NE", "label": "FOV"}) + elt3 = Element( + "field", + {"type": "integer2", + "scaling-factor": "10^2,10^2,10^4,10^2,10^2"}) + elt.append(elt2) + elt2.append(elt3) + monkeypatch.setattr(xmlformat, "VARIABLES", {"NE": 10}) + (name, tp, shp, scl) = process_array(elt, False) + assert name == "SCENE_RADIANCES" + assert tp == ">i2" + assert shp == (5, 10) + np.testing.assert_allclose( + scl, + np.array([0.01, 0.01, 0.0001, 0.01, 0.01])) From 5e7f32ce1ed4b209ae0bb219f0ac9c8051cb5ca6 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Wed, 5 Jul 2023 12:35:09 +0200 Subject: [PATCH 07/23] Collect dimensions in xmlformat In xmlformat, collect a mapping of variable names against dimension names. --- satpy/readers/xmlformat.py | 22 +++++++++++----------- satpy/tests/reader_tests/test_xmlformat.py | 6 ++++-- 2 files changed, 15 insertions(+), 13 deletions(-) diff --git a/satpy/readers/xmlformat.py b/satpy/readers/xmlformat.py index 1327a2cae9..15277bf5ba 100644 --- a/satpy/readers/xmlformat.py +++ b/satpy/readers/xmlformat.py @@ -36,12 +36,12 @@ "vuinteger4": ">u4"} -def process_delimiter(elt, ascii=False): +def process_delimiter(elt, ascii=False, dims={}): # noqa: B006 """Process a 'delimiter' tag.""" del elt, ascii -def process_field(elt, ascii=False): +def process_field(elt, ascii=False, dims={}): # noqa: B006 """Process a 'field' tag.""" # NOTE: if there is a variable defined in this field and it is different # from the default, we could change the value and restart. @@ -69,28 +69,28 @@ def process_field(elt, ascii=False): return ((elt.get("name"), current_type, scale)) -def process_array(elt, ascii=False): +def process_array(elt, ascii=False, dims={}): # noqa: B006 """Process an 'array' tag.""" - del ascii chld = list(elt) if len(chld) > 1: raise ValueError() chld = chld[0] try: - name, current_type, scale = CASES[chld.tag](chld) + name, current_type, scale = CASES[chld.tag](chld, ascii, dims) size = None except ValueError: - name, current_type, size, scale = CASES[chld.tag](chld) + name, current_type, size, scale = CASES[chld.tag](chld, ascii, dims) del name myname = elt.get("name") or elt.get("label") if elt.get("length").startswith("$"): - length = int(VARIABLES[elt.get("length")[1:]]) + dimname = elt.get("length")[1:] + length = int(VARIABLES[dimname]) + dims[myname] = dimname else: length = int(elt.get("length")) if size is not None: return (myname, current_type, (length, ) + size, scale) - else: - return (myname, current_type, (length, ), scale) + return (myname, current_type, (length, ), scale) CASES = {"delimiter": process_delimiter, @@ -143,7 +143,7 @@ def to_scales(val): return scales -def parse_format(xml_file): +def parse_format(xml_file, dims={}): # noqa: B006 """Parse the xml file to create types, scaling factor types, and scales.""" tree = ElementTree() tree.parse(xml_file) @@ -156,7 +156,7 @@ def parse_format(xml_file): ascii = (prod.tag in ["mphr", "sphr"]) res = [] for i in prod: - lres = CASES[i.tag](i, ascii) + lres = CASES[i.tag](i, ascii, dims) if lres is not None: res.append(lres) types_scales[(prod.tag, int(prod.get("subclass")))] = res diff --git a/satpy/tests/reader_tests/test_xmlformat.py b/satpy/tests/reader_tests/test_xmlformat.py index 0340de1a8a..39a649dbf0 100644 --- a/satpy/tests/reader_tests/test_xmlformat.py +++ b/satpy/tests/reader_tests/test_xmlformat.py @@ -27,7 +27,7 @@ def test_parse_format(form): from satpy._config import get_config_path from satpy.readers.xmlformat import parse_format filename = get_config_path(form) - parse_format(filename) + (_, _, _) = parse_format(filename) def test_process_array(monkeypatch): @@ -49,10 +49,12 @@ def test_process_array(monkeypatch): elt.append(elt2) elt2.append(elt3) monkeypatch.setattr(xmlformat, "VARIABLES", {"NE": 10}) - (name, tp, shp, scl) = process_array(elt, False) + dims = {} + (name, tp, shp, scl) = process_array(elt, False, dims) assert name == "SCENE_RADIANCES" assert tp == ">i2" assert shp == (5, 10) np.testing.assert_allclose( scl, np.array([0.01, 0.01, 0.0001, 0.01, 0.01])) + assert dims == {"FOV": "NE"} From 7cf3c3d7e406171bd83ae68f26e3b93df488d937 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Wed, 5 Jul 2023 14:50:26 +0200 Subject: [PATCH 08/23] Make xmlformat.parse_format a method Turn xmlformat.parse_format from a module level function to a method on XMLFormat. Remove the module level dictionary VARIABLES and make this an instance attribute on the XMLFormat class. --- satpy/readers/xmlformat.py | 84 ++++++++++------------ satpy/tests/reader_tests/test_xmlformat.py | 32 +++++++-- 2 files changed, 63 insertions(+), 53 deletions(-) diff --git a/satpy/readers/xmlformat.py b/satpy/readers/xmlformat.py index 15277bf5ba..4b14bbd7a7 100644 --- a/satpy/readers/xmlformat.py +++ b/satpy/readers/xmlformat.py @@ -23,8 +23,6 @@ import numpy as np -VARIABLES: dict[str, str] = {} - TYPEC = {"boolean": ">i1", "integer2": ">i2", "integer4": ">i4", @@ -36,12 +34,11 @@ "vuinteger4": ">u4"} -def process_delimiter(elt, ascii=False, dims={}): # noqa: B006 +def process_delimiter(elt, variables, ascii=False): """Process a 'delimiter' tag.""" - del elt, ascii -def process_field(elt, ascii=False, dims={}): # noqa: B006 +def process_field(elt, variables, ascii=False): """Process a 'field' tag.""" # NOTE: if there is a variable defined in this field and it is different # from the default, we could change the value and restart. @@ -69,23 +66,21 @@ def process_field(elt, ascii=False, dims={}): # noqa: B006 return ((elt.get("name"), current_type, scale)) -def process_array(elt, ascii=False, dims={}): # noqa: B006 +def process_array(elt, variables, ascii=False): """Process an 'array' tag.""" chld = list(elt) if len(chld) > 1: raise ValueError() chld = chld[0] try: - name, current_type, scale = CASES[chld.tag](chld, ascii, dims) + name, current_type, scale = CASES[chld.tag](chld, variables, ascii) size = None except ValueError: - name, current_type, size, scale = CASES[chld.tag](chld, ascii, dims) - del name + name, current_type, size, scale = CASES[chld.tag](chld, variables, ascii) myname = elt.get("name") or elt.get("label") if elt.get("length").startswith("$"): dimname = elt.get("length")[1:] - length = int(VARIABLES[dimname]) - dims[myname] = dimname + length = int(variables[dimname]) else: length = int(elt.get("length")) if size is not None: @@ -143,36 +138,6 @@ def to_scales(val): return scales -def parse_format(xml_file, dims={}): # noqa: B006 - """Parse the xml file to create types, scaling factor types, and scales.""" - tree = ElementTree() - tree.parse(xml_file) - for param in tree.find("parameters"): - VARIABLES[param.get("name")] = param.get("value") - - types_scales = {} - - for prod in tree.find("product"): - ascii = (prod.tag in ["mphr", "sphr"]) - res = [] - for i in prod: - lres = CASES[i.tag](i, ascii, dims) - if lres is not None: - res.append(lres) - types_scales[(prod.tag, int(prod.get("subclass")))] = res - - types = {} - stypes = {} - scales = {} - - for key, val in types_scales.items(): - types[key] = to_dtype(val) - stypes[key] = to_scaled_dtype(val) - scales[key] = to_scales(val) - - return types, stypes, scales - - def _apply_scales(array, scales, dtype): """Apply scales to the array.""" new_array = np.empty(array.shape, dtype) @@ -187,12 +152,14 @@ def _apply_scales(array, scales, dtype): return new_array -class XMLFormat(object): +class XMLFormat: """XMLFormat object.""" def __init__(self, filename): """Init the format reader.""" - self.types, self.stypes, self.scales = parse_format(filename) + self.dims = {} + self.variables = {} + self.types, self.stypes, self.scales = self.parse_format(filename) self.translator = {} @@ -207,6 +174,31 @@ def apply_scales(self, array): """Apply scales to *array*.""" return _apply_scales(array, *self.translator[array.dtype]) - -if __name__ == '__main__': - pass + def parse_format(self, xml_file): + """Parse the xml file to create types, scaling factor types, and scales.""" + tree = ElementTree() + tree.parse(xml_file) + for param in tree.find("parameters"): + self.variables[param.get("name")] = param.get("value") + + types_scales = {} + + for prod in tree.find("product"): + ascii = (prod.tag in ["mphr", "sphr"]) + res = [] + for i in prod: + lres = CASES[i.tag](i, self.variables, ascii) + if lres is not None: + res.append(lres) + types_scales[(prod.tag, int(prod.get("subclass")))] = res + + types = {} + stypes = {} + scales = {} + + for key, val in types_scales.items(): + types[key] = to_dtype(val) + stypes[key] = to_scaled_dtype(val) + scales[key] = to_scales(val) + + return types, stypes, scales diff --git a/satpy/tests/reader_tests/test_xmlformat.py b/satpy/tests/reader_tests/test_xmlformat.py index 39a649dbf0..1b22089ea5 100644 --- a/satpy/tests/reader_tests/test_xmlformat.py +++ b/satpy/tests/reader_tests/test_xmlformat.py @@ -25,14 +25,13 @@ def test_parse_format(form): """Test parsing the XML format.""" from satpy._config import get_config_path - from satpy.readers.xmlformat import parse_format + from satpy.readers.xmlformat import XMLFormat filename = get_config_path(form) - (_, _, _) = parse_format(filename) + XMLFormat(filename) def test_process_array(monkeypatch): """Test processing an array tag.""" - from satpy.readers import xmlformat from satpy.readers.xmlformat import process_array elt = Element( "array", @@ -48,13 +47,32 @@ def test_process_array(monkeypatch): "scaling-factor": "10^2,10^2,10^4,10^2,10^2"}) elt.append(elt2) elt2.append(elt3) - monkeypatch.setattr(xmlformat, "VARIABLES", {"NE": 10}) - dims = {} - (name, tp, shp, scl) = process_array(elt, False, dims) + variables = {"NE": 10} +# dims = {} + (name, tp, shp, scl) = process_array(elt, variables, False) assert name == "SCENE_RADIANCES" assert tp == ">i2" assert shp == (5, 10) np.testing.assert_allclose( scl, np.array([0.01, 0.01, 0.0001, 0.01, 0.01])) - assert dims == {"FOV": "NE"} +# assert dims == {"FOV": "NE"} + + elt = Element( + "array", + {"name": "ATMOSPHERIC_WATER_VAPOUR", + "length": "120"}) + elt2 = Element( + "array", + {"length": "$NLQ"}) + elt3 = Element( + "field", + {"type": "uinteger4", + "scaling-factor": "10^7", + "units": "kg/kg"}) + elt.append(elt2) + elt2.append(elt3) + variables = {"NLQ": 125} +# dims = {} + (name, tp, shp, scl) = process_array(elt, variables, False) +# assert dims == {"ATMOSPHERIC_WATER_VAPOUR": "NLQ"} From 540fdcec73e51dbc8d83b4f0e179edb845f1ecef Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Thu, 6 Jul 2023 11:07:45 +0200 Subject: [PATCH 09/23] IASI L2 EPS reader from EUMETSAT Take the IASI L2 EPS reader from EUMETSAT and integrate into Satpy. The code was obtained from https://gitlab.eumetsat.int/open-source/data-tailor-plugins/epct_plugin_gis where it is licensed under the Apache License Version 2.0, January 2004. It was developed by B-Open Solutions srl for EUMETSAT under contract EUM/C0/17/4600001943/0PN and is Copyright (c) 2017-2022 EUMETSAT. To the best of my knowledge, the Apache License allows for its integration into GPL code. --- satpy/etc/readers/iasi_l2_eps.yaml | 2 +- satpy/readers/epsnative_reader.py | 285 +++++++++++ satpy/readers/iasi_l2_eps.py | 749 ++++++++++++++++++++++++++++- 3 files changed, 1030 insertions(+), 6 deletions(-) create mode 100644 satpy/readers/epsnative_reader.py diff --git a/satpy/etc/readers/iasi_l2_eps.yaml b/satpy/etc/readers/iasi_l2_eps.yaml index 315770a36f..160e060c70 100644 --- a/satpy/etc/readers/iasi_l2_eps.yaml +++ b/satpy/etc/readers/iasi_l2_eps.yaml @@ -15,7 +15,7 @@ reader: file_types: iasi_l2_eps: - file_reader: !!python/name:satpy.readers.eps_l2.EPSIASIL2FileHandler + file_reader: !!python/name:satpy.readers.iasi_l2_eps.EPSIASIL2FileHandler # see https://www.eumetsat.int/media/40048 §7 for file pattern file_patterns: - IASI_SND_02_{platform_short_name}_{start_time:%Y%m%d%H%M%SZ}_{end_time:%Y%m%d%H%M%SZ}_{processing_mode}_{disposition_mode}_{creation_time:%Y%m%d%H%M%SZ} diff --git a/satpy/readers/epsnative_reader.py b/satpy/readers/epsnative_reader.py new file mode 100644 index 0000000000..cb171e42fb --- /dev/null +++ b/satpy/readers/epsnative_reader.py @@ -0,0 +1,285 @@ +# Copyright 2017-2022, European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) +# 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 . + +# This module is based on source code obtained from the +# epct_plugin_gis package developed by B-Open Solutions srl for EUMETSAT under +# contract EUM/C0/17/4600001943/0PN and released under Apache License +# Version 2.0, January 2004, http://www.apache.org/licenses/. The original +# source including revision history and details on authorship can be found at +# https://gitlab.eumetsat.int/open-source/data-tailor-plugins/epct_plugin_gis + +"""FIXME DOC. + +The ``epsnative_reader`` module provides the support to read from EPS native products records info +and metadata. +""" + +import collections +import datetime +import os + +import numpy as np +import pandas as pd +import yaml + +from satpy._config import get_config_path + +dtype_to_gdal_type = { + ">i1": "Byte", + ">i2": "Int16", + ">i4": "Int32", + ">u1": "Byte", + ">u2": "UInt16", + ">u4": "UInt32", + ">f4": "Float32", +} +eps_classes = ["reserved", "mphr", "sphr", "ipr", "geadr", "giadr", "veadr", "viadr", "mdr"] +eps_type_to_dtype = { + "boolean": ">i1", + "ubyte": ">u1", + "integer": ">i", + "integer1": ">i1", + "integer2": ">i2", + "integer4": ">i4", + "uinteger": ">u1", + "uinteger1": ">u1", + "uinteger2": ">u2", + "uinteger4": ">u4", + "vuinteger2": "byte,>u2", + "vinteger2": "byte,>i2", + "vinteger4": "byte,>i4", + "short_cds_time": ">u2,>u4", + "enumerated": ">i1", + "bitst": ">ux", +} +grh_type = [ + ("RECORD_CLASS", "|i1"), + ("INSTRUMENT_GROUP", "|i1"), + ("RECORD_SUBCLASS", "|i1"), + ("RECORD_SUBCLASS_VERSION", "|i1"), + ("RECORD_SIZE", ">u4"), + ("RECORD_START_TIME", eps_type_to_dtype["short_cds_time"]), + ("RECORD_STOP_TIME", eps_type_to_dtype["short_cds_time"]), +] +ipr_type = [ + ("TARGET_RECORD_CLASS", "|i1"), + ("TARGET_INSTRUMENT_GROUP", "|i1"), + ("TARGET_RECORD_SUBCLASS", "|i1"), + ("TARGET_RECORD_OFFSET", ">u4"), +] +NODATA_VALUE = -9999 +UINTEGER1_BIT = 8 +FORMAT_TO_MDR_CLASS = { + "ASCATL1SZF": {13: ("mdr", 3, 5), 12: ("mdr", 3, 4), 11: ("mdr", 3, 4), 10: ("mdr", 3, 4)}, + "ASCATL1SZO": {13: ("mdr", 2, 4), 12: ("mdr", 2, 3), 11: ("mdr", 2, 3), 10: ("mdr", 2, 3)}, + "ASCATL1SZR": {13: ("mdr", 1, 4), 12: ("mdr", 1, 3), 11: ("mdr", 1, 3), 10: ("mdr", 1, 3)}, + "ASCATL2SMO": {12: ("mdr", 5, 1), 11: ("mdr", 5, 1), 10: ("mdr", 5, 0)}, + "ASCATL2SMR": {12: ("mdr", 4, 1), 11: ("mdr", 4, 1), 10: ("mdr", 4, 0)}, +} + + +def get_class_tuple(class_string): + """FIXME DOC. + + :param class_string: + :return: + """ + class_name = class_string.split("_") + try: + class_name[1:] = map(int, class_name[1:]) + except ValueError: + class_name = (class_name[0], "_".join(class_name[1:])) + return tuple(class_name) + + +def assemble_descriptor(product): + """FIXME DOC. + + :param product: + :return: + """ + csv_list = [ + get_config_path("readers/eps_native_format/IASISND02/mdr_1_4.csv"), + get_config_path("readers/eps_native_format/IASISND02/giadr_1_4.csv"), + get_config_path("readers/eps_native_format/mphr_0_2.csv")] + descriptor = {} + for csv in csv_list: + class_name = get_class_tuple(os.path.basename(csv).split(".")[0]) + ds = pd.read_csv(csv, sep=",") + ds.columns = ds.columns.str.replace(" ", "_") + descriptor[tuple(class_name)] = ds + return descriptor + + +def scds_to_datetime(days, milliseconds): + """FIXME DOC. + + :param int days: + :param int milliseconds: + :return datetime.datetime: + """ + epoch = datetime.datetime(2000, 1, 1) + return epoch + datetime.timedelta(days=int(days), milliseconds=int(milliseconds)) + + +def grh_reader(eps_fileobj): + """FIXME DOC. + + :param eps_fileobj: + :return: + """ + grh_array = np.fromfile(eps_fileobj, np.dtype(grh_type), 1) + # When the EOF is reached + if grh_array.size == 0: + return () + rec_class = eps_classes[grh_array["RECORD_CLASS"][0]] + rec_subclass = grh_array["RECORD_SUBCLASS"][0] + rec_subclass_version = grh_array["RECORD_SUBCLASS_VERSION"][0] + rec_size = grh_array["RECORD_SIZE"][0] + rec_start_time = scds_to_datetime(*grh_array["RECORD_START_TIME"][0]) + rec_stop_time = scds_to_datetime(*grh_array["RECORD_STOP_TIME"][0]) + return rec_class, rec_subclass, rec_subclass_version, rec_size, rec_start_time, rec_stop_time + + +def mphr_reader(input_product): + """FIXME DOC. + + :param input_product: + :return: + """ + mphr_descriptor = pd.read_csv( + get_config_path("readers/eps_native_format/mphr_0_2.csv"), + sep=",") + mphr_content = collections.OrderedDict() + with open(input_product, "rb") as eps_fileobj: + eps_fileobj.seek(20) + for _, row in mphr_descriptor.iterrows(): + # empty row or general header + if np.isnan(row["TYPE_SIZE"]) or row["FIELD"] == "RECORD_HEADER": + continue + else: + dtype = f"S30,S2,S{int(row['TYPE_SIZE'])},S1" + row_content = np.fromfile(eps_fileobj, dtype, 1) + key = row_content[0][0].decode().strip() + value = row_content[0][2].decode().strip(" x") + if value.lstrip("-").isdigit(): + value = int(value) + if np.isfinite(row["SF"]): + value /= 10 ** row["SF"] + mphr_content[key] = value + return mphr_content + + +def first_class_occurrence(input_product, class_name): + """FIXME DOC. + + :param input_product: + :param class_name: + :return: + + """ + header_size = 20 # bytes + offset = None + with open(input_product, "rb") as eps_fileobj: + while True: + grh = grh_reader(eps_fileobj) + if len(grh) == 0: + break + elif grh[0] != class_name: + class_size = grh[3] + eps_fileobj.seek(class_size - header_size, 1) + continue + else: + offset = eps_fileobj.tell() - header_size + break + return grh, offset + + +def ipr_reader(eps_fileobj): + """FIXME DOC. + + :param eps_fileobj: + :return: + """ + ipr_array = np.fromfile(eps_fileobj, np.dtype(ipr_type), 1) + data = list(ipr_array[0]) + data[0] = eps_classes[data[0]] + return data + + +def read_ipr_sequence(input_product): + """FIXME DOC. + + :param input_product: + :return: + """ + dummy_id = 13 + _, ipr_offset = first_class_occurrence(input_product, "ipr") + ipr_sequence = [] + with open(input_product, "rb") as eps_fileobj: + eps_fileobj.seek(ipr_offset) + while True: + grh = grh_reader(eps_fileobj) + if grh[0] != "ipr": + break + ipr_content = ipr_reader(eps_fileobj) + info = { + "class": (ipr_content[0], ipr_content[2]), + "instrument": ipr_content[1], + "offset": ipr_content[-1], + "is_dummy": ipr_content[1] == dummy_id, + } + ipr_sequence.append(info) + return ipr_sequence + + +def reckon_dtype(eps_type): + """FIXME DOC. + + Convert EPS type into a suitable string which represents a ``numpy dtype`` object, according + to whether the tag belongs to an ASCII record or not. For ASCII records see: "EPS Generic + Product Format Specification", par. 4.3.1. + + :param str eps_type: + :return string: string which represents a ``numpy dtype`` object + """ + # remove "-" character, e.g. "u-integer2" + eps_type = eps_type.replace("-", "") + if eps_type.startswith("bitst"): + eps_type, length = eps_type.replace("(", " ").split() + length = int(length.split(")")[0]) + dtype = eps_type_to_dtype.get(eps_type) + # if the length in bytes is an odd number the numpy dtype is a sequence of uinteger1 + if ((length // UINTEGER1_BIT) % 2 != 0) or (length // UINTEGER1_BIT) > UINTEGER1_BIT: + dtype = ",".join([">u1"] * (length // UINTEGER1_BIT)) + else: + dtype = dtype.replace("x", str(length // UINTEGER1_BIT)) + else: + dtype = eps_type_to_dtype.get(eps_type) + return dtype + + +def bands_to_records_reader(product): + """FIXME DOC. + + :param product: + :return: + """ + band_to_record_path = get_config_path("readers/eps_native_format/IASISND02/band_to_record.yaml") + with open(band_to_record_path, "rb") as fp: + contents = fp.read() + return yaml.safe_load(contents) diff --git a/satpy/readers/iasi_l2_eps.py b/satpy/readers/iasi_l2_eps.py index dccaf9d0c4..2c7bf57d23 100644 --- a/satpy/readers/iasi_l2_eps.py +++ b/satpy/readers/iasi_l2_eps.py @@ -1,4 +1,6 @@ # Copyright (c) 2023 Satpy developers +# Copyright (c) 2017-2022, European Organisation for the Exploitation of +# Meteorological Satellites (EUMETSAT) # # This file is part of satpy. # @@ -14,12 +16,27 @@ # You should have received a copy of the GNU General Public License along with # satpy. If not, see . +# Parts of this module are based on source code obtained from the +# epct_plugin_gis package developed by B-Open Solutions srl for EUMETSAT under +# contract EUM/C0/17/4600001943/0PN and released under Apache License +# Version 2.0, January 2004, http://www.apache.org/licenses/. The original +# source including revision history and details on authorship can be found at +# https://gitlab.eumetsat.int/open-source/data-tailor-plugins/epct_plugin_gis + """Reader for EPS level 2 IASI data. Uses xml files as a format description.""" -from .eps_base import EPSFile +import collections +import datetime +import itertools + +import numpy as np +import xarray as xr +from . import epsnative_reader +from .file_handlers import BaseFileHandler -class EPSIASIFile(EPSFile): + +class EPSIASIL2FileHandler(BaseFileHandler): """EPS level 2 reader for IASI data. Reader for the IASI Level 2 combined sounding products in native format. @@ -33,6 +50,728 @@ class EPSIASIFile(EPSFile): Generic EPS product format specification: https://www.eumetsat.int/media/40048 """ - def __init__(self, filename, filename_info, filetype_info): - """Initialise Filehandler.""" - raise NotImplementedError() + # TODO: + # - make dask-friendly / lazy + # - only read variables that are requested + + _nc = None + + def get_dataset(self, dataid, dataset_info): + """Get dataset.""" + if self._nc is None: + self._nc = self._get_netcdf_dataset() + da = self._nc[dataid["name"]] + # scale factor is not relevant and gets in the way + da.attrs.pop("scale_factor", None) + return da + + def _get_netcdf_dataset(self): + """Get full NetCDF dataset.""" + input_product = self.filename + descriptor = epsnative_reader.assemble_descriptor("IASISND02") + ipr_sequence = epsnative_reader.read_ipr_sequence(input_product) + first_mdr_class = [cl for cl in ipr_sequence if cl["class"] == ("mdr", 1)][0] + first_mdr_offset = first_mdr_class["offset"] + with open(input_product, "rb") as epsfile_obj: + data_before_errors, algo_data, errors_data = read_product_data( + epsfile_obj, descriptor, first_mdr_offset, self.start_time, + self.end_time) + giadr = read_giadr(input_product, descriptor) + ds = create_netcdf_dataset(data_before_errors, algo_data, errors_data, giadr) + return ds + + def available_datasets(self, configured_datasets): + """Get available datasets.""" + # FIXME: do this without converting/reading the file — maybe hardcode + # still? + common = {"file_type": "iasi_l2_eps", "resolution": 12000} + if self._nc is None: + self._nc = self._get_netcdf_dataset() + for var in self._nc.data_vars: + yield (True, {"name": var} | common | self._nc[var].attrs) + + +missing_values = { + "CO_CP_AIR": -2, + "CO_CP_CO_A": -2, + "CO_X_CO": 9.96921e36, + "CO_H_EIGENVALUES": 9.96921e36, + "CO_H_EIGENVECTORS": 9.96921e36, +} +var_2nd_dim = { + "co_cp_air": "co_nbr", + "co_cp_co_a": "co_nbr", + "co_x_co": "co_nbr", + "co_h_eigenvalues": "co_nbr", + "co_h_eigenvectors": "co_nbr", + "ozone_error": "nerr", + "temperature_error": "nerr", + "water_vapour_error": "nerr", +} +var_3rd_dim = { + "fg_atmospheric_temperature": "nlt", + "fg_atmospheric_water_vapour": "nlq", + "fg_atmospheric_ozone": "nlo", + "atmospheric_temperature": "nlt", + "atmospheric_water_vapour": "nlq", + "atmospheric_ozone": "nlo", +} +# supported major version format +MAJOR_VERSION = 11 +NR_FOV = 120 + + +def read_values(eps_fileobj, row, reshape=False): + """FIXME DOC. + + :param _io.BinaryIO eps_fileobj: + :param pandas.core.series.Series row: + :param bool reshape: + :return numpy.ndarray: + """ + dtype = epsnative_reader.reckon_dtype(row["TYPE"].strip()) + shape = [int(row[k]) for k in row.keys() if k.lower().startswith("dim")] + shape = [k for k in shape if k > 1] + count = int(np.product(shape)) + values = np.fromfile(eps_fileobj, dtype, count) + # if np.isfinite(row['SF']): + # values = values / 10 ** row['SF'] + if row["TYPE"].startswith("v-"): + values = values["f1"] * 10.0 ** -values["f0"] + if row["TYPE"].strip().lower() == "short_cds_time": + values = [epsnative_reader.scds_to_datetime(*v) for v in values] + values = values[0] if len(values) == 1 else values + if reshape: + try: + values = values.reshape(tuple(reversed(shape))).T + except AttributeError: + values = values + return values + + +def read_giadr(input_product, descriptor, **kwargs): + """FIXME DOC. + + :param input_product: + :param descriptor: + :param kwargs: + :return: + """ + grh_size = 20 # bytes + ipr_sequence = epsnative_reader.read_ipr_sequence(input_product) + giadr_class = [cl for cl in ipr_sequence if "giadr" in cl["class"]][0] + with open(input_product, "rb") as eps_fileobj: + eps_fileobj.seek(giadr_class["offset"] + grh_size) + class_data = collections.OrderedDict() + giadr_descriptor = descriptor.get(("giadr", 1, 4)) + for _, row in giadr_descriptor.iterrows(): + # empty row or general header + if np.isnan(row["FIELD_SIZE"]) or row["FIELD"] == "RECORD_HEADER": + continue + else: + class_data[row["FIELD"].strip()] = { + "description": row["DESCRIPTION"], + "units": row["UNITS"], + "values": read_values(eps_fileobj, row, **kwargs), + } + return class_data + + +def set_values_in_mdr_descriptor(epsfile_obj, mdr_descriptor, mdr_class_offset, field_name): + """FIXME DOC. + + :param _io.BinaryIO epsfile_obj: + :param pandas.core.frame.DataFrame mdr_descriptor: + :param int mdr_class_offset: + :param str field_name: + :return list: + """ + row = mdr_descriptor.loc[mdr_descriptor["FIELD"] == field_name] + row_idx = row.index[0] + row = row.squeeze() + # read values + epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) + value = read_values(epsfile_obj, row).astype(int)[0] + # set the read values in MDR-descriptor + mdr_descriptor.loc[mdr_descriptor["DIM2"] == field_name.lower(), "DIM2"] = value + return value, row_idx + + +def update_mdr_descriptor(mdr_descriptor, start_row_idx): + """FIXME DOC. + + :param pandas.core.frame.DataFrame mdr_descriptor: + :param int start_row_idx: + :return pandas.core.frame.DataFrame: + """ + for idx in range(start_row_idx, len(mdr_descriptor)): + row = mdr_descriptor.iloc[idx] + if np.isnan(row["FIELD_SIZE"]) or np.isnan(row["OFFSET"]): + previous_row = mdr_descriptor.iloc[idx - 1] + try: + if np.isnan(row["FIELD_SIZE"]): + dims = np.int_(row.loc[["DIM1", "DIM2", "DIM3"]].values) + field_size = np.prod(dims) * int(row["TYPE_SIZE"]) + # set the FILED_SIZE value + mdr_descriptor.at[idx, "FIELD_SIZE"] = field_size + if np.isnan(row["OFFSET"]): + # set the OFFSET value + offset = previous_row["OFFSET"] + previous_row["FIELD_SIZE"] + mdr_descriptor.at[idx, "OFFSET"] = int(offset) + except ValueError: + continue + return mdr_descriptor + + +def update_mdr_descriptor_for_all_parameters(epsfile_obj, mdr_descriptor, mdr_class_offset): + """FIXME DOC. + + :param _io.BinaryIo epsfile_obj: + :param pandas.core.frame.DataFrame mdr_descriptor: + :param int mdr_class_offset: + :return (dict, pandas.core.frame.DataFrame): + """ + params = {} + for field_name in ["NERR", "CO_NBR", "HNO3_NBR", "O3_NBR"]: + value, row_idx = set_values_in_mdr_descriptor( + epsfile_obj, mdr_descriptor, mdr_class_offset, field_name + ) + params[field_name] = value + update_mdr_descriptor(mdr_descriptor, row_idx) + return params, mdr_descriptor + + +def read_errors(epsfile_obj, mdr_descriptor, mdr_class_offset, max_nerr): + """FIXME DOC. + + :param _io.BinaryIo epsfile_obj: + :param pandas.core.frame.DataFrame mdr_descriptor: + :param int mdr_class_offset: + :param int max_nerr: + :return dict: + """ + missing_value = -2147483646 + fields_to_read = [ + "TEMPERATURE_ERROR", + "WATER_VAPOUR_ERROR", + "OZONE_ERROR", + "SURFACE_Z", + ] + errors = {} + for field in fields_to_read: + row = mdr_descriptor.loc[mdr_descriptor["FIELD"] == field].squeeze() + epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) + values = read_values(epsfile_obj, row, reshape=True) + if field != "SURFACE_Z": + if row["DIM2"] == 0: + values = np.full((row["DIM1"], max_nerr), missing_value) + else: + values = values.reshape(-1, 1) if values.ndim == 1 else values + values = np.pad( + values, + ((0, 0), (0, max_nerr - values.shape[1])), + mode="constant", + constant_values=missing_value, + ) + errors[field] = values + return errors + + +def read_algorithm_sections( + epsfile_obj, mdr_descriptor, mdr_class_offset, vars_must_be_extend_to_fov, max_nerr +): + """FIXME DOC. + + Read the following algorithm sections from an MDR record (i.e. row): FORLI-CO, FORLI-HNO3, + FORLI-O3, BRESCIA-SO2 + + :param _io.BinaryIO epsfile_obj: + :param pandas.core.frame.DataFrame mdr_descriptor: + :param int mdr_class_offset: + :param list vars_must_be_extend_to_fov: + :param int max_nerr: + :return (dict, dict): + """ + params, mdr_descriptor = update_mdr_descriptor_for_all_parameters( + epsfile_obj, mdr_descriptor.copy(), mdr_class_offset + ) + errors = read_errors(epsfile_obj, mdr_descriptor, mdr_class_offset, max_nerr) + res = {"params": params} + # data about HNO3 and O3 are skipped because the NetCDF file does not include those data + for prefix in ["CO_"]: # ["CO", "HNO3_", "O3_"] + param_name = [k for k in params.keys() if k.startswith(prefix)][0] + if params[param_name] != 0: + section_rows = mdr_descriptor.loc[mdr_descriptor["FIELD"].str.startswith(prefix)] + section = {} + for _, row in section_rows.iterrows(): + epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) + values = read_values(epsfile_obj, row, reshape=True) + if row["FIELD"] in vars_must_be_extend_to_fov: + values = values.reshape(-1, 1) if values.ndim == 1 else values + values = np.pad( + values, + ((0, 0), (0, NR_FOV - values.shape[1])), + mode="constant", + constant_values=missing_values[row["FIELD"]], + ) + section[row["FIELD"]] = values + res[prefix[:-1]] = section + return res, errors + + +def reckon_valid_algorithm_sections(algorithms_data): + """FIXME DOC. + + Return the value for each MDR-descriptor parameter for each row. An algorithm section is + valid only if at least 1 correlated parameter (e.g 'CO_NBR' for FORLI-CO algorithm) is not + equal to zero + + :param list algorithms_data: + :return dict: + """ + # data about HNO3 and O3 are skipped because the NetCDF file does not include those data + params_data = { + "CO_NBR": {"values": [], "is_valid": None}, + # "HNO3_NBR": {"values": [], "is_valid": None}, + # "O3_NBR": {"values": [], "is_valid": None}, + } + for row in algorithms_data: + for param_name in params_data: + params_data[param_name]["values"].append(row["params"][param_name]) + for param_name in params_data: + total = sum(params_data[param_name]["values"]) + params_data[param_name]["is_valid"] = total > 0 + return params_data + + +def fill_invalid_rows(algorithms_data): + """FIXME DOC. + + :param list algorithms_data: + :return list: + """ + params_data = reckon_valid_algorithm_sections(algorithms_data) + for param, data in params_data.items(): + # only valid algorithms will be saved in the NetCDF file + if data["is_valid"]: + first_valid_row_idx = np.where(np.array(data["values"]) != 0)[0][0] + section_name = param.split("_")[0] + dummy_section = algorithms_data[first_valid_row_idx][section_name].copy() + for key, values in dummy_section.items(): + dummy_section[key] = np.full(values.shape, np.nan) + invalid_rows = np.where(np.array(data["values"]) == 0)[0] + for row_idx in invalid_rows: + algorithms_data[row_idx][section_name] = dummy_section + return algorithms_data + + +def initialize_stacked_algorithms_output(algorithms_data): + """FIXME DOC. + + :param list algorithms_data: + :return dict: + """ + stacked_data = {} + for algorithm in algorithms_data[0]: + if algorithm == "params": + continue + stacked_data[algorithm] = algorithms_data[0][algorithm].copy() + for var in stacked_data[algorithm]: + stacked_data[algorithm][var] = [] + return stacked_data + + +def stack_algorithm_sections_along_rows(algorithms_data, stacked_data): + """FIXME DOC. + + :param list algorithms_data: + :param dict stacked_data: + :return: + """ + # collect values from each row + for idx, row in enumerate(algorithms_data): + for algorithm in row.keys(): + if algorithm == "params": + continue + for var, values in row[algorithm].items(): + stacked_data[algorithm][var].append(values) + # free memory deleting read values + algorithms_data[idx] = None + # stack data along rows + for algotirthm in stacked_data: + for var, values in stacked_data[algotirthm].items(): + stacked_data[algotirthm][var] = np.stack(values, axis=0) + return stacked_data + + +def get_vars_must_be_extend_to_fov(mdr_descriptor): + """FIXME DOC. + + :param pandas.core.frame.DataFrame mdr_descriptor: + :return list: + """ + co_nbr_rows = list(mdr_descriptor.loc[mdr_descriptor["DIM2"] == "co_nbr", "FIELD"]) + hno3_nbr_rows = list(mdr_descriptor.loc[mdr_descriptor["DIM2"] == "hno3_nbr", "FIELD"]) + o3_nbr_rows = list(mdr_descriptor.loc[mdr_descriptor["DIM2"] == "o3_nbr", "FIELD"]) + vars_must_be_extend_to_fov = list( + itertools.chain.from_iterable([co_nbr_rows, hno3_nbr_rows, o3_nbr_rows]) + ) + return vars_must_be_extend_to_fov + + +def stack_non_algorithm_data(data): + """FIXME DOC. + + :param list data: + :return dict: + """ + stacked_data = {k: [] for k in data[0]} + for _, row in enumerate(data): + for key in stacked_data: + stacked_data[key].append(row[key]) + for key, content in stacked_data.items(): + stacked_data[key] = np.stack(content, axis=0) + return stacked_data + + +def read_records_before_error_section(epsfile_obj, mdr_descriptor, mdr_class_offset): + """FIXME DOC. + + :param _io.BinaryIO epsfile_obj: + :param mdr_descriptor: + :param mdr_class_offset: + :return dict: + """ + data = {} + for _, row in mdr_descriptor.iterrows(): + epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) + reshape_flag = False if row["FIELD"] in ["EARTH_LOCATION", "ANGULAR_RELATION"] else True + values = read_values(epsfile_obj, row, reshape=reshape_flag) + data[row["FIELD"]] = values + return data + + +def datetime_to_second_since_2000(date): + """FIXME DOC. + + :param datetime.datetime date: + :return float: + """ + era = datetime.datetime(2000, 1, 1) + return (date - era).total_seconds() + + +def read_nerr_values(epsfile_obj, mdr_descriptor, mdr_class_offset): + """FIXME DOC. + + Return a numpy array of all NERR values as red from each MDR class (i.e. data row). The + maximum values of NERR is used as second dimension of the errors variables (with the exclusion + of SURFACE_Z variable), so it is necessary read all values to reckon the max. + + :param epsfile_obj: + :param mdr_descriptor: + :param mdr_class_offset: + :return: + """ + nerr_row = mdr_descriptor.loc[mdr_descriptor["FIELD"] == "NERR"].iloc[0] + nerr_values = [] + epsfile_obj.seek(mdr_class_offset) + while True: + grh = epsnative_reader.grh_reader(epsfile_obj) + if grh: + if grh[0:3] == ("mdr", 1, 4): + epsfile_obj.seek(mdr_class_offset + int(nerr_row["OFFSET"])) + nerr_values.append(read_values(epsfile_obj, nerr_row, reshape=False)[0]) + mdr_class_offset += grh[3] + epsfile_obj.seek(mdr_class_offset) + else: + break + return np.array(nerr_values) + + +def read_all_rows(epsfile_obj, descriptor, mdr_class_offset, sensing_start, sensing_stop): + """FIXME DOC. + + :param _io.BinaryIO epsfile_obj: + :param dict descriptor: + :param int mdr_class_offset: + :param datetime.datetime sensing_start: + :param datetime.datetime sensing_stop: + :return (list, list, list): + """ + mdr_descriptor = descriptor[("mdr", 1, 4)] + last_constant_row = mdr_descriptor.loc[mdr_descriptor["FIELD"] == "ERROR_DATA_INDEX"].index[0] + mdr_descr_constant_offsets = mdr_descriptor[1:last_constant_row + 1] + vars_must_be_extend_to_fov = get_vars_must_be_extend_to_fov(mdr_descriptor) + max_nerr = read_nerr_values(epsfile_obj, mdr_descriptor, mdr_class_offset).max() + algorithms_data = [] + data_before_errors_section = [] + errors_data = [] + epsfile_obj.seek(mdr_class_offset) + while True: + grh = epsnative_reader.grh_reader(epsfile_obj) + if grh: + if grh[0:3] == ("mdr", 1, 4): + overlap = reckon_overlap( + sensing_start, sensing_stop, grh[-2], grh[-1] + ) + if overlap > 0: + record_start_time = datetime_to_second_since_2000(grh[-2]) + record_stop_time = datetime_to_second_since_2000(grh[-1]) + + records_before_error_section = read_records_before_error_section( + epsfile_obj, mdr_descr_constant_offsets, mdr_class_offset + ) + records_before_error_section["record_start_time"] = record_start_time + records_before_error_section["record_stop_time"] = record_stop_time + data_before_errors_section.append(records_before_error_section) + algorithm_data, errors = read_algorithm_sections( + epsfile_obj, + mdr_descriptor, + mdr_class_offset, + vars_must_be_extend_to_fov, + max_nerr, + ) + algorithms_data.append(algorithm_data) + errors_data.append(errors) + mdr_class_offset += grh[3] + else: + mdr_class_offset += grh[3] + else: + algorithms_data.append("dummy_mdr") + epsfile_obj.seek(mdr_class_offset) + else: + break + return data_before_errors_section, algorithms_data, errors_data + + +def read_product_data(epsfile_obj, descriptor, mdr_class_offset, sensing_start, sensing_stop): + """FIXME DOC. + + :param _io.BinaryIO epsfile_obj: + :param descriptor: + :param mdr_class_offset: + :param datetime.datetime sensing_start: + :param datetime.datetime sensing_stop: + :return (dict, dict, dict): + """ + data_before_errors, algorithms_data, errors = read_all_rows( + epsfile_obj, descriptor, mdr_class_offset, sensing_start, sensing_stop + ) + algorithms_data = fill_invalid_rows(algorithms_data) + stacked_algo_data = initialize_stacked_algorithms_output(algorithms_data) + stacked_algo_data = stack_algorithm_sections_along_rows(algorithms_data, stacked_algo_data) + stacked_data_before_errors = stack_non_algorithm_data(data_before_errors) + stacked_errors_data = stack_non_algorithm_data(errors) + return stacked_data_before_errors, stacked_algo_data, stacked_errors_data + + +def add_angular_relations(data_before_errors_section): + """FIXME DOC. + + :param data_before_errors_section: + :return: + """ + angular_relation = data_before_errors_section.pop("ANGULAR_RELATION") + solar_zenith = angular_relation[:, ::4] + satellite_zenith = angular_relation[:, 1::4] + solar_azimuth = angular_relation[:, 2::4] + satellite_azimuth = angular_relation[:, 3::4] + data_before_errors_section["solar_zenith"] = solar_zenith + data_before_errors_section["satellite_zenith"] = satellite_zenith + data_before_errors_section["solar_azimuth"] = solar_azimuth + data_before_errors_section["satellite_azimuth"] = satellite_azimuth + return data_before_errors_section + + +def add_latitude_longitude(data_before_errors_section): + """FIXME DOC. + + :param data_before_errors_section: + :return: + """ + earth_location = data_before_errors_section.pop("EARTH_LOCATION") + data_before_errors_section["lat"] = earth_location[:, ::2] + data_before_errors_section["lon"] = earth_location[:, 1::2] + return data_before_errors_section + + +def get_var_dimensions(var_name, values, dimensions): + """FIXME DOC. + + :param str var_name: + :param numpy.ndarray values: + :param dict dimensions: + :return dict: + """ + # nlt, nlq, nlo are all equals to 101, so the automatic detection fails + dims_keys = list(dimensions.keys()) + sizes = list(dimensions.values()) + dims = {} + for idx, k in enumerate(values.shape): + if idx == 0: + dims["along_track"] = dimensions["along_track"] + elif idx == 1: + if var_name in var_2nd_dim: + dims[var_2nd_dim[var_name]] = k + else: + dims["across_track"] = k + else: + if var_name in var_3rd_dim: + dims[var_3rd_dim[var_name]] = k + else: + dim_name = dims_keys[sizes.index(k)] + dims[dim_name] = k + return dims + + +def add_giadr_variables(dataset, giadr, band_to_record, dimensions, coordinates): + """FIXME DOC.""" + var_to_read = { + "pressure_levels_humidity": "nlq", + "pressure_levels_ozone": "nlo", + "pressure_levels_temp": "nlt", + "surface_emissivity_wavelengths": "new", + } + for var_name, dim_name in var_to_read.items(): + dims = {dim_name: dimensions[dim_name]} + coords = {k: coordinates[k] for k in dims if k in coordinates} + dataset[var_name] = xr.DataArray( + giadr[var_name.upper()]["values"], dims=dims, coords=coords + ) + metadata = band_to_record.get(var_name, {}).get("metadata", {}) + if metadata.get("scale_factor"): + metadata["scale_factor"] = 10 ** -metadata["scale_factor"] + dataset[var_name].attrs = metadata + return + + +def add_errors_variables(dataset, errors_data, band_to_record, dimensions, coordinates): + """FIXME DOC.""" + for var, values in errors_data.items(): + var_name = var.lower() + values = np.rollaxis(values, -1, 1) + dims = get_var_dimensions(var_name, values, dimensions) + coords = {k: coordinates[k] for k in dims if k in coordinates} + dataset[var_name] = xr.DataArray(values, dims=dims, coords=coords) + metadata = band_to_record.get(var_name, {}).get("metadata", {}) + if metadata.get("scale_factor"): + metadata["scale_factor"] = 10 ** -metadata["scale_factor"] + dataset[var_name].attrs = metadata + return + + +def assemble_dimensions(giadr, nr_rows, max_nerr): + """See Parameter Table (8.6) in EUM/OPS-EPS/MAN/04/0033. + + :param dict giadr: + :param int nr_rows: + :param int max_nerr: + :return dict: + """ + dimensions = { + "across_track": NR_FOV, + "along_track": nr_rows, + "cloud_formations": 3, + "co_nbr": NR_FOV, # pad to FOV number; real value for CO_NBR is provided in MDR + "nerr": max_nerr, # max value of NERR values provided in MDR + "nerro": giadr["NUM_OZONE_PCS"]["values"][0] + * (giadr["NUM_OZONE_PCS"]["values"][0] + 1) + / 2, + "nerrt": giadr["NUM_TEMPERATURE_PCS"]["values"][0] + * (giadr["NUM_TEMPERATURE_PCS"]["values"][0] + 1) + / 2, + "nerrw": giadr["NUM_WATER_VAPOUR_PCS"]["values"][0] + * (giadr["NUM_WATER_VAPOUR_PCS"]["values"][0] + 1) + / 2, + "neva_co": round(giadr["FORLI_NUM_LAYERS_CO"]["values"][0] / 2), + "neve_co": round(giadr["FORLI_NUM_LAYERS_CO"]["values"][0] / 2) + * giadr["FORLI_NUM_LAYERS_CO"]["values"][0], + "new": giadr["NUM_SURFACE_EMISSIVITY_WAVELENGTHS"]["values"][0], + "nl_co": giadr["FORLI_NUM_LAYERS_CO"]["values"][0], + "nl_hno3": giadr["FORLI_NUM_LAYERS_HNO3"]["values"][0], + "nl_o3": giadr["FORLI_NUM_LAYERS_O3"]["values"][0], + "nl_so2": giadr["BRESCIA_NUM_ALTITUDES_SO2"]["values"][0], + "nlo": giadr["NUM_PRESSURE_LEVELS_OZONE"]["values"][0], + "nlq": giadr["NUM_PRESSURE_LEVELS_HUMIDITY"]["values"][0], + "nlt": giadr["NUM_PRESSURE_LEVELS_TEMP"]["values"][0], + "npco": giadr["NUM_OZONE_PCS"]["values"][0], + "npct": giadr["NUM_TEMPERATURE_PCS"]["values"][0], + "npcw": giadr["NUM_WATER_VAPOUR_PCS"]["values"][0], + } + return dimensions + + +def add_variable_to_nc(ds, var_name, values, dimensions, coordinates, band_to_record): + """FIXME DOC.""" + dims = get_var_dimensions(var_name, values, dimensions) + coords = {k: coordinates[k] for k in dims if k in coordinates} + ds[var_name] = xr.DataArray(values, dims=dims, coords=coords) + metadata = band_to_record.get(var_name, {}).get("metadata", {}) + if metadata.get("scale_factor") is not None: + metadata["scale_factor"] = 10.0 ** -metadata["scale_factor"] + ds[var_name].attrs = metadata + return ds + + +def create_netcdf_dataset(data_before_errors_section, algorithms_data, errors_data, giadr): + """FIXME DOC. + + :param data_before_errors_section: + :param algorithms_data: + :param errors_data: + :param giadr: + :return: + """ + band_to_record = epsnative_reader.bands_to_records_reader("IASISND02") + add_latitude_longitude(data_before_errors_section) + add_angular_relations(data_before_errors_section) + nr_rows = data_before_errors_section["lat"].shape[0] + max_nerr = data_before_errors_section["NERR"].max() + dimensions = assemble_dimensions(giadr, nr_rows, max_nerr) + coordinates = { + "new": np.arange(giadr["NUM_SURFACE_EMISSIVITY_WAVELENGTHS"]["values"]), + "nl_co": np.arange(giadr["FORLI_NUM_LAYERS_CO"]["values"]), + "nl_hno3": np.arange(giadr["FORLI_NUM_LAYERS_HNO3"]["values"]), + "nl_o3": np.arange(giadr["FORLI_NUM_LAYERS_O3"]["values"]), + "nl_so2": np.arange(giadr["BRESCIA_NUM_ALTITUDES_SO2"]["values"]), + "nlo": np.arange(giadr["NUM_PRESSURE_LEVELS_OZONE"]["values"]), + "nlq": np.arange(giadr["NUM_PRESSURE_LEVELS_HUMIDITY"]["values"]), + "nlt": np.arange(giadr["NUM_PRESSURE_LEVELS_TEMP"]["values"]), + "npco": np.arange(giadr["NUM_OZONE_PCS"]["values"]), + "npct": np.arange(giadr["NUM_TEMPERATURE_PCS"]["values"]), + "npcw": np.arange(giadr["NUM_WATER_VAPOUR_PCS"]["values"]), + } + ds = xr.Dataset() + for var, values in data_before_errors_section.items(): + var_name = var.lower().replace("flg_", "flag_").replace("nerr", "nerr_values") + if values.ndim > 2: + values = np.rollaxis(values, -1, 1) + add_variable_to_nc(ds, var_name, values, dimensions, coordinates, band_to_record) + for _, content in algorithms_data.items(): + for var, values in content.items(): + var_name = f"{var.lower()}_values" if var == "CO_NBR" else var.lower() + if values.ndim > 2: + values = np.rollaxis(values, -1, 1) + add_variable_to_nc(ds, var_name, values, dimensions, coordinates, band_to_record) + add_giadr_variables(ds, giadr, band_to_record, dimensions, coordinates) + add_errors_variables(ds, errors_data, band_to_record, dimensions, coordinates) + ds["cloud_formation"] = xr.DataArray(np.arange(3), dims={"cloud_formation": 3}) + ds["cloud_formation"].attrs = band_to_record["cloud_formation"]["metadata"] + return ds + + +def reckon_overlap(sensing_start, sensing_end, class_start, class_stop): + """FIXME DOC. + + :param datetime.datetime sensing_start: + :param datetime.datetime sensing_end: + :param datetime.datetime class_start: + :param datetime.datetime class_stop: + :return float: + """ + latest_start = max(sensing_start, class_start) + earliest_end = min(sensing_end, class_stop) + delta = (earliest_end - latest_start).total_seconds() + overlap = max(0.0, delta) + return overlap From 78c2bb533b361cb284aeab3d0893ee1a71809f98 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Thu, 6 Jul 2023 15:15:13 +0200 Subject: [PATCH 10/23] Added tests for IASI L2 SND in EPS format --- .../IASISND02/band_to_record.yaml | 790 ++++++++++++++++++ .../IASISND02/general_metadata.yaml | 25 + .../eps_native_format/IASISND02/giadr_1_4.csv | 26 + .../eps_native_format/IASISND02/mdr_1_4.csv | 94 +++ .../readers/eps_native_format/mphr_0_2.csv | 75 ++ satpy/readers/epsnative_reader.py | 12 + satpy/tests/reader_tests/test_iasisnd02.py | 131 +++ 7 files changed, 1153 insertions(+) create mode 100644 satpy/etc/readers/eps_native_format/IASISND02/band_to_record.yaml create mode 100644 satpy/etc/readers/eps_native_format/IASISND02/general_metadata.yaml create mode 100644 satpy/etc/readers/eps_native_format/IASISND02/giadr_1_4.csv create mode 100644 satpy/etc/readers/eps_native_format/IASISND02/mdr_1_4.csv create mode 100644 satpy/etc/readers/eps_native_format/mphr_0_2.csv create mode 100644 satpy/tests/reader_tests/test_iasisnd02.py diff --git a/satpy/etc/readers/eps_native_format/IASISND02/band_to_record.yaml b/satpy/etc/readers/eps_native_format/IASISND02/band_to_record.yaml new file mode 100644 index 0000000000..a2579991eb --- /dev/null +++ b/satpy/etc/readers/eps_native_format/IASISND02/band_to_record.yaml @@ -0,0 +1,790 @@ +atmospheric_ozone: + name: atmospheric_ozone + record_name: ATMOSPHERIC_OZONE + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + standard_name: "mass_concentration_of_ozone_in_air" + # coordinates: "lat lon pressure_levels_ozone" + comment: "Ozone (for 120 IFOV with up to 101 vertical levels)" + units: "kg/kg" + scale_factor: 8 +atmospheric_temperature: + name: atmospheric_temperature + record_name: ATMOSPHERIC_TEMPERATURE + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + standard_name: "air_temperature" + # coordinates: "lat lon pressure_levels_temp" + comment: "Temperature (for 120 IFOV with up to 101 vertical levels)" + units: "degree(K)" + scale_factor: 2 +atmospheric_water_vapour: + name: atmospheric_water_vapour + record_name: ATMOSPHERIC_WATER_VAPOUR + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "water vapour" + # coordinates: "lat lon pressure_levels_humidity" + units: "kg/kg" + scale_factor: 7 +cloud_formation: + name: cloud_formation + record_name: NUMBER_CLOUD_FORMATIONS + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "number of cloud formations" + # coordinates: "lat lon pressure_levels_humidity" + positive: "down" + definition: "0: CLOUD_FORMATION_1, 1: CLOUD_FORMATION_2, 2: CLOUD_FORMATION_3," +cloud_phase: + name: cloud_phase + record_name: CLOUD_PHASE + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "cloud phase" + # coordinates: "lat lon cloud_formation" + comment: "Cloud Phase (for 120 IFOV with up to 3 cloud formations) (0 = no cloud, 1 = liquid, 2 = ice, 3 = mixed, 255 = undefined)" + positive: "down" + definition: "0: CLOUD_FORMATION_1, 1: CLOUD_FORMATION_2, 2: CLOUD_FORMATION_3," + units: "1" + flag_meaning: "no cloud, liquid, ice, mixed, undefined" + scale_factor: 0 +cloud_top_pressure: + name: cloud_top_pressure + record_name: CLOUD_TOP_PRESSURE + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "cloud top pressure" + # coordinates: "lat lon cloud_formation" + comment: "Cloud top pressure (for 120 IFOV with up to 3 cloud formations)" + units: "Pa" + scale_factor: 0 +cloud_top_temperature: + name: cloud_top_temperature + record_name: CLOUD_TOP_TEMPERATURE + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "cloud top temperature" + # coordinates: "lat lon cloud_formation" + comment: "Cloud top temperature (for 120 IFOV with up to 3 cloud formations)" + units: "K" + scale_factor: 2 +co_bdiv: + name: co_bdiv + record_name: CO_BDIV + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Retrieval flags " + comment: "Retrieval flags " + units: "1" +co_cp_air: + name: co_cp_air + record_name: CO_CP_AIR + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Air partial columns on each retrieved layer " + comment: "Air partial columns on each retrieved layer" + units: "molecules/cm2" + scale_factor: -20 + missing_values: -2 +co_cp_co_a: + name: co_cp_co_a + record_name: CO_CP_CO_A + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "A-priori partial columns for CO en each retrieved layer " + comment: "A-priori partial columns for CO en each retrieved layer " + units: "molecules/cm2" + scale_factor: -13 + missing_value: -2 +co_h_eigenvalues: + name: co_h_eigenvalues + record_name: CO_H_EIGENVALUES + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Main eigenvalues of the sensitivity matrix " + comment: "Main eigenvalues of the sensitivity matrix " + units: "1" + missing_value: 9.9999998E12 +co_h_eigenvectors: + name: co_h_eigenvectors + record_name: CO_H_EIGENVECTORS + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Main eigenvectors of the sensitivity matrix " + comment: "Main eigenvectors of the sensitivity matrix " + units: "1" + missing_value: 9.9999998E12 +co_nbr_values: + name: co_nbr_values + record_name: CO_NBR + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Number of CO profiles retrieved in scanline " + units: "1" +co_nfitlayers: + name: co_nfitlayers + record_name: CO_NFITLAYERS + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Number of CO profiles retrieved in scanline" + units: "1" + comment: "Number of CO profiles retrieved in scanline" +co_npca: + name: co_npca + record_name: CO_NPCA + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Number of vectors describing the characterization matrices " + units: "1" + comment: "Number of vectors describing the characterization matrices " +co_qflag: + name: co_qflag + record_name: CO_QFLAG + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Air partial columns on each retrieved layer " + units: "1" + comment: "General retrieval quality flag" +co_x_co: + name: co_x_co + record_name: CO_QFLAG + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Scaling vector multiplying the a-priori CO vector in order to define the retrieved CO vector. " + units: "1" + comment: "Scaling vector multiplying the a-priori CO vector in order to define the retrieved CO vector. " + missing_value: 9.9999998E12 +degraded_inst_mdr: + name: degraded_ins_MDR + record_name: DEGRADED_INST_MDR + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Quality of MDR has been degraded due to an instrument degradation " + units: "1" +degraded_proc_mdr: + name: degraded_proc_MDR + record_name: DEGRADED_PROC_MDR + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Quality of MDR has been degraded due to an instrument degradation " + units: "1" +error_data_index: + name: error_data_index + record_name: ERROR_DATA_INDEX + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Index of the error data record corresponding to the IFOVs in the line " + comment: "Index of the error data record corresponding to the IFOVs in the line " + units: "1" +fg_atmospheric_ozone: + name: fg_atmospheric_ozone + record_name: FG_QI_ATMOSPHERIC_OZONE + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "apriori ozone" + units: "kg.m^-2" + # coordinates: "lat lon pressure_levels_ozone" + comment: "A-priori ozone profile (for 30 EFOV with up to 101 vertical levels)" + scale_factor: 8 +fg_atmospheric_temperature: + name: fg_atmospheric_temperature + record_name: FG_ATMOSPHERIC_TEMPERATURE + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "apriori ozone" + units: "kg.m^-2" + # coordinates: "lat lon pressure_levels_ozone" + comment: "A-priori ozone profile (for 30 EFOV with up to 101 vertical levels)" + scale_factor: 2 +fg_atmospheric_water_vapour: + name: fg_atmospheric_water_vapour + record_name: FG_ATMOSPHERIC_WATER_VAPOUR + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "apriori water vapour" + units: "kg/kg" + # coordinates: "lat lon pressure_levels_humidity" + scale_factor: 7 +fg_surface_temperature: + name: fg_surface_temperature + record_name: FG_ATMOSPHERIC_WATER_VAPOUR + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "surface temperature" + units: "K" + # coordinates: "lat lon" + comment: "A-priori surface skin temperature" + scale_factor: 2 +flag_amsubad: + name: flag_amsubad + record_name: FLG_AMSUBAD + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Availability and quality of AMSU measurements" + units: "1" + flag_meaning: "enumerated flag" +flag_avhrrbad: + name: flag_avhrrbad + record_name: FLG_AVHRRBAD + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Availability and quality of AVHRR measurements" + units: "1" + flag_meaning: "enumerated flag" +flag_cldfrm: + name: flag_cldfrm + record_name: FLG_CLDFRM + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Origin of characterisation of the cloud formations" + units: "1" + flag_meaning: "bit field flag" +flag_cldtst: + name: flag_cldtst + record_name: FLG_CLDTST + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Details of cloud tests exectued and their results" + units: "1" + flag_meaning: "bit field flag" +flag_cldnes: + name: flag_cldnes + record_name: FLG_CLDNES + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Cloudiness assessment summary" + units: "1" + flag_meaning: "enumerated flag" +flag_daynit: + name: flag_daynit + record_name: FLG_DAYNIT + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + flag_values: "0b, 1b, 2b" + flag_meaning: "0:day, 1:night, 2:twilight" + long_name: "Discrimination between day and night" + units: "1" +flag_dustcld: + name: flag_dustcld + record_name: FLG_DUSTCLD + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Indicates presence of dust clouds in the IFOV" + units: "1" +flag_fgcheck: + name: flag_fgcheck + record_name: FLG_FGCHECK + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + flag_meaning: "bit field flag" + long_name: "Check that geophysical parameters from the first guess are within bounds" + units: "1" +flag_iasibad: + name: flag_iasibad + record_name: FLG_IASIBAD + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + flag_meaning: "bit field flag" + long_name: "Validation flag for IASI Level 1 product" + units: "1" +flag_initia: + name: flag_initia + record_name: FLG_INITIA + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + flag_meaning: "bit field flag" + long_name: "Indicates the measurements used in the first guess retrieval" + units: "1" +flag_itconv: + name: flag_itconv + record_name: FLG_ITCONV + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + flag_meaning: "0:Iteration did not converge-sounding rejected, 1:Iteration did not converge-sounding accepted, 2:Iteration converged-sounding accepted" + long_name: "Convergence of the iterative retrieval" + units: "1" +flag_lansea: + name: flag_lansea + record_name: FLG_LANSEA + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + flag_values: "0b, 1b, 2b, 3b, 4b" + flag_meaning: "0:water, 1:land low, 2:land high, 3:land water low, 4:land water high" + long_name: "Specifies surface type" + units: "1" +flag_mhsbad: + name: flag_mhsbad + record_name: FLG_MHSBAD + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Availability and quality of MHS measurements" + units: "1" +flag_numit: + name: flag_numit + record_name: FLG_NUMIT + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Number of iterations in the OEM" + units: "1" +flag_nwpbad: + name: flag_nwpbad + record_name: FLG_NWPBAD + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Availability and quality of NWP measurements" + units: "1" +flag_physcheck: + name: flag_physcheck + record_name: FLG_PHYSCHECK + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Indicates potential corrections for superadiabatic and supersaturation conditions" + units: "1" +flag_retcheck: + name: flag_retcheck + record_name: FLG_RETCHECK + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Check that geophysical parameters from the OEM are within bounds" + units: "1" +flag_satman: + name: flag_satman + record_name: FLG_SATMAN + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Indication of satellite manouevre" + units: "1" +flag_sunglnt: + name: flag_sunglnt + record_name: FLG_SUNGLNT + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Identification of sun glint" + units: "1" +flag_thicir: + name: flag_thicir + record_name: FLG_THICIR + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "Thin cirrus cloud test" + units: "1" +fractional_cloud_cover: + name: fractional_cloud_cover + record_name: FRACTIONAL_CLOUD_COVER + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "fractional cloud cover" + units: "%" + # coordinates: "lat lon cloud_formation" + standard_name: "cloud_area_fraction" + comment: "Fractional cloud cover (for 120 IFOV with up to 3 cloud formations)" + scale_factor: 2 +instrument_mode: + name: instrument_mode + record_name: INSTRUMENT_MODE + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "GEPSIasiMode Instrument Mode Not Used Flag" + units: "1" + flag_meaning: "false true" + comment: "Instrument mode. This is a copy of the MDR-1C flag GEPSIasiMode as defined in the IASI L1 PFS." +integrated_ch4: + name: integrated_ch4 + record_name: INTEGRATED_CH4 + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "integrated CH4" + standard_name: "atmosphere_mass_content_of_methane" + units: "kg.m^-2" + # coordinates: "lat lon" + comment: "Integrated CH4 (for 120 IFOV)" + scale_factor: 6 +integrated_co: + name: integrated_co + record_name: INTEGRATED_CO + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "integrated CO" + standard_name: "atmosphere_mass_content_of_carbon_monoxide" + units: "kg.m^-2" + # coordinates: "lat lon" + comment: "Integrated CO (for 120 IFOV)" + scale_factor: 7 +integrated_co2: + name: integrated_co2 + record_name: INTEGRATED_CO2 + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "integrated CO2" + standard_name: "atmosphere_mass_content_of_carbon_monoxide" + units: "kg.m^-2" + # coordinates: "lat lon" + comment: "Integrated CO2 (for 120 IFOV)" + scale_factor: 3 +integrated_n2o: + name: integrated_n2o + record_name: INTEGRATED_N2O + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "integrated N2O" + standard_name: "atmosphere_mass_content_of_carbon_monoxide" + units: "kg.m^-2" + # coordinates: "lat lon" + comment: "Integrated N2O (for 120 IFOV)" + scale_factor: 6 +integrated_ozone: + name: integrated_ozone + record_name: INTEGRATED_OZONE + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "integrated ozone" + standard_name: "atmosphere_mass_content_of_ozone" + units: "kg.m^-2" + # coordinates: "lat lon" + comment: "Integrated ozone (for 120 IFOV)" + scale_factor: 6 +integrated_water_vapour: + name: integrated_water_vapour + record_name: INTEGRATED_WATER_VAPOUR + mdr_id: mdr_1 + band_position: 0 + convention: BIL + metadata: + long_name: "integrated water vapour" + standard_name: "atmosphere_mass_content_of_water_vapour" + units: "kg.m^-2" + # coordinates: "lat lon" + comment: "Integrated ozone (for 120 IFOV)" + scale_factor: 2 +nerr_values: + name: nerr_values + record_name: NERR + mdr_id: mdr_1 + band_position: 0 + convention: BIP + metadata: + long_name: "Number of error data records for current scan line" + units: "1" +number_cloud_formations: + name: number_cloud_formations + record_name: NUMBER_CLOUD_FORMATIONS + mdr_id: mdr_1 + band_position: 0 + convention: BIP + metadata: + long_name: "number of cloud formations" + units: "1" + comment: "Number of cloud formations in IFOV" +ozone_error: + name: ozone_error + record_name: OZONE_ERROR + mdr_id: mdr_1 + band_position: 0 + convention: BIP + metadata: + long_name: "Retrieval error covariance matrix for ozone in principal component domain" + units: "1" + comment: "Retrieval error covariance matrix for ozone in principal component domain " +pressure_levels_humidity: + name: pressure_levels_humidity + record_name: pressure_levels_humidity + mdr_id: giadr_1_4 + band_position: 0 + convention: BIP + metadata: + long_name: "pressure levels on which humidity profiles are retrieved" + units: "Pa" + comment: "Pressure levels on which retrieved humidity profiles are given" + scale_factor: 2 +pressure_levels_ozone: + name: pressure_levels_ozone + record_name: PRESSURE_LEVELS_OZONE + mdr_id: giadr_1_4 + band_position: 1 + convention: BIP + metadata: + long_name: "pressure levels on which humidity profiles are retrieved" + units: "Pa" + comment: "Pressure levels on which retrieved humidity profiles are given" + scale_factor: 2 +pressure_levels_temp: + name: pressure_levels_temp + record_name: ressure_levels_temp + mdr_id: giadr_1_4 + band_position: 1 + convention: BIP + metadata: + long_name: "pressure levels on which humidity profiles are retrieved" + units: "Pa" + comment: "Pressure levels on which retrieved humidity profiles are given" + scale_factor: 2 +record_start_time: + metadata: + long_name: "Record start time from the Generic Record Header of the MDR" + units: "s" + scale_factor: 0 +record_stop_time: + metadata: + long_name: "Record stop time from the Generic Record Header of the MDR" + units: "s" + scale_factor: 0 +satellite_azimuth: + name: satellite_azimuth + record_name: ANGULAR_RELATION + mdr_id: mdr_1 + band_position: 4 + convention: BIL + metadata: + long_name: "angular relation for the earth view: satellite azimuth" + units: "degrees" + scale_factor: 2 +satellite_zenith: + name: satellite_zenith + record_name: ANGULAR_RELATION + mdr_id: mdr_1 + band_position: 2 + convention: BIL + metadata: + long_name: "angular relation for the earth view: satellite zenith" + units: "degrees" + scale_factor: 2 +solar_azimuth: + name: solar_azimuth + record_name: ANGULAR_RELATION + mdr_id: mdr_1 + band_position: 3 + convention: BIL + metadata: + long_name: "angular relation for the earth view: solar azimuth" + units: "degrees" + scale_factor: 2 +solar_zenith: + name: solar_zenith + record_name: ANGULAR_RELATION + mdr_id: mdr_1 + band_position: 1 + convention: BIL + metadata: + long_name: "angular relation for the earth view: solar zenith" + units: "degrees" + scale_factor: 2 +spacecraft_altitude: + name: spacecraft_altitude + record_name: SPACECRAFT_ALTITUDE + mdr_id: mdr_1 + band_position: 1 + convention: BIL + metadata: + comment: "Spacecraft Altitude Above Reference Geoid (MSL) " + long_name: "spacecraft altitude" + units: "km" + scale_factor: 1 +surface_emissivity: + name: surface_emissivity + record_name: SURFACE_EMISSIVITY + mdr_id: mdr_1 + band_position: 1 + convention: BIL + metadata: + long_name: "surface emissivity" + units: "1" + standard_name: "surface_longwave_emissivity" + # coordinates: "lat lon surface_emissivity_wavelengths" + comment: "Surface emissivity (for 120 IFOV with up to 20 wavelengths)" + scale_factor: 4 +surface_emissivity_wavelengths: + name: surface_emissivity_wavelengths + record_name: SURFACE_EMISSIVITY_WAVELENGTHS + mdr_id: giadr_1_4 + band_position: 1 + convention: BIL + metadata: + long_name: "wavelengths for surface emissivity" + units: "um" + standard_name: "surface_longwave_emissivity" + comment: "Wavelengths for surface emissivity" + scale_factor: 4 +surface_pressure: + name: surface_pressure + record_name: SURFACE_PRESSURE + mdr_id: mdr_1 + band_position: 1 + convention: BIL + metadata: + long_name: "surface pressure" + units: "Pa" + # coordinates: "lat lon" + standard_name: "surface_air_pressure" + comment: "Surface pressure" + scale_factor: 0 +surface_temperature: + name: surface_temperature + record_name: SURFACE_TEMPERATURE + mdr_id: mdr_1 + band_position: 1 + convention: BIL + metadata: + long_name: "surface temperature" + units: "K" + standard_name: "surface_temperature" + # coordinates: "lat lon" + comment: "Surface temperature (for 120 IFOV and up to 2 temperatures)" + scale_factor: 2 +surface_z: + name: surface_z + record_name: SURFACE_Z + mdr_id: mdr_1 + band_position: 1 + convention: BIL + metadata: + missing_value: 32767S + long_name: "Altitude of surface " + units: "m" + comment: "Altitude of surface " +temperature_error: + name: temperature_error + record_name: TEMPERATURE_ERROR + mdr_id: mdr_1 + band_position: 1 + convention: BIL + metadata: + missing_value: -2147483646 + long_name: "Retrieval error covariance matrix for temperature in principal component domain. " + units: "1" + comment: "Retrieval error covariance matrix for temperature in principal component domain" +water_vapour_error: + name: water_vapour_error + record_name: TEMPERATURE_ERROR + mdr_id: mdr_1 + band_position: 1 + convention: BIL + metadata: + missing_value: -2147483646 + long_name: "Retrieval error covariance matrix for water-vapour in principal component domain " + units: "1" + comment: "Retrieval error covariance matrix for water-vapour in principal component domain" +lat: + name: latitude + record_name: EARTH_LOCATION + mdr_id: mdr_1 + band_position: 1 + convention: BIP + metadata: + long_name: "grid_latitude" + standard_name: "latitude" + units: "degrees_north" + scale_factor: 4 +lon: + name: longitude + record_name: EARTH_LOCATION + mdr_id: mdr_1 + band_position: 2 + convention: BIP + metadata: + long_name: "grid_longitude" + standard_name: "longitude" + units: "degrees_east" + scale_factor: 4 diff --git a/satpy/etc/readers/eps_native_format/IASISND02/general_metadata.yaml b/satpy/etc/readers/eps_native_format/IASISND02/general_metadata.yaml new file mode 100644 index 0000000000..aa85672433 --- /dev/null +++ b/satpy/etc/readers/eps_native_format/IASISND02/general_metadata.yaml @@ -0,0 +1,25 @@ +# Copyright 2017-2022, European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +'title': "IASI 2b the Infrared Atmospheric Sounding Interferometer" +'title_short_name': "IASI L2" +'summary': "NA" +'comment': "Search for Infrared Atmospheric Sounding Interferometer in the references URL" +'keywords': "EUMETSAT, DATA CENTRE, EPS, IASI1C, NetCDF" +'reference_url': "https://navigator.eumetsat.int/product/EO:EUM:DAT:METOP:IASSND02" +'wmo_filename': "" +'data_format_type': "" +'producer_agency': "EUMETSAT" +'platform_type': "spacecraft" +"contents": "IASI Level 2 Measurements" diff --git a/satpy/etc/readers/eps_native_format/IASISND02/giadr_1_4.csv b/satpy/etc/readers/eps_native_format/IASISND02/giadr_1_4.csv new file mode 100644 index 0000000000..341e5a8fc4 --- /dev/null +++ b/satpy/etc/readers/eps_native_format/IASISND02/giadr_1_4.csv @@ -0,0 +1,26 @@ +FIELD,DESCRIPTION,SF,UNITS,DIM1,DIM2,DIM3,TYPE,TYPE SIZE,FIELD SIZE,OFFSET +RECORD_HEADER,Generic record header,0,,1,1,1,REC_HEAD,20,20,0 +GIADR CONTENTS,,,,,,,,,,20 +NUM_PRESSURE_LEVELS_TEMP,Number of pressure levels for temperature profile retrieval (NLT),0,NA,1,1,1,u-byte,1,1,20 +PRESSURE_LEVELS_TEMP,Pressure levels on which retrieved temperature profiles are given,2,Pa,101,1,1,u-integer4,4,404,21 +NUM_PRESSURE_LEVELS_HUMIDITY,Number of pressure levels for humidity profile retrieval (NLQ),0,NA,1,1,1,u-byte,1,1,425 +PRESSURE_LEVELS_HUMIDITY,Pressure levels on which retrieved humidity profiles are given,2,Pa,101,1,1,u-integer4,4,404,426 +NUM_PRESSURE_LEVELS_OZONE,Number of pressure levels for ozone profile retrieval (NLO),0,NA,1,1,1,u-byte,1,1,830 +PRESSURE_LEVELS_OZONE,Pressure levels on which retrieved ozone profiles are given,2,Pa,101,1,1,u-integer4,4,404,831 +NUM_SURFACE_EMISSIVITY_WAVELENGTHS,Number of wavelengths for surface emissivity retrieval (NEW),0,NA,1,1,1,u-byte,1,1,1235 +SURFACE_EMISSIVITY_WAVELENGTHS,Wavelengths for surface emissivity,4,μm,12,1,1,u-integer4,4,48,1236 +ERROR_DATA,,,,,,,,,,1284 +NUM_TEMPERATURE_PCS,Number of principal components for temperature in the ERROR_DATA (NPCT),0,NA,1,1,1,u-byte,1,1,1284 +NUM_WATER_VAPOUR_PCS,Number of principal components for water-vapour in the ERROR_DATA (NPCW),0,NA,1,1,1,u-byte,1,1,1285 +NUM_OZONE_PCS,Number of principal components for ozone in the ERROR_DATA (NPCO),0,NA,1,1,1,u-byte,1,1,1286 +FORLI,,,,,,,,,,1287 +FORLI_NUM_LAYERS_CO,Number of partial layers for CO (NL_CO),0,NA,1,1,1,u-byte,1,1,1287 +FORLI_LAYER_HEIGHTS_CO,CO partial layer heights ,0,m,19,1,1,u-integer2,2,38,1288 +FORLI_NUM_LAYERS_HNO3,Number of partial layers for HNO3 (NL_HNO3),0,NA,1,1,1,u-byte,1,1,1326 +FORLI_LAYER_HEIGHTS_HNO3,HNO3 partial layer heights ,0,m,19,1,1,u-integer2,2,38,1327 +FORLI_NUM_LAYERS_O3,Number of partial layers for O3 (NL_O3),0,NA,1,1,1,u-byte,1,1,1365 +FORLI_LAYER_HEIGHTS_O3,O3 partial layer heights ,0,m,40,1,1,u-integer2,2,80,1366 +BRESCIA,,,,,,,,,,1446 +BRESCIA_NUM_ALTITUDES_SO2,Number of estimated SO2 plume heights (NL_SO2),0,NA,1,1,1,u-byte,1,1,1446 +BRESCIA_ALTITUDES_SO2,Estimated SO2 plume heights,0,m,5,1,1,u-integer2,2,10,1447 +,,,,,,,,,,1457 diff --git a/satpy/etc/readers/eps_native_format/IASISND02/mdr_1_4.csv b/satpy/etc/readers/eps_native_format/IASISND02/mdr_1_4.csv new file mode 100644 index 0000000000..77b3ea10a6 --- /dev/null +++ b/satpy/etc/readers/eps_native_format/IASISND02/mdr_1_4.csv @@ -0,0 +1,94 @@ +FIELD,DESCRIPTION,SF,UNITS,DIM1,DIM2,DIM3,TYPE,TYPE SIZE,FIELD SIZE,OFFSET +RECORD_HEADER,Generic record header,,,1,1,1,REC_HEAD,20,20,0 +DEGRADED_INST_MDR,Quality of MDR has been degraded from nominal due to an instrument degradation,,,1,1,1,boolean,1,1,20 +DEGRADED_PROC_MDR,Quality of MDR has been degraded from nominal due to a processing degradation,,,1,1,1,boolean,1,1,21 +FG_ATMOSPHERIC_TEMPERATURE,A-priori temperature profile (for 120 FOV with up to 101 vertical levels),2,K,101,120,1,u-integer2,2,24240,22 +FG_ATMOSPHERIC_WATER_VAPOUR,A-priori water vapour profile (for 30 EFOV with up to 101 vertical levels),7,kg/kg,101,120,1,u-integer4,4,48480,24262 +FG_ATMOSPHERIC_OZONE,A-priori ozone profile (for 30 EFOV with up to 101 vertical levels),8,kg/kg,101,120,1,u-integer2,2,24240,72742 +FG_SURFACE_TEMPERATURE,A-priori surface skin temperature,2,K,120,1,1,u-integer2,2,240,96982 +FG_QI_ATMOSPHERIC_TEMPERATURE,Quality indicator for a-priori temperature profile,1,,120,1,1,u-byte,1,120,97222 +FG_QI_ATMOSPHERIC_WATER_VAPOUR,Quality indicator for a-priori water vapour profile,1,,120,1,1,u-byte,1,120,97342 +FG_QI_ATMOSPHERIC_OZONE,Quality indicator for a-priori ozone profile,1,,120,1,1,u-byte,1,120,97462 +FG_QI_SURFACE_TEMPERATURE,Quality indicator for a-priori surface skin temperature,1,,120,1,1,u-byte,1,120,97582 +ATMOSPHERIC_TEMPERATURE,Temperature (for 120 IFOV with up to 101 vertical levels),2,K,101,120,1,u-integer2,2,24240,97702 +ATMOSPHERIC_WATER_VAPOUR,Water vapour (for 120 IFOV with up to 101 vertical levels),7,kg/kg,101,120,1,u-integer4,4,48480,121942 +ATMOSPHERIC_OZONE,Ozone (for 120 IFOV with up to 101 vertical levels),8,kg/kg,101,120,1,u-integer2,2,24240,170422 +SURFACE_TEMPERATURE,Surface temperature (for 120 IFOV),2,K,120,1,1,u-integer2,2,240,194662 +INTEGRATED_WATER_VAPOUR,Integrated water vapour (for 120 IFOV),2,kg.m^-2,120,1,1,u-integer2,2,240,194902 +INTEGRATED_OZONE,Integrated ozone (for 120 IFOV),6,kg.m^-2,120,1,1,u-integer2,2,240,195142 +INTEGRATED_N2O,Integrated N20 (for 120 IFOV),6,kg.m^-2,120,1,1,u-integer2,2,240,195382 +INTEGRATED_CO,Integrated CO (for 120 IFOV),7,kg.m^-2,120,1,1,u-integer2,2,240,195622 +INTEGRATED_CH4,Integrated CH4 (for 120 IFOV),6,kg.m^-2,120,1,1,u-integer2,2,240,195862 +INTEGRATED_CO2,Integrated CO2 ( for 120 IFOV),3,kg.m^-2,120,1,1,u-integer2,2,240,196102 +SURFACE_EMISSIVITY,Surface emissivity (for 120 IFOV with up to 20 wavelengths),4,,12,120,1,u-integer2,2,2880,196342 +NUMBER_CLOUD_FORMATIONS,Number of cloud formations in IFOV,0,,120,1,1,u-byte,1,120,199222 +FRACTIONAL_CLOUD_COVER,Fractional cloud cover (for 120 IFOV with up to 3 cloud formations),2,%,3,120,1,u-integer2,2,720,199342 +CLOUD_TOP_TEMPERATURE,Cloud top temperature (for 120 IFOV with up to 3 cloud formations),2,K,3,120,1,u-integer2,2,720,200062 +CLOUD_TOP_PRESSURE,Cloud top pressure (for 120 IFOV with up to 3 cloud formations),0,Pa,3,120,1,u-integer4,4,1440,200782 +CLOUD_PHASE,"Cloud Phase (for 120 IFOV with up to 3 cloud formations) (0 = no cloud, 1 = liquid, 2 = ice, 3 = mixed, 255 = undefined)",0,,3,120,1,enumerated,1,360,202222 +SURFACE_PRESSURE,Surface pressure,0,Pa,120,1,1,u-integer4,4,480,202582 +INSTRUMENT_MODE,Instrument mode. This is a copy of the MDR-1C flag GEPSIasiMode as defined in the IASI L1 PFS.,0,,1,1,1,enumerated,1,1,203062 +SPACECRAFT_ALTITUDE,Spacecraft Altitude Above Reference Geoid (MSL) ,1,km,1,1,1,u-integer4,4,4,203063 +ANGULAR_RELATION,"Angular relationships: solar zenith angle, satellite zenith angle, solar azimuth angle, satellite azimuth angle for 120 IFOV ",2,deg,4,120,1,integer2,2,960,203067 +EARTH_LOCATION,"Earth Location: latitude, longitude of surface footprint (for 120 IFOV)",4,deg,2,120,1,integer4,4,960,204027 +FLG_AMSUBAD,Availability and quality of AMSU measurements,,,120,1,1,enumerated,1,120,204987 +FLG_AVHRRBAD,Availability and quality of AVHRR measurements,,,120,1,1,enumerated,1,120,205107 +FLG_CLDFRM,Origin of characterisation of the cloud formations,,,120,1,1,bitst(8),1,120,205227 +FLG_CLDNES,Cloudiness assessment summary,,,120,1,1,enumerated,1,120,205347 +FLG_CLDTST,Details of cloud tests exectued and their results,,,120,1,1,bitst(16),2,240,205467 +FLG_DAYNIT,Discrimination between day and night,,,120,1,1,enumerated,1,120,205707 +FLG_DUSTCLD,Indicates presence of dust clouds in the IFOV,,,120,1,1,enumerated,1,120,205827 +FLG_FGCHECK,Check that geophysical parameters from the first guess are within bounds,,,120,1,1,bitst(16),2,240,205947 +FLG_IASIBAD,Availability and quality of IASI L1 measurements,,,120,1,1,enumerated,1,120,206187 +FLG_INITIA,Indicates the measurements used in the first guess retrieval,,,120,1,1,bitst(8),1,120,206307 +FLG_ITCONV,Convergence and acceptance of the OEM result,,,120,1,1,enumerated,1,120,206427 +FLG_LANSEA,Specifies surface type,,,120,1,1,enumerated,1,120,206547 +FLG_MHSBAD,Availability and quality of MHS measurements,,,120,1,1,enumerated,1,120,206667 +FLG_NUMIT,Number of iterations in the OEM,0,,120,1,1,u-byte,1,120,206787 +FLG_NWPBAD,Availability and quality of NWP data,,,120,1,1,enumerated,1,120,206907 +FLG_PHYSCHECK,Indicates potential corrections for superadiabatic and supersaturation conditions,,,120,1,1,bitst(8),1,120,207027 +FLG_RETCHECK,Check that geophysical parameters from the OEM are within bounds,,,120,1,1,bitst(16),2,240,207147 +FLG_SATMAN,Indication of satellite manouevre,,,120,1,1,enumerated,1,120,207387 +FLG_SUNGLNT,Identification of sun glint,,,120,1,1,enumerated,1,120,207507 +FLG_THICIR,Thin cirrus cloud test,,,120,1,1,enumerated,1,120,207627 +NERR,Number of error data records for current scan line,0,,1,1,1,u-byte,1,1,207747 +ERROR_DATA_INDEX,Index of the error data record corresponding to the IFOVs in the line (=255 if N/A),,,120,1,1,u-byte,1,120,207748 +TEMPERATURE_ERROR,Retrieval error covariance matrix for temperature in principal component domain,,,406,nerr,1,bitst(32),4,, +WATER_VAPOUR_ERROR,Retrieval error covariance matrix for water-vapour in principal component domain,,,171,nerr,1,bitst(32),4,, +OZONE_ERROR,Retrieval error covariance matrix for ozone in principal component domain,,,55,nerr,1,bitst(32),4,, +SURFACE_Z,Altitude of surface ,0,m,120,1,1,integer2,2,240, +CO_QFLAG,General retrieval quality flag,,,120,1,1,enumerated,1,120, +CO_BDIV,Retrieval flags,,,120,1,1,bitst(32),4,480, +CO_NPCA,Number of vectors describing the characterization matrices,0,,120,1,1,u-byte,1,120, +CO_NFITLAYERS,Number of layers actually retrieved,0,,120,1,1,u-byte,1,120, +CO_NBR,Number of CO profiles retrieved in scanline,0,,1,1,1,u-byte,1,1, +CO_CP_AIR,Air partial columns on each retrieved layer,-20,molecules/cm2,19,co_nbr,1,u-integer2,2,, +CO_CP_CO_A,A-priori partial columns for CO en each retrieved layer,-13,molecules/cm2,19,co_nbr,1,u-integer2,2,, +CO_X_CO,Scaling vector multiplying the a-priori CO vector in order to define the retrieved CO vector. ,,,19,co_nbr,1,v-uinteger2,3,, +CO_H_EIGENVALUES,Main eigenvalues of the sensitivity matrix,,,10,co_nbr,1,v-integer4,5,, +CO_H_EIGENVECTORS,Main eigenvectors of the sensitivity matrix,,,190,co_nbr,1,v-integer4,5,, +HNO3_QFLAG,General retrieval quality flag,,,120,1,1,enumerated,1,120, +HNO3_BDIV,Retrieval flags,,,120,1,1,bitst(32),4,480, +HNO3_NPCA,Number of vectors describing the characterization matrices,0,,120,1,1,u-byte,1,120, +HNO3_NFITLAYERS,Number of layers actually retrieved,0,,120,1,1,u-byte,1,120, +HNO3_NBR,Number of HNO3 profiles retrieved in scanline,0,,1,1,1,u-byte,1,1, +HNO3_CP_AIR,Air partial columns on each retrieved layer,-20,molecules/cm2,19,hno3_nbr,1,u-integer2,2,, +HNO3_CP_HNO3_A,A-priori partial columns for HNO3 in each retrieved layer,-11,molecules/cm2,19,hno3_nbr,1,u-integer2,2,, +HNO3_X_HNO3,Scaling vector multiplying the a-priori HNO3 vector in order to define the retrieved HNO3 vector. ,,,19,hno3_nbr,1,v-uinteger2,3,, +HNO3_H_EIGENVALUES,Main eigenvalues of the sensitivity matrix,,,10,hno3_nbr,1,v-integer4,5,, +HNO3_H_EIGENVECTORS,Main eigenvectors of the sensitivity matrix,,,190,hno3_nbr,1,v-integer4,5,, +O3_QFLAG,General retrieval quality flag,,,120,1,1,enumerated,1,120, +O3_BDIV,Retrieval flags,,,120,1,1,bitst(32),4,480, +O3_NPCA,Number of vectors describing the characterization matrices,0,,120,1,1,u-byte,1,120, +O3_NFITLAYERS,Number of layers actually retrieved,0,,120,1,1,u-byte,1,120, +O3_NBR,Number of O3 profiles retrieved in scanline,0,,1,1,1,u-byte,1,1, +O3_CP_AIR,Air partial columns on each retrieved layer,-20,molecules/cm2,40,o3_nbr,1,u-integer2,2,, +O3_CP_O3_A,A-priori partial columns for O3 en each retrieved layer,-14,molecules/cm2,40,o3_nbr,1,u-integer2,2,, +O3_X_O3,Scaling vector multiplying the a-priori O3 vector in order to define the retrieved O3 vector. ,,,40,o3_nbr,1,v-uinteger2,3,, +O3_H_EIGENVALUES,Main eigenvalues of the sensitivity matrix,,,20,o3_nbr,1,v-integer4,5,, +O3_H_EIGENVECTORS,Main eigenvectors of the sensitivity matrix,,,800,o3_nbr,1,v-integer4,5,, +SO2_QFLAG,General retrieval quality flag,,,120,1,1,enumerated,1,120, +SO2_COL_AT_ALTITUDES,SO2 column for a plume at different estimated altitudes,1,DU,5,120,1,u-integer2,2,1200, +SO2_ALTITUDE,Retrieved plume altitude ,0,m,120,1,1,u-integer2,2,240, +SO2_COL,SO2 column at the retrieved plume altitude from an OEM approach,1,DU,120,1,1,u-integer2,2,240, +SO2_BT_DIFFERENCE,Indicative brightness temperature difference,2,K,120,1,1,integer2,2,240, diff --git a/satpy/etc/readers/eps_native_format/mphr_0_2.csv b/satpy/etc/readers/eps_native_format/mphr_0_2.csv new file mode 100644 index 0000000000..b636644c01 --- /dev/null +++ b/satpy/etc/readers/eps_native_format/mphr_0_2.csv @@ -0,0 +1,75 @@ +FIELD,DESCRIPTION,SF,UNITS,DIM1,DIM2,DIM3,DIM4,TYPE,TYPE_SIZE,FIELD_SIZE,OFFSET +PRODUCT_NAME ,Complete name of the product ,,,,,,,CHAR ,67,100,20 +PARENT_PRODUCT_NAME_1 ,"Name of the parent product from which this product has been produced. For Level 0 products, this field is filled with lower case x's. ",,,1,1,1,1,CHAR ,67,100,120 +PARENT_PRODUCT_NAME_2 ,"Name of the parent product from which this product has been produced. For Level 0 products or products for which this is not appropriate, this field is filled with lower case x's. ",,,1,1,1,1,CHAR ,67,100,220 +PARENT_PRODUCT_NAME_3 ,"Name of the parent product from which this product has been produced. For Level 0 products or products for which this is not appropriate, this field is filled with lower case x's. ",,,1,1,1,1,CHAR ,67,100,320 +PARENT_PRODUCT_NAME_4 ,"Name of the parent product from which this product has been produced. For Level 0 products or products for which this is not appropriate, this field is filled with lower case x's. ",,,1,1,1,1,CHAR ,67,100,420 +INSTRUMENT_ID ,Instrument identification ,,,1,1,1,1,E-CHAR ,4,37,520 +INSTRUMENT_MODEL ,Instrument Model identification ,,,1,1,1,1,ENUMERATED ,3,36,557 +PRODUCT_TYPE ,Product Type ,,,1,1,1,1,E-CHAR ,3,36,593 +PROCESSING_LEVEL ,Processing Level Identification ,,,1,1,1,1,E-CHAR ,2,35,629 +SPACECRAFT_ID ,Spacecraft identification ,,,1,1,1,1,E-CHAR ,3,36,664 +SENSING_START ,"UTC start time of sensing data in this object (PDU, ROI or Full Product) of original product as planned by the PROCESSING CENTRE ",,,1,1,1,1,GENERAL TIME ,15,48,700 +SENSING_END ,"UTC end time of sensing data in this object (PDU, ROI or Full Product) of original product as planned by the PROCESSING CENTRE ",,,1,1,1,1,GENERAL TIME ,15,48,748 +SENSING_START_THEORETICAL ,Theoretical UTC Time of start of sensing data in the dump from which this object is derived. This data is the predicted start time at the MPF level. ,,,1,1,1,1,GENERAL TIME ,15,48,796 +SENSING_END_THEORETICAL ,Theoretical UTC Time of end of sensing data in the dump from which this object is derived. This data is the predicted end time at the MPF level. ,,,1,1,1,1,GENERAL TIME ,15,48,844 +PROCESSING_CENTRE ,Processing Centre Identification ,,,1,1,1,1,E-CHAR ,4,37,892 +PROCESSOR_MAJOR_VERSION ,Processing chain major version number ,,,1,1,1,1,U-INTEGER ,5,38,929 +PROCESSOR_MINOR_VERSION ,Processing chain minor version number ,,,1,1,1,1,U-INTEGER ,5,38,967 +FORMAT_MAJOR_VERSION ,Dataset Format Major Version number ,,,1,1,1,1,U-INTEGER ,5,38,1005 +FORMAT_MINOR_VERSION ,Dataset Format Minor Version number ,,,1,1,1,1,U-INTEGER ,5,38,1043 +PROCESSING_TIME_START ,UTC time of the processing at start of processing for the product ,,,1,1,1,1,GENERAL TIME ,15,48,1081 +PROCESSING_TIME_END ,UTC time of the processing at end of processing for the product ,,,1,1,1,1,GENERAL TIME ,15,48,1129 +PROCESSING_MODE ,Identification of the mode of processing ,,,1,1,1,1,E-CHAR ,1,34,1177 +DISPOSITION_MODE ,Identification of the diposition mode ,,,1,1,1,1,E-CHAR ,1,34,1211 +RECEIVING_GROUND_STATION ,Acquisition Station Identification ,,,1,1,1,1,E-CHAR ,3,36,1245 +RECEIVE_TIME_START ,UTC time of the reception at CDA for first Data Item ,,,1,1,1,1,GENERAL TIME ,15,48,1281 +RECEIVE_TIME_END ,UTC time of the reception at CDA for last Data Item ,,,1,1,1,1,GENERAL TIME ,15,48,1329 +ORBIT_START ,"Start Orbit Number, counted incrementally since launch. Determined at time of SENSING START THEORETICAL ",,,1,1,1,1,U-INTEGER ,5,38,1377 +ORBIT_END ,Stop Orbit Number. Determined at time of SENSING_END_THEORETICAL ,,,1,1,1,1,U-INTEGER ,5,38,1415 +ACTUAL_PRODUCT_SIZE ,Size of the complete product ,,bytes ,1,1,1,1,U-INTEGER ,11,44,1453 +STATE_VECTOR_TIME ,Epoch time (in UTC) of the orbital elements and the orbit state vector. This corresponds to the time of crossing the ascending node for ORBIT START ,,UTC ,1,1,1,1,LONG GENERAL TIME ,18,51,1497 +SEMI_MAJOR_AXIS ,Semi major axis of orbit at time of the ascending node crossing [TRUE-OF-DATE] ,,mm ,1,1,1,1,INTEGER ,11,44,1548 +ECCENTRICITY ,Orbit eccentricity at time of the ascending node crossing [TRUE-OF-DATE] ,6,,1,1,1,1,INTEGER ,11,44,1592 +INCLINATION ,Orbit inclination at time of the ascending node crossing [TRUE-OF-DATE] ,3,deg ,1,1,1,1,INTEGER ,11,44,1636 +PERIGEE_ARGUMENT ,Argument of perigee at time of the ascending node crossing [TRUE-OF-DATE] ,3,deg ,1,1,1,1,INTEGER ,11,44,1680 +RIGHT_ASCENSION ,Right ascension at time of the ascending node crossing [TRUE-OF-DATE] ,3,deg ,1,1,1,1,INTEGER ,11,44,1724 +MEAN_ANOMALY ,Mean anomaly at time of the ascending node crossing [TRUE-OF-DATE] ,3,deg ,1,1,1,1,INTEGER ,11,44,1768 +X_POSITION ,X position of the orbit state vector in the orbit frame at ascending node [EARTH-FIXED] ,3,m ,1,1,1,1,INTEGER ,11,44,1812 +Y_POSITION ,Y position of the orbit state vector in the orbit frame at ascending node [EARTH-FIXED] ,3,m ,1,1,1,1,INTEGER ,11,44,1856 +Z_POSITION ,Z position of the orbit state vector in the orbit frame at ascending node [EARTH-FIXED] ,3,m ,1,1,1,1,INTEGER ,11,44,1900 +X_VELOCITY ,X velocity of the orbit state vector in the orbit frame at ascending node [EARTH-FIXED] ,3,m/s ,1,1,1,1,INTEGER ,11,44,1944 +Y_VELOCITY ,Y velocity of the orbit state vector in the orbit frame at ascending node [EARTH-FIXED] ,3,m/s ,1,1,1,1,INTEGER ,11,44,1988 +Z_VELOCITY ,Z velocity of the orbit state vector in the orbit frame at ascending node [EARTH-FIXED] ,3,m/s ,1,1,1,1,INTEGER ,11,44,2032 +EARTH_SUN_DISTANCE_RATIO ,"Earth-Sun distance ratio - ratio of current Earth-Sun distance to the Mean Earth-Sun distance which is 1Astronomical Unit or AU. +1 AU has a value of 1.495 978 706 91 x 10^11 m as defined by the ""International System of Units"", Bureau International des Poids et Mesures ",6,,1,1,1,1,INTEGER ,11,44,2076 +LOCATION_TOLERANCE_RADIAL ,Nadir Earth location tolerance radial ,,m ,1,1,1,1,INTEGER ,11,44,2120 +LOCATION_TOLERANCE_CROSSTRACK ,Nadir Earth location tolerance cross-track ,,m ,1,1,1,1,INTEGER ,11,44,2164 +LOCATION_TOLERANCE_ALONGTRACK ,Nadir Earth location tolerance along-track ,,m ,1,1,1,1,INTEGER ,11,44,2208 +YAW_ERROR ,Constant Yaw attitude error ,3,deg ,1,1,1,1,INTEGER ,11,44,2252 +ROLL_ERROR ,Constant Roll attitude error ,3,deg ,1,1,1,1,INTEGER ,11,44,2296 +PITCH_ERROR ,Constant Pitch attitude error ,3,deg ,1,1,1,1,INTEGER ,11,44,2340 +SUBSAT_LATITUDE_START ,Latitude of sub-satellite point at start of the data set ,3,deg ,1,1,1,1,INTEGER ,11,44,2384 +SUBSAT_LONGITUDE_START ,Longitude of sub-satellite point at start of the data set ,3,deg ,1,1,1,1,INTEGER ,11,44,2428 +SUBSAT_LATITUDE_END ,Latitude of sub-satellite point at end of the data set ,3,deg ,1,1,1,1,INTEGER ,11,44,2472 +SUBSAT_LONGITUDE_END ,Longitude of sub-satellite point at end of the data set ,3,deg ,1,1,1,1,INTEGER ,11,44,2516 +LEAP_SECOND ,"Occurence of Leap second within the product. Field is set to -1, 0 or +1 dependent upon occurrence of leap second and direction. ",,,1,1,1,1,INTEGER ,2,35,2560 +LEAP_SECOND_UTC ,"UTC time of occurrence of the Leap Second (If no leap second in the product, value is null) ",,,1,1,1,1,GENERAL TIME ,15,48,2595 +TOTAL_RECORDS ,Total count of all records in the product ,,,1,1,1,1,U-INTEGER ,6,39,2643 +TOTAL_MPHR ,Total count of all MPHRs in product (should always be 1!) ,,,1,1,1,1,U-INTEGER ,6,39,2682 +TOTAL_SPHR ,Total count of all SPHRs in product (should be 0 or 1 only) ,,,1,1,1,1,U-INTEGER ,6,39,2721 +TOTAL_IPR ,Total count of all IPRs in the product ,,,1,1,1,1,U-INTEGER ,6,39,2760 +TOTAL_GEADR ,Total count of all GEADRs in the product ,,,1,1,1,1,U-INTEGER ,6,39,2799 +TOTAL_GIADR ,Total count of all GIADRs in the product ,,,1,1,1,1,U-INTEGER ,6,39,2838 +TOTAL_VEADR ,Total count of all VEADRs in the product ,,,1,1,1,1,U-INTEGER ,6,39,2877 +TOTAL_VIADR ,Total count of all VIADRs in the product ,,,1,1,1,1,U-INTEGER ,6,39,2916 +TOTAL_MDR ,Total count of all MDRs in the product ,,,1,1,1,1,U-INTEGER ,6,39,2955 +COUNT_DEGRADED_INST_MDR ,Count of MDRs with degradation due to instrument problems ,,,1,1,1,1,U-INTEGER ,6,39,2994 +COUNT_DEGRADED_PROC_MDR ,Count of MDRs with degradation due to processing problems ,,,1,1,1,1,U-INTEGER ,6,39,3033 +COUNT_DEGRADED_INST_MDR_BLOCKS ,Count of the number of blocks of MDRs degraded due to degraded instrument ,,,1,1,1,1,U-INTEGER ,6,39,3072 +COUNT_DEGRADED_PROC_MDR_BLOCKS ,Count of the number of blocks of MDRs degraded due to degraded processing ,,,1,1,1,1,U-INTEGER ,6,39,3111 +DURATION_OF_PRODUCT ,The duration of the product in milliseconds ,,ms ,1,1,1,1,U-INTEGER ,8,41,3150 +MILLISECONDS_OF_DATA_PRESENT ,The total amount of data present in the product ,,ms ,1,1,1,1,U-INTEGER ,8,41,3191 +MILLISECONDS_OF_DATA_MISSING ,The total amount of data missing from the prodcut ,,ms ,1,1,1,1,U-INTEGER ,8,41,3232 +SUBSETTED_PRODUCT ,Set when product has been subsetted (e.g. geographically subsetted using a region of interest filter). Implies the presence of one or more UMARF GIADRs in GAD section for product retrieved from UMARF. ,,,1,1,1,1,BOOLEAN ,1,34,3273 +Size of the Record ,,,,,,,,,,,3307 diff --git a/satpy/readers/epsnative_reader.py b/satpy/readers/epsnative_reader.py index cb171e42fb..c0e7e401e8 100644 --- a/satpy/readers/epsnative_reader.py +++ b/satpy/readers/epsnative_reader.py @@ -22,6 +22,18 @@ # source including revision history and details on authorship can be found at # https://gitlab.eumetsat.int/open-source/data-tailor-plugins/epct_plugin_gis +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + """FIXME DOC. The ``epsnative_reader`` module provides the support to read from EPS native products records info diff --git a/satpy/tests/reader_tests/test_iasisnd02.py b/satpy/tests/reader_tests/test_iasisnd02.py new file mode 100644 index 0000000000..55c2b295dc --- /dev/null +++ b/satpy/tests/reader_tests/test_iasisnd02.py @@ -0,0 +1,131 @@ +# Copyright 2017-2022, European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) +# 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 . + +# This module is based on source code obtained from the +# epct_plugin_gis package developed by B-Open Solutions srl for EUMETSAT under +# contract EUM/C0/17/4600001943/0PN and released under Apache License +# Version 2.0, January 2004, http://www.apache.org/licenses/. The original +# source including revision history and details on authorship can be found at +# https://gitlab.eumetsat.int/open-source/data-tailor-plugins/epct_plugin_gis + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Test reading IASI L2 SND.""" + +import datetime +import shutil + +import pytest +import requests + +_url_sample_file = ("https://go.dwd-nextcloud.de/index.php/s/z87KfL72b9dM5xm/download/" + "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") + + +@pytest.fixture(scope="session") +def sample_file(tmp_path_factory): + """Obtain sample file.""" + fn = tmp_path_factory.mktemp("data") / "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat" + data = requests.get(_url_sample_file, stream=True) + with fn.open(mode="wb") as fp: + shutil.copyfileobj(data.raw, fp) + return fn + + +def test_read_giadr(sample_file): + """Test reading GIADR.""" + from satpy.readers.epsnative_reader import assemble_descriptor + from satpy.readers.iasi_l2_eps import read_giadr + descriptor = assemble_descriptor("IASISND02") + class_data = read_giadr(sample_file, descriptor) + + assert len(class_data) == 19 + assert class_data["PRESSURE_LEVELS_OZONE"]["units"] == "Pa" + assert class_data["PRESSURE_LEVELS_OZONE"]["values"].size == 101 + assert class_data["PRESSURE_LEVELS_OZONE"]["values"].min() == 50 + assert class_data["PRESSURE_LEVELS_OZONE"]["values"].max() == 11000000 + + +def test_datetime_to_second_since_2000(): + """Test counting seconds since 2000.""" + from satpy.readers.iasi_l2_eps import datetime_to_second_since_2000 + date = datetime.datetime(2000, 1, 1, 0, 1) + sec = datetime_to_second_since_2000(date) + + assert sec == 60 + + +def test_read_all_rows(sample_file): + """Test reading all rows from data.""" + from satpy.readers.epsnative_reader import assemble_descriptor + from satpy.readers.iasi_l2_eps import read_all_rows + sensing_start = datetime.datetime.strptime("20190605002352Z", "%Y%m%d%H%M%SZ") + sensing_stop = datetime.datetime.strptime("20190605020856Z", "%Y%m%d%H%M%SZ") + descriptor = assemble_descriptor("IASISND02") + mdr_class_offset = 271719118 # last row + with open(sample_file, "rb") as epsfile_obj: + data_before_errors_section, algorithms_data, errors_data = read_all_rows( + epsfile_obj, descriptor, mdr_class_offset, sensing_start, sensing_stop + ) + + assert len(data_before_errors_section) == 1 + assert data_before_errors_section[0]["INTEGRATED_CO"].min() == 7205 + + assert len(algorithms_data) == 1 + assert algorithms_data[0]["params"] == {"NERR": 69, "CO_NBR": 79, "HNO3_NBR": 0, "O3_NBR": 0} + + assert len(errors_data) == 1 + assert errors_data[0]["SURFACE_Z"].min() == 30 + + +def test_read_product_data(sample_file): + """Test reading product data.""" + from satpy.readers.epsnative_reader import assemble_descriptor + from satpy.readers.iasi_l2_eps import read_product_data + sensing_start = datetime.datetime.strptime("20190605002352Z", "%Y%m%d%H%M%SZ") + sensing_stop = datetime.datetime.strptime("20190605020856Z", "%Y%m%d%H%M%SZ") + descriptor = assemble_descriptor("IASISND02") + mdr_class_offset = 271279642 # last row + with open(sample_file, "rb") as epsfile_obj: + ( + stacked_data_before_errors, + stacked_algo_data, + stacked_errors_data, + ) = read_product_data( + epsfile_obj, descriptor, mdr_class_offset, sensing_start, sensing_stop + ) + + assert stacked_data_before_errors["INTEGRATED_CO"].shape == (2, 120) + assert list(stacked_algo_data.keys()) == ["CO"] + assert stacked_errors_data["SURFACE_Z"].shape == (2, 120) + + +def test_load(sample_file, tmp_path): + """Test read at a scene and write again.""" + from satpy import Scene + sc = Scene(filenames=[sample_file], reader=["iasi_l2_eps"]) + sc.load(["atmospheric_temperature"]) From 15d5fdd5689abe5794c9d011a6d4e4efb64f051b Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Mon, 10 Jul 2023 14:22:53 +0200 Subject: [PATCH 11/23] Try to combine the approaches --- satpy/etc/readers/iasi_l2_eps.yaml | 2 +- satpy/readers/eps_base.py | 6 +- satpy/readers/eps_l2.py | 219 +++++++++++++++++++++ satpy/tests/reader_tests/test_iasisnd02.py | 11 +- 4 files changed, 232 insertions(+), 6 deletions(-) create mode 100644 satpy/readers/eps_l2.py diff --git a/satpy/etc/readers/iasi_l2_eps.yaml b/satpy/etc/readers/iasi_l2_eps.yaml index 160e060c70..315770a36f 100644 --- a/satpy/etc/readers/iasi_l2_eps.yaml +++ b/satpy/etc/readers/iasi_l2_eps.yaml @@ -15,7 +15,7 @@ reader: file_types: iasi_l2_eps: - file_reader: !!python/name:satpy.readers.iasi_l2_eps.EPSIASIL2FileHandler + file_reader: !!python/name:satpy.readers.eps_l2.EPSIASIL2FileHandler # see https://www.eumetsat.int/media/40048 §7 for file pattern file_patterns: - IASI_SND_02_{platform_short_name}_{start_time:%Y%m%d%H%M%SZ}_{end_time:%Y%m%d%H%M%SZ}_{processing_mode}_{disposition_mode}_{creation_time:%Y%m%d%H%M%SZ} diff --git a/satpy/readers/eps_base.py b/satpy/readers/eps_base.py index d33e13ebeb..6ffcaaefdb 100644 --- a/satpy/readers/eps_base.py +++ b/satpy/readers/eps_base.py @@ -118,10 +118,12 @@ def read_records(filename, format_definition): return sections, form -def create_xarray(arr): +def create_xarray(arr, dims=("y", "x"), attrs=None): """Create xarray with correct dimensions.""" res = arr - res = xr.DataArray(res, dims=['y', 'x']) + if attrs is None: + attrs = {} + res = xr.DataArray(res, dims=dims, attrs=attrs) return res diff --git a/satpy/readers/eps_l2.py b/satpy/readers/eps_l2.py new file mode 100644 index 0000000000..7a61c9ecb2 --- /dev/null +++ b/satpy/readers/eps_l2.py @@ -0,0 +1,219 @@ +# 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 . + +"""Reader for EPS level 2 data. Uses xml files as a format description.""" + +import logging + +import numpy as np + +from satpy.readers.eps_base import EPSBaseFileHandler, create_xarray + +from . import epsnative_reader + +logger = logging.getLogger(__name__) + +var_2nd_dim = { + "co_cp_air": "co_nbr", + "co_cp_co_a": "co_nbr", + "co_x_co": "co_nbr", + "co_h_eigenvalues": "co_nbr", + "co_h_eigenvectors": "co_nbr", + "ozone_error": "nerr", + "temperature_error": "nerr", + "water_vapour_error": "nerr", +} +var_3rd_dim = { + "fg_atmospheric_temperature": "nlt", + "fg_atmospheric_water_vapour": "nlq", + "fg_atmospheric_ozone": "nlo", + "atmospheric_temperature": "nlt", + "atmospheric_water_vapour": "nlq", + "atmospheric_ozone": "nlo", +} + +NR_FOV = 120 + + +class EPSIASIL2FileHandler(EPSBaseFileHandler): + """EPS level 2 reader for IASI data. + + Reader for the IASI Level 2 combined sounding products in native format. + + Overview of the data including links to the product user guide, product format + specification, validation reports, and other documents, can be found at the + EUMETSAT Data Services at https://data.eumetsat.int/product/EO:EUM:DAT:METOP:IASSND02 + + """ + + xml_conf = "eps_iasil2_9.0.xml" + mdr_subclass = 1 + + def __init__(self, filename, filename_info, filetype_info): + """Initialise handler.""" + super().__init__(filename, filename_info, filetype_info) + self.descr = epsnative_reader.assemble_descriptor("IASISND02") + + def get_dataset(self, key, info): + """Get calibrated data.""" + dsname = key["name"].upper() + vals = self[dsname] + dims = get_var_dimensions(dsname, vals, self.dimensions) + # FIXME: get variable attributes + return create_xarray(vals, dims=list(dims.keys())) + + def available_datasets(self, configured_datasets=None): + """Get available datasets.""" + if self.sections is None: + self._read_all() + common = {"file_type": "iasi_l2_eps", "resolution": 12000} + for key in self.keys(): + yield (True, {"name": key.lower()} | common) + + def _read_all(self): + """Read all variables, lazily.""" + super()._read_all() + self.giadr = read_giadr(self.filename, self.descr) + self.dimensions = assemble_dimensions( + self.giadr, + self.scanlines, + self["NUM_ERROR_DATA"].max()) # Correct? + + +def read_giadr(input_product, descriptor, **kwargs): + """FIXME DOC. + + :param input_product: + :param descriptor: + :param kwargs: + :return: + """ + grh_size = 20 # bytes + ipr_sequence = epsnative_reader.read_ipr_sequence(input_product) + giadr_class = [cl for cl in ipr_sequence if "giadr" in cl["class"]][0] + with open(input_product, "rb") as eps_fileobj: + eps_fileobj.seek(giadr_class["offset"] + grh_size) + class_data = {} + giadr_descriptor = descriptor.get(("giadr", 1, 4)) + for _, row in giadr_descriptor.iterrows(): + # empty row or general header + if np.isnan(row["FIELD_SIZE"]) or row["FIELD"] == "RECORD_HEADER": + continue + else: + class_data[row["FIELD"].strip()] = { + "description": row["DESCRIPTION"], + "units": row["UNITS"], + "values": read_values(eps_fileobj, row, **kwargs), + } + return class_data + + +def read_values(eps_fileobj, row, reshape=False): + """FIXME DOC. + + :param _io.BinaryIO eps_fileobj: + :param pandas.core.series.Series row: + :param bool reshape: + :return numpy.ndarray: + """ + dtype = epsnative_reader.reckon_dtype(row["TYPE"].strip()) + shape = [int(row[k]) for k in row.keys() if k.lower().startswith("dim")] + shape = [k for k in shape if k > 1] + count = int(np.product(shape)) + values = np.fromfile(eps_fileobj, dtype, count) + # if np.isfinite(row['SF']): + # values = values / 10 ** row['SF'] + if row["TYPE"].startswith("v-"): + values = values["f1"] * 10.0 ** -values["f0"] + if row["TYPE"].strip().lower() == "short_cds_time": + values = [epsnative_reader.scds_to_datetime(*v) for v in values] + values = values[0] if len(values) == 1 else values + if reshape: + try: + values = values.reshape(tuple(reversed(shape))).T + except AttributeError: + values = values + return values + + +def assemble_dimensions(giadr, nr_rows, max_nerr): + """See Parameter Table (8.6) in EUM/OPS-EPS/MAN/04/0033. + + :param dict giadr: + :param int nr_rows: + :param int max_nerr: + :return dict: + """ + dimensions = { + "across_track": NR_FOV, + "along_track": nr_rows, + "cloud_formations": 3, + "co_nbr": NR_FOV, # pad to FOV number; real value for CO_NBR is provided in MDR + "nerr": max_nerr, # max value of NERR values provided in MDR + "nerro": giadr["NUM_OZONE_PCS"]["values"][0] + * (giadr["NUM_OZONE_PCS"]["values"][0] + 1) + / 2, + "nerrt": giadr["NUM_TEMPERATURE_PCS"]["values"][0] + * (giadr["NUM_TEMPERATURE_PCS"]["values"][0] + 1) + / 2, + "nerrw": giadr["NUM_WATER_VAPOUR_PCS"]["values"][0] + * (giadr["NUM_WATER_VAPOUR_PCS"]["values"][0] + 1) + / 2, + "neva_co": round(giadr["FORLI_NUM_LAYERS_CO"]["values"][0] / 2), + "neve_co": round(giadr["FORLI_NUM_LAYERS_CO"]["values"][0] / 2) + * giadr["FORLI_NUM_LAYERS_CO"]["values"][0], + "new": giadr["NUM_SURFACE_EMISSIVITY_WAVELENGTHS"]["values"][0], + "nl_co": giadr["FORLI_NUM_LAYERS_CO"]["values"][0], + "nl_hno3": giadr["FORLI_NUM_LAYERS_HNO3"]["values"][0], + "nl_o3": giadr["FORLI_NUM_LAYERS_O3"]["values"][0], + "nl_so2": giadr["BRESCIA_NUM_ALTITUDES_SO2"]["values"][0], + "nlo": giadr["NUM_PRESSURE_LEVELS_OZONE"]["values"][0], + "nlq": giadr["NUM_PRESSURE_LEVELS_HUMIDITY"]["values"][0], + "nlt": giadr["NUM_PRESSURE_LEVELS_TEMP"]["values"][0], + "npco": giadr["NUM_OZONE_PCS"]["values"][0], + "npct": giadr["NUM_TEMPERATURE_PCS"]["values"][0], + "npcw": giadr["NUM_WATER_VAPOUR_PCS"]["values"][0], + } + return dimensions + + +def get_var_dimensions(var_name, values, dimensions): + """FIXME DOC. + + :param str var_name: + :param numpy.ndarray values: + :param dict dimensions: + :return dict: + """ + # nlt, nlq, nlo are all equals to 101, so the automatic detection fails + dims_keys = list(dimensions.keys()) + sizes = list(dimensions.values()) + dims = {} + for idx, k in enumerate(values.shape): + if idx == 0: + dims["y"] = dimensions["along_track"] + elif idx == 1: + if var_name in var_2nd_dim: + dims[var_2nd_dim[var_name]] = k + else: + dims["x"] = k + else: + if var_name in var_3rd_dim: + dims[var_3rd_dim[var_name]] = k + else: + dim_name = dims_keys[sizes.index(k)] + dims[dim_name] = k + return dims diff --git a/satpy/tests/reader_tests/test_iasisnd02.py b/satpy/tests/reader_tests/test_iasisnd02.py index 55c2b295dc..c4cb63278e 100644 --- a/satpy/tests/reader_tests/test_iasisnd02.py +++ b/satpy/tests/reader_tests/test_iasisnd02.py @@ -39,14 +39,17 @@ import datetime import shutil +import dask import pytest import requests +from ..utils import CustomScheduler + _url_sample_file = ("https://go.dwd-nextcloud.de/index.php/s/z87KfL72b9dM5xm/download/" "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") -@pytest.fixture(scope="session") +@pytest.fixture(scope="module") def sample_file(tmp_path_factory): """Obtain sample file.""" fn = tmp_path_factory.mktemp("data") / "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat" @@ -127,5 +130,7 @@ def test_read_product_data(sample_file): def test_load(sample_file, tmp_path): """Test read at a scene and write again.""" from satpy import Scene - sc = Scene(filenames=[sample_file], reader=["iasi_l2_eps"]) - sc.load(["atmospheric_temperature"]) + with dask.config.set(scheduler=CustomScheduler(max_computes=0)): + sc = Scene(filenames=[sample_file], reader=["iasi_l2_eps"]) + sc.load(["atmospheric_temperature"]) + assert isinstance(sc["atmospheric_temperature"], dask.array.Array) From f902144ee417b9b32dffe5329976adc2aafeda46 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Tue, 11 Jul 2023 10:51:49 +0200 Subject: [PATCH 12/23] Improved tests for IASI reading Prior to speeding up the reading, improve tests to test actual values to be read. --- satpy/readers/iasi_l2_eps.py | 11 ++-- satpy/tests/reader_tests/test_iasisnd02.py | 58 +++++++++++++++------- 2 files changed, 46 insertions(+), 23 deletions(-) diff --git a/satpy/readers/iasi_l2_eps.py b/satpy/readers/iasi_l2_eps.py index 2c7bf57d23..797c5bf216 100644 --- a/satpy/readers/iasi_l2_eps.py +++ b/satpy/readers/iasi_l2_eps.py @@ -61,8 +61,7 @@ def get_dataset(self, dataid, dataset_info): if self._nc is None: self._nc = self._get_netcdf_dataset() da = self._nc[dataid["name"]] - # scale factor is not relevant and gets in the way - da.attrs.pop("scale_factor", None) + da = da * da.attrs.pop("scale_factor", 1) return da def _get_netcdf_dataset(self): @@ -610,12 +609,12 @@ def get_var_dimensions(var_name, values, dimensions): dims = {} for idx, k in enumerate(values.shape): if idx == 0: - dims["along_track"] = dimensions["along_track"] + dims["y"] = dimensions["y"] elif idx == 1: if var_name in var_2nd_dim: dims[var_2nd_dim[var_name]] = k else: - dims["across_track"] = k + dims["x"] = k else: if var_name in var_3rd_dim: dims[var_3rd_dim[var_name]] = k @@ -670,8 +669,8 @@ def assemble_dimensions(giadr, nr_rows, max_nerr): :return dict: """ dimensions = { - "across_track": NR_FOV, - "along_track": nr_rows, + "x": NR_FOV, + "y": nr_rows, "cloud_formations": 3, "co_nbr": NR_FOV, # pad to FOV number; real value for CO_NBR is provided in MDR "nerr": max_nerr, # max value of NERR values provided in MDR diff --git a/satpy/tests/reader_tests/test_iasisnd02.py b/satpy/tests/reader_tests/test_iasisnd02.py index c4cb63278e..ad38188d77 100644 --- a/satpy/tests/reader_tests/test_iasisnd02.py +++ b/satpy/tests/reader_tests/test_iasisnd02.py @@ -37,9 +37,11 @@ """Test reading IASI L2 SND.""" import datetime +import pathlib import shutil import dask +import numpy as np import pytest import requests @@ -52,6 +54,9 @@ @pytest.fixture(scope="module") def sample_file(tmp_path_factory): """Obtain sample file.""" + fn = pathlib.Path("/media/nas/x21308/IASI/IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") + if fn.exists(): + return fn fn = tmp_path_factory.mktemp("data") / "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat" data = requests.get(_url_sample_file, stream=True) with fn.open(mode="wb") as fp: @@ -86,9 +91,11 @@ def test_read_all_rows(sample_file): """Test reading all rows from data.""" from satpy.readers.epsnative_reader import assemble_descriptor from satpy.readers.iasi_l2_eps import read_all_rows + sensing_start = datetime.datetime.strptime("20190605002352Z", "%Y%m%d%H%M%SZ") sensing_stop = datetime.datetime.strptime("20190605020856Z", "%Y%m%d%H%M%SZ") descriptor = assemble_descriptor("IASISND02") + mdr_class_offset = 271719118 # last row with open(sample_file, "rb") as epsfile_obj: data_before_errors_section, algorithms_data, errors_data = read_all_rows( @@ -109,28 +116,45 @@ def test_read_product_data(sample_file): """Test reading product data.""" from satpy.readers.epsnative_reader import assemble_descriptor from satpy.readers.iasi_l2_eps import read_product_data - sensing_start = datetime.datetime.strptime("20190605002352Z", "%Y%m%d%H%M%SZ") - sensing_stop = datetime.datetime.strptime("20190605020856Z", "%Y%m%d%H%M%SZ") descriptor = assemble_descriptor("IASISND02") - mdr_class_offset = 271279642 # last row - with open(sample_file, "rb") as epsfile_obj: - ( - stacked_data_before_errors, - stacked_algo_data, - stacked_errors_data, - ) = read_product_data( - epsfile_obj, descriptor, mdr_class_offset, sensing_start, sensing_stop - ) + mdr_class_offset = 260002007 # not quite the last row - assert stacked_data_before_errors["INTEGRATED_CO"].shape == (2, 120) - assert list(stacked_algo_data.keys()) == ["CO"] - assert stacked_errors_data["SURFACE_Z"].shape == (2, 120) + with open(sample_file, "rb") as epsfile_obj: + (sdbs, sad, sed) = read_product_data( + epsfile_obj, + descriptor, + mdr_class_offset, + datetime.datetime(2019, 6, 5, 0, 23, 52), + datetime.datetime(2019, 6, 5, 2, 8, 56)) + assert sdbs["SURFACE_TEMPERATURE"].shape == (37, 120) + assert sdbs["ATMOSPHERIC_TEMPERATURE"].shape == (37, 101, 120) + np.testing.assert_array_equal( + sdbs["SURFACE_TEMPERATURE"][0, 58:64], + np.array([27252, 27223, 27188, 27136, 27227, 65535], dtype=np.uint16)) + np.testing.assert_array_equal( + sdbs["SURFACE_TEMPERATURE"][30, 80:86], + np.array([27366, 65535, 27368, 27369, 27350, 27406], dtype=np.uint16)) + assert isinstance(sdbs["ATMOSPHERIC_TEMPERATURE"], dask.array.Array) def test_load(sample_file, tmp_path): """Test read at a scene and write again.""" from satpy import Scene with dask.config.set(scheduler=CustomScheduler(max_computes=0)): - sc = Scene(filenames=[sample_file], reader=["iasi_l2_eps"]) - sc.load(["atmospheric_temperature"]) - assert isinstance(sc["atmospheric_temperature"], dask.array.Array) + sc = Scene(filenames=[sample_file], reader=["iasi_l2_eps2"]) + sc.load(["surface_temperature"]) + assert sc["surface_temperature"].dims == ("y", "x") + np.testing.assert_allclose( + sc["surface_temperature"][0, 104:110], + np.array([262.55, 261.83, 265.85, 265.81, 256.15, 255.23])) + np.testing.assert_allclose( + sc["surface_temperature"][1, 104:110], + np.array([261.9, 266.88, 268.47, 266.76, 256.63, 261.32])) + np.testing.assert_allclose( + sc["surface_temperature"][2, 108:112], + np.array([262.33, 265.42, 268.32, 266.21])) + np.testing.assert_allclose( + sc["surface_temperature"][10, 83:89], + np.array([271.11, 280.07, 278.65, 278.84, 278.37, 281.54])) + + assert isinstance(sc["surface_temperature"], dask.array.Array) From 8eb657cb888cd5a0cd2022baf5fb0e5ea1168306 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Tue, 11 Jul 2023 17:44:49 +0200 Subject: [PATCH 13/23] Copy EPCT test file for epsnative reader --- .../reader_tests/test_epsnative_reader.py | 195 ++++++++++++++++++ 1 file changed, 195 insertions(+) create mode 100644 satpy/tests/reader_tests/test_epsnative_reader.py diff --git a/satpy/tests/reader_tests/test_epsnative_reader.py b/satpy/tests/reader_tests/test_epsnative_reader.py new file mode 100644 index 0000000000..5eeb8091ed --- /dev/null +++ b/satpy/tests/reader_tests/test_epsnative_reader.py @@ -0,0 +1,195 @@ +# Copyright 2017-2022, European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# -*- coding: utf-8 -*- +import os + +from epct_plugin_gis import epsnative_reader + + +TEST_DATA_PATH = os.environ.get("EPCT_TEST_DATA_DIR", "") +TEST_DATA = os.path.join( + TEST_DATA_PATH, + "EPS", + "AVHRRL1", + "AVHR_xxx_1B_M02_20181212113403Z_20181212131303Z_N_O_20181212130905Z.nat", +) +PDUS_DATA = [ + os.path.join( + TEST_DATA_PATH, + "EPS", + "AVHRRL1", + "AVHR_xxx_1B_M01_20180120003103Z_20180120003403Z_N_O_20180120004248Z", + ), + os.path.join( + TEST_DATA_PATH, + "EPS", + "AVHRRL1", + "AVHR_xxx_1B_M01_20180120003403Z_20180120003703Z_N_O_20180120013050Z", + ), + os.path.join( + TEST_DATA_PATH, + "EPS", + "AVHRRL1", + "AVHR_xxx_1B_M01_20180120003703Z_20180120004003Z_N_O_20180120013202Z", + ), +] + +SAMPLE_PATH = os.path.join(os.path.dirname(__file__), "sample_data") +SAMPLE_DATA = os.path.join(SAMPLE_PATH, "AVHRRL1.nat") + + +def test_get_class_tuple(): + class_string = "class_1_2_3" + assert epsnative_reader.get_class_tuple(class_string) == ("class", 1, 2, 3) + + class_string = "class_dummy_dummy2" + assert epsnative_reader.get_class_tuple(class_string) == ("class", "dummy_dummy2") + + +def test_csv_list_by_product(): + csv_list = epsnative_reader.csv_list_by_product("AVHRRL1") + exp_csv = ["giadr_1_3.csv", "giadr_2_2.csv", "mdr_2_5.csv", "mphr_0_2.csv"] + + assert len(csv_list) == 4 + for csv_path in csv_list: + assert os.path.isfile(csv_path) + assert os.path.basename(csv_path) in exp_csv + + +def test_assemble_descriptor(): + descriptor = epsnative_reader.assemble_descriptor("AVHRRL1") + exp_keys = {("giadr", 1, 3), ("giadr", 2, 2), ("mdr", 2, 5), ("mphr", 0, 2)} + + assert set(descriptor.keys()) == exp_keys + + +def test_grh_reader(): + grh = epsnative_reader.grh_reader(SAMPLE_DATA) + + assert len(grh) == 6 + assert grh[0] == "mphr" + assert grh[1] == 0 + assert grh[2] == 2 + assert grh[3] == 3307 + assert grh[4].isoformat() == "2018-01-22T13:28:03.130000" + assert grh[5].isoformat() == "2018-01-22T13:28:04.797000" + + +def test_find_mphr_csv(): + mphr_csv = epsnative_reader.find_mphr_csv() + + assert os.path.isfile(mphr_csv) + + +def test_mphr_reader(): + mphr_content = epsnative_reader.mphr_reader(SAMPLE_DATA) + + assert len(mphr_content) == 72 + assert mphr_content["INSTRUMENT_ID"] == "AVHR" + assert mphr_content["ORBIT_START"] == 58431 + assert mphr_content["ORBIT_END"] == 58432 + + +def test_first_class_occurrence(): + class_name = "mphr" + grh, offset = epsnative_reader.first_class_occurrence(SAMPLE_DATA, class_name) + assert offset == 0 + + class_name = "giadr" + grh, offset = epsnative_reader.first_class_occurrence(SAMPLE_DATA, class_name) + assert offset == 4587 + assert grh[0] == class_name + assert grh[1] == 1 + assert grh[2] == 3 + + class_name = "mdr" + grh, offset = epsnative_reader.first_class_occurrence(SAMPLE_DATA, class_name) + assert offset == 5077 + assert grh[0] == class_name + assert grh[1] == 2 + assert grh[2] == 4 + + +def test_read_ipr_sequence(): + ipr_sequence = epsnative_reader.read_ipr_sequence(SAMPLE_DATA) + + assert len(ipr_sequence) == 11 + assert ipr_sequence[-1]["class"] == ("mdr", 2) + assert ipr_sequence[-1]["offset"] == 5077 + + +def test_read_grh_of_target_class(): + target_offset = 0 + with open(SAMPLE_DATA, "rb") as eps_fileobj: + target_grh = epsnative_reader.read_grh_of_target_class(eps_fileobj, target_offset) + assert target_grh[0] == "mphr" + assert target_grh[1] == 0 + assert target_grh[2] == 2 + assert target_grh[3] == 3307 + + target_offset = 5077 + with open(SAMPLE_DATA, "rb") as eps_fileobj: + target_grh = epsnative_reader.read_grh_of_target_class(eps_fileobj, target_offset) + assert target_grh[0] == "mdr" + assert target_grh[1] == 2 + assert target_grh[2] == 4 + assert target_grh[3] == 26660 + + +def test_add_info_about_class(): + mphr = epsnative_reader.mphr_reader(SAMPLE_DATA) + current_class = {"class": "", "offset": 5077} + with open(SAMPLE_DATA, "rb") as eps_fileobj: + current_class = epsnative_reader.add_info_about_class(eps_fileobj, current_class, {}, mphr) + + assert current_class["class_id"] == ("mdr", 2, 4) + assert current_class["class_size"] == 26660 + assert current_class["nr_records"] == 10.0 + + +def test_reckon_dtype(): + assert epsnative_reader.reckon_dtype("boolean") == ">i1" + assert epsnative_reader.reckon_dtype("u-integer1") == ">u1" + assert epsnative_reader.reckon_dtype("vinteger2") == "byte,>i2" + assert epsnative_reader.reckon_dtype("bitst(16)") == ">u2" + assert epsnative_reader.reckon_dtype("bitst(24)") == ">u1,>u1,>u1" + + +def test_bands_to_records_reader(): + records = epsnative_reader.bands_to_records_reader("AVHRRL1") + + assert len(records) == 7 + assert records["channel_5"]["record_name"] == "SCENE_RADIANCES" + assert records["channel_5"]["band_position"] == 4 + assert records["channel_5"]["convention"] == "BIL" + + +def test_add_record_info_to_band(): + record_info = epsnative_reader.add_record_info_to_band("AVHRRL1") + + assert len(record_info) == 7 + assert record_info[0]["band_id"] == "channel_1" + assert record_info[0]["shape"] == [2048, 5] + assert record_info[-1]["band_id"] == "longitude" + assert record_info[-1]["shape"] == [2, 103] + + +def test_create_toc(): + toc = epsnative_reader.create_toc(SAMPLE_DATA) + + assert len(toc) == 11 + assert toc[0]["class_id"] == ("geadr", 1, 1) + assert toc[-1]["class_id"] == ("mdr", 2, 4) + assert toc[-1]["y_offset"] == 0 From 74d0081d986e56c05d3bd8ffde489a198268169e Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Tue, 11 Jul 2023 18:05:17 +0200 Subject: [PATCH 14/23] Support memmap in grh_reader When reading grh, support interpreting a memmap. --- satpy/readers/epsnative_reader.py | 6 +- satpy/tests/conftest.py | 21 +++++++ .../reader_tests/test_epsnative_reader.py | 60 +++++++++++++++++-- satpy/tests/reader_tests/test_iasisnd02.py | 20 ------- 4 files changed, 80 insertions(+), 27 deletions(-) diff --git a/satpy/readers/epsnative_reader.py b/satpy/readers/epsnative_reader.py index c0e7e401e8..fd53d67e62 100644 --- a/satpy/readers/epsnative_reader.py +++ b/satpy/readers/epsnative_reader.py @@ -154,7 +154,11 @@ def grh_reader(eps_fileobj): :param eps_fileobj: :return: """ - grh_array = np.fromfile(eps_fileobj, np.dtype(grh_type), 1) + dtp = np.dtype(grh_type) + if isinstance(eps_fileobj, np.memmap): + grh_array = eps_fileobj[:dtp.itemsize].view(dtp) + else: + grh_array = np.fromfile(eps_fileobj, dtp, 1) # When the EOF is reached if grh_array.size == 0: return () diff --git a/satpy/tests/conftest.py b/satpy/tests/conftest.py index 8bcbea2093..501fb49530 100644 --- a/satpy/tests/conftest.py +++ b/satpy/tests/conftest.py @@ -21,8 +21,11 @@ """ import os +import pathlib +import shutil import pytest +import requests import satpy @@ -55,3 +58,21 @@ def include_test_etc(): """Tell Satpy to use the config 'etc' directory from the tests directory.""" with satpy.config.set(config_path=[TEST_ETC_DIR]): yield TEST_ETC_DIR + + +_url_sample_file = ( + "https://go.dwd-nextcloud.de/index.php/s/z87KfL72b9dM5xm/download/" + "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") + + +@pytest.fixture(scope="session") +def sample_file(tmp_path_factory): + """Obtain sample file.""" + fn = pathlib.Path("/media/nas/x21308/IASI/IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") + if fn.exists(): + return fn + fn = tmp_path_factory.mktemp("data") / "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat" + data = requests.get(_url_sample_file, stream=True) + with fn.open(mode="wb") as fp: + shutil.copyfileobj(data.raw, fp) + return fn diff --git a/satpy/tests/reader_tests/test_epsnative_reader.py b/satpy/tests/reader_tests/test_epsnative_reader.py index 5eeb8091ed..6acf9b57a0 100644 --- a/satpy/tests/reader_tests/test_epsnative_reader.py +++ b/satpy/tests/reader_tests/test_epsnative_reader.py @@ -1,4 +1,27 @@ # Copyright 2017-2022, European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) +# 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 . + +# This module is based on source code obtained from the +# epct_plugin_gis package developed by B-Open Solutions srl for EUMETSAT under +# contract EUM/C0/17/4600001943/0PN and released under Apache License +# Version 2.0, January 2004, http://www.apache.org/licenses/. The original +# source including revision history and details on authorship can be found at +# https://gitlab.eumetsat.int/open-source/data-tailor-plugins/epct_plugin_gis + # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -12,11 +35,15 @@ # See the License for the specific language governing permissions and # limitations under the License. -# -*- coding: utf-8 -*- +"""FIXME DOC.""" + +import datetime import os -from epct_plugin_gis import epsnative_reader +import numpy as np +import pytest +from satpy.readers import epsnative_reader TEST_DATA_PATH = os.environ.get("EPCT_TEST_DATA_DIR", "") TEST_DATA = os.path.join( @@ -51,6 +78,7 @@ def test_get_class_tuple(): + """FIXME DOC.""" class_string = "class_1_2_3" assert epsnative_reader.get_class_tuple(class_string) == ("class", 1, 2, 3) @@ -59,6 +87,7 @@ def test_get_class_tuple(): def test_csv_list_by_product(): + """FIXME DOC.""" csv_list = epsnative_reader.csv_list_by_product("AVHRRL1") exp_csv = ["giadr_1_3.csv", "giadr_2_2.csv", "mdr_2_5.csv", "mphr_0_2.csv"] @@ -69,31 +98,42 @@ def test_csv_list_by_product(): def test_assemble_descriptor(): + """FIXME DOC.""" descriptor = epsnative_reader.assemble_descriptor("AVHRRL1") exp_keys = {("giadr", 1, 3), ("giadr", 2, 2), ("mdr", 2, 5), ("mphr", 0, 2)} assert set(descriptor.keys()) == exp_keys -def test_grh_reader(): - grh = epsnative_reader.grh_reader(SAMPLE_DATA) +@pytest.mark.parametrize("how", ["string", "file", "mmap"]) +def test_grh_reader(sample_file, how): + """FIXME DOC.""" + if how == "file": + sample = open(sample_file, mode="rb") + elif how == "mmap": + sample = np.memmap(sample_file, mode="r", offset=0) + else: + sample = sample_file + grh = epsnative_reader.grh_reader(sample) assert len(grh) == 6 assert grh[0] == "mphr" assert grh[1] == 0 assert grh[2] == 2 assert grh[3] == 3307 - assert grh[4].isoformat() == "2018-01-22T13:28:03.130000" - assert grh[5].isoformat() == "2018-01-22T13:28:04.797000" + assert grh[4] == datetime.datetime(2019, 6, 5, 0, 23, 52, 653000) + assert grh[5] == datetime.datetime(2019, 6, 5, 2, 8, 56, 225000) def test_find_mphr_csv(): + """FIXME DOC.""" mphr_csv = epsnative_reader.find_mphr_csv() assert os.path.isfile(mphr_csv) def test_mphr_reader(): + """FIXME DOC.""" mphr_content = epsnative_reader.mphr_reader(SAMPLE_DATA) assert len(mphr_content) == 72 @@ -103,6 +143,7 @@ def test_mphr_reader(): def test_first_class_occurrence(): + """FIXME DOC.""" class_name = "mphr" grh, offset = epsnative_reader.first_class_occurrence(SAMPLE_DATA, class_name) assert offset == 0 @@ -123,6 +164,7 @@ def test_first_class_occurrence(): def test_read_ipr_sequence(): + """FIXME DOC.""" ipr_sequence = epsnative_reader.read_ipr_sequence(SAMPLE_DATA) assert len(ipr_sequence) == 11 @@ -131,6 +173,7 @@ def test_read_ipr_sequence(): def test_read_grh_of_target_class(): + """FIXME DOC.""" target_offset = 0 with open(SAMPLE_DATA, "rb") as eps_fileobj: target_grh = epsnative_reader.read_grh_of_target_class(eps_fileobj, target_offset) @@ -149,6 +192,7 @@ def test_read_grh_of_target_class(): def test_add_info_about_class(): + """FIXME DOC.""" mphr = epsnative_reader.mphr_reader(SAMPLE_DATA) current_class = {"class": "", "offset": 5077} with open(SAMPLE_DATA, "rb") as eps_fileobj: @@ -160,6 +204,7 @@ def test_add_info_about_class(): def test_reckon_dtype(): + """FIXME DOC.""" assert epsnative_reader.reckon_dtype("boolean") == ">i1" assert epsnative_reader.reckon_dtype("u-integer1") == ">u1" assert epsnative_reader.reckon_dtype("vinteger2") == "byte,>i2" @@ -168,6 +213,7 @@ def test_reckon_dtype(): def test_bands_to_records_reader(): + """FIXME DOC.""" records = epsnative_reader.bands_to_records_reader("AVHRRL1") assert len(records) == 7 @@ -177,6 +223,7 @@ def test_bands_to_records_reader(): def test_add_record_info_to_band(): + """FIXME DOC.""" record_info = epsnative_reader.add_record_info_to_band("AVHRRL1") assert len(record_info) == 7 @@ -187,6 +234,7 @@ def test_add_record_info_to_band(): def test_create_toc(): + """FIXME DOC.""" toc = epsnative_reader.create_toc(SAMPLE_DATA) assert len(toc) == 11 diff --git a/satpy/tests/reader_tests/test_iasisnd02.py b/satpy/tests/reader_tests/test_iasisnd02.py index ad38188d77..d4ebf9205c 100644 --- a/satpy/tests/reader_tests/test_iasisnd02.py +++ b/satpy/tests/reader_tests/test_iasisnd02.py @@ -37,32 +37,12 @@ """Test reading IASI L2 SND.""" import datetime -import pathlib -import shutil import dask import numpy as np -import pytest -import requests from ..utils import CustomScheduler -_url_sample_file = ("https://go.dwd-nextcloud.de/index.php/s/z87KfL72b9dM5xm/download/" - "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") - - -@pytest.fixture(scope="module") -def sample_file(tmp_path_factory): - """Obtain sample file.""" - fn = pathlib.Path("/media/nas/x21308/IASI/IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") - if fn.exists(): - return fn - fn = tmp_path_factory.mktemp("data") / "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat" - data = requests.get(_url_sample_file, stream=True) - with fn.open(mode="wb") as fp: - shutil.copyfileobj(data.raw, fp) - return fn - def test_read_giadr(sample_file): """Test reading GIADR.""" From 7f557c4342fb9600765c401cc2d86fc5e3775726 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Wed, 12 Jul 2023 15:25:25 +0200 Subject: [PATCH 15/23] Use dask and mmap for IASI L2 SND reading --- satpy/readers/epsnative_reader.py | 4 +- satpy/readers/iasi_l2_eps.py | 63 +++++---- satpy/tests/conftest.py | 19 ++- satpy/tests/reader_tests/conftest.py | 28 ++++ .../reader_tests/test_epsnative_reader.py | 17 +-- satpy/tests/reader_tests/test_iasisnd02.py | 130 ++++++++++++++---- 6 files changed, 191 insertions(+), 70 deletions(-) diff --git a/satpy/readers/epsnative_reader.py b/satpy/readers/epsnative_reader.py index fd53d67e62..095147d255 100644 --- a/satpy/readers/epsnative_reader.py +++ b/satpy/readers/epsnative_reader.py @@ -155,7 +155,9 @@ def grh_reader(eps_fileobj): :return: """ dtp = np.dtype(grh_type) - if isinstance(eps_fileobj, np.memmap): + # test for ndarray instead of memmap, because an empty slice from a memmap + # becomes an ndarray, and this happens at EOF + if isinstance(eps_fileobj, np.ndarray): grh_array = eps_fileobj[:dtp.itemsize].view(dtp) else: grh_array = np.fromfile(eps_fileobj, dtp, 1) diff --git a/satpy/readers/iasi_l2_eps.py b/satpy/readers/iasi_l2_eps.py index 797c5bf216..e21abec2d5 100644 --- a/satpy/readers/iasi_l2_eps.py +++ b/satpy/readers/iasi_l2_eps.py @@ -29,6 +29,7 @@ import datetime import itertools +import dask.array as da import numpy as np import xarray as xr @@ -132,7 +133,10 @@ def read_values(eps_fileobj, row, reshape=False): shape = [int(row[k]) for k in row.keys() if k.lower().startswith("dim")] shape = [k for k in shape if k > 1] count = int(np.product(shape)) - values = np.fromfile(eps_fileobj, dtype, count) + if isinstance(eps_fileobj, np.memmap): + values = eps_fileobj[:(np.dtype(dtype).itemsize*count)].view(dtype) + else: + values = np.fromfile(eps_fileobj, dtype, count) # if np.isfinite(row['SF']): # values = values / 10 ** row['SF'] if row["TYPE"].startswith("v-"): @@ -189,7 +193,8 @@ def set_values_in_mdr_descriptor(epsfile_obj, mdr_descriptor, mdr_class_offset, row_idx = row.index[0] row = row.squeeze() # read values - epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) +# epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) + epsfile_obj = epsfile_obj[mdr_class_offset + int(row["OFFSET"]):] value = read_values(epsfile_obj, row).astype(int)[0] # set the read values in MDR-descriptor mdr_descriptor.loc[mdr_descriptor["DIM2"] == field_name.lower(), "DIM2"] = value @@ -259,8 +264,10 @@ def read_errors(epsfile_obj, mdr_descriptor, mdr_class_offset, max_nerr): errors = {} for field in fields_to_read: row = mdr_descriptor.loc[mdr_descriptor["FIELD"] == field].squeeze() - epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) - values = read_values(epsfile_obj, row, reshape=True) + epsfile_mmap_subset = epsfile_obj[mdr_class_offset + int(row["OFFSET"]):] +# epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) +# values = read_values(epsfile_obj, row, reshape=True) + values = read_values(epsfile_mmap_subset, row, reshape=True) if field != "SURFACE_Z": if row["DIM2"] == 0: values = np.full((row["DIM1"], max_nerr), missing_value) @@ -303,8 +310,10 @@ def read_algorithm_sections( section_rows = mdr_descriptor.loc[mdr_descriptor["FIELD"].str.startswith(prefix)] section = {} for _, row in section_rows.iterrows(): - epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) - values = read_values(epsfile_obj, row, reshape=True) + epsfile_mmap_subset = epsfile_obj[mdr_class_offset + int(row["OFFSET"]):] +# epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) + values = read_values(epsfile_mmap_subset, row, reshape=True) +# values = read_values(epsfile_obj, row, reshape=True) if row["FIELD"] in vars_must_be_extend_to_fov: values = values.reshape(-1, 1) if values.ndim == 1 else values values = np.pad( @@ -429,11 +438,14 @@ def stack_non_algorithm_data(data): for key in stacked_data: stacked_data[key].append(row[key]) for key, content in stacked_data.items(): - stacked_data[key] = np.stack(content, axis=0) + if getattr(content[0], "size", 0) > 10: + stacked_data[key] = da.stack(content, axis=0) + else: + stacked_data[key] = np.stack(content, axis=0) return stacked_data -def read_records_before_error_section(epsfile_obj, mdr_descriptor, mdr_class_offset): +def read_records_before_error_section(epsfile_mmap, mdr_descriptor, mdr_class_offset): """FIXME DOC. :param _io.BinaryIO epsfile_obj: @@ -443,9 +455,9 @@ def read_records_before_error_section(epsfile_obj, mdr_descriptor, mdr_class_off """ data = {} for _, row in mdr_descriptor.iterrows(): - epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) + epsfile_mmap_subset = epsfile_mmap[mdr_class_offset + int(row["OFFSET"]):] reshape_flag = False if row["FIELD"] in ["EARTH_LOCATION", "ANGULAR_RELATION"] else True - values = read_values(epsfile_obj, row, reshape=reshape_flag) + values = read_values(epsfile_mmap_subset, row, reshape=reshape_flag) data[row["FIELD"]] = values return data @@ -460,7 +472,7 @@ def datetime_to_second_since_2000(date): return (date - era).total_seconds() -def read_nerr_values(epsfile_obj, mdr_descriptor, mdr_class_offset): +def read_nerr_values(epsfile_mmap, mdr_descriptor, mdr_class_offset): """FIXME DOC. Return a numpy array of all NERR values as red from each MDR class (i.e. data row). The @@ -474,21 +486,23 @@ def read_nerr_values(epsfile_obj, mdr_descriptor, mdr_class_offset): """ nerr_row = mdr_descriptor.loc[mdr_descriptor["FIELD"] == "NERR"].iloc[0] nerr_values = [] - epsfile_obj.seek(mdr_class_offset) + epsfile_mmap_subset = epsfile_mmap[mdr_class_offset:] while True: - grh = epsnative_reader.grh_reader(epsfile_obj) + grh = epsnative_reader.grh_reader(epsfile_mmap_subset) if grh: if grh[0:3] == ("mdr", 1, 4): - epsfile_obj.seek(mdr_class_offset + int(nerr_row["OFFSET"])) - nerr_values.append(read_values(epsfile_obj, nerr_row, reshape=False)[0]) + new_offset = mdr_class_offset + int(nerr_row["OFFSET"]) + epsfile_mmap_subset = epsfile_mmap[new_offset:] + nerr_values.append(read_values( + epsfile_mmap_subset, nerr_row, reshape=False)[0]) mdr_class_offset += grh[3] - epsfile_obj.seek(mdr_class_offset) + epsfile_mmap_subset = epsfile_mmap[mdr_class_offset:] else: break return np.array(nerr_values) -def read_all_rows(epsfile_obj, descriptor, mdr_class_offset, sensing_start, sensing_stop): +def read_all_rows(epsfile_mmap, descriptor, mdr_class_offset, sensing_start, sensing_stop): """FIXME DOC. :param _io.BinaryIO epsfile_obj: @@ -502,13 +516,13 @@ def read_all_rows(epsfile_obj, descriptor, mdr_class_offset, sensing_start, sens last_constant_row = mdr_descriptor.loc[mdr_descriptor["FIELD"] == "ERROR_DATA_INDEX"].index[0] mdr_descr_constant_offsets = mdr_descriptor[1:last_constant_row + 1] vars_must_be_extend_to_fov = get_vars_must_be_extend_to_fov(mdr_descriptor) - max_nerr = read_nerr_values(epsfile_obj, mdr_descriptor, mdr_class_offset).max() + max_nerr = read_nerr_values(epsfile_mmap, mdr_descriptor, mdr_class_offset).max() algorithms_data = [] data_before_errors_section = [] errors_data = [] - epsfile_obj.seek(mdr_class_offset) + epsfile_mmap_subset = epsfile_mmap[mdr_class_offset:] while True: - grh = epsnative_reader.grh_reader(epsfile_obj) + grh = epsnative_reader.grh_reader(epsfile_mmap_subset) if grh: if grh[0:3] == ("mdr", 1, 4): overlap = reckon_overlap( @@ -519,13 +533,13 @@ def read_all_rows(epsfile_obj, descriptor, mdr_class_offset, sensing_start, sens record_stop_time = datetime_to_second_since_2000(grh[-1]) records_before_error_section = read_records_before_error_section( - epsfile_obj, mdr_descr_constant_offsets, mdr_class_offset + epsfile_mmap, mdr_descr_constant_offsets, mdr_class_offset ) records_before_error_section["record_start_time"] = record_start_time records_before_error_section["record_stop_time"] = record_stop_time data_before_errors_section.append(records_before_error_section) algorithm_data, errors = read_algorithm_sections( - epsfile_obj, + epsfile_mmap, mdr_descriptor, mdr_class_offset, vars_must_be_extend_to_fov, @@ -538,7 +552,7 @@ def read_all_rows(epsfile_obj, descriptor, mdr_class_offset, sensing_start, sens mdr_class_offset += grh[3] else: algorithms_data.append("dummy_mdr") - epsfile_obj.seek(mdr_class_offset) + epsfile_mmap_subset = epsfile_mmap[mdr_class_offset:] else: break return data_before_errors_section, algorithms_data, errors_data @@ -554,8 +568,9 @@ def read_product_data(epsfile_obj, descriptor, mdr_class_offset, sensing_start, :param datetime.datetime sensing_stop: :return (dict, dict, dict): """ + epsfile_mmap = np.memmap(epsfile_obj, mode="r") data_before_errors, algorithms_data, errors = read_all_rows( - epsfile_obj, descriptor, mdr_class_offset, sensing_start, sensing_stop + epsfile_mmap, descriptor, mdr_class_offset, sensing_start, sensing_stop ) algorithms_data = fill_invalid_rows(algorithms_data) stacked_algo_data = initialize_stacked_algorithms_output(algorithms_data) diff --git a/satpy/tests/conftest.py b/satpy/tests/conftest.py index 501fb49530..b2d061a833 100644 --- a/satpy/tests/conftest.py +++ b/satpy/tests/conftest.py @@ -24,6 +24,7 @@ import pathlib import shutil +import numpy as np import pytest import requests @@ -66,13 +67,17 @@ def include_test_etc(): @pytest.fixture(scope="session") -def sample_file(tmp_path_factory): +def sample_iasisndl2(tmp_path_factory, request): """Obtain sample file.""" fn = pathlib.Path("/media/nas/x21308/IASI/IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") - if fn.exists(): + if not fn.exists(): + fn = tmp_path_factory.mktemp("data") / "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat" + data = requests.get(_url_sample_file, stream=True) + with fn.open(mode="wb") as fp: + shutil.copyfileobj(data.raw, fp) + if request.param == "string": return fn - fn = tmp_path_factory.mktemp("data") / "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat" - data = requests.get(_url_sample_file, stream=True) - with fn.open(mode="wb") as fp: - shutil.copyfileobj(data.raw, fp) - return fn + if request.param == "file": + return open(fn, mode="rb") + if request.param == "mmap": + return np.memmap(fn, mode="r", offset=0) diff --git a/satpy/tests/reader_tests/conftest.py b/satpy/tests/reader_tests/conftest.py index ca9e3fc66e..bee8be394b 100644 --- a/satpy/tests/reader_tests/conftest.py +++ b/satpy/tests/reader_tests/conftest.py @@ -17,6 +17,13 @@ # satpy. If not, see . """Setup and configuration for all reader tests.""" +import pathlib +import shutil + +import numpy as np +import pytest +import requests + from ._modis_fixtures import ( modis_l1b_imapp_1000m_file, modis_l1b_imapp_geo_file, @@ -33,3 +40,24 @@ modis_l2_nasa_mod35_file, modis_l2_nasa_mod35_mod03_files, ) + +_url_sample_file = ( + "https://go.dwd-nextcloud.de/index.php/s/z87KfL72b9dM5xm/download/" + "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") + + +@pytest.fixture(scope="session") +def iasisndl2_file(tmp_path_factory, request): + """Obtain sample file.""" + fn = pathlib.Path("/media/nas/x21308/IASI/IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") + if not fn.exists(): + fn = tmp_path_factory.mktemp("data") / "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat" + data = requests.get(_url_sample_file, stream=True) + with fn.open(mode="wb") as fp: + shutil.copyfileobj(data.raw, fp) + if request.param == "string": + return fn + if request.param == "file": + return open(fn, mode="rb") + if request.param == "mmap": + return np.memmap(fn, mode="r", offset=0) diff --git a/satpy/tests/reader_tests/test_epsnative_reader.py b/satpy/tests/reader_tests/test_epsnative_reader.py index 6acf9b57a0..2e32f28d71 100644 --- a/satpy/tests/reader_tests/test_epsnative_reader.py +++ b/satpy/tests/reader_tests/test_epsnative_reader.py @@ -40,7 +40,6 @@ import datetime import os -import numpy as np import pytest from satpy.readers import epsnative_reader @@ -105,17 +104,13 @@ def test_assemble_descriptor(): assert set(descriptor.keys()) == exp_keys -@pytest.mark.parametrize("how", ["string", "file", "mmap"]) -def test_grh_reader(sample_file, how): +@pytest.mark.parametrize( + "iasisndl2_file", + ["string", "file", "mmap"], + indirect=["iasisndl2_file"]) +def test_grh_reader(iasisndl2_file): """FIXME DOC.""" - if how == "file": - sample = open(sample_file, mode="rb") - elif how == "mmap": - sample = np.memmap(sample_file, mode="r", offset=0) - else: - sample = sample_file - - grh = epsnative_reader.grh_reader(sample) + grh = epsnative_reader.grh_reader(iasisndl2_file) assert len(grh) == 6 assert grh[0] == "mphr" assert grh[1] == 0 diff --git a/satpy/tests/reader_tests/test_iasisnd02.py b/satpy/tests/reader_tests/test_iasisnd02.py index d4ebf9205c..c11159653c 100644 --- a/satpy/tests/reader_tests/test_iasisnd02.py +++ b/satpy/tests/reader_tests/test_iasisnd02.py @@ -40,16 +40,34 @@ import dask import numpy as np +import pandas as pd +import pytest from ..utils import CustomScheduler +sample_file_str = pytest.mark.parametrize( + "iasisndl2_file", + ["string"], + indirect=["iasisndl2_file"]) -def test_read_giadr(sample_file): +sample_file_file = pytest.mark.parametrize( + "iasisndl2_file", + ["file"], + indirect=["iasisndl2_file"]) + +sample_file_mmap = pytest.mark.parametrize( + "iasisndl2_file", + ["mmap"], + indirect=["iasisndl2_file"]) + + +@sample_file_str +def test_read_giadr(iasisndl2_file): """Test reading GIADR.""" from satpy.readers.epsnative_reader import assemble_descriptor from satpy.readers.iasi_l2_eps import read_giadr descriptor = assemble_descriptor("IASISND02") - class_data = read_giadr(sample_file, descriptor) + class_data = read_giadr(iasisndl2_file, descriptor) assert len(class_data) == 19 assert class_data["PRESSURE_LEVELS_OZONE"]["units"] == "Pa" @@ -67,7 +85,8 @@ def test_datetime_to_second_since_2000(): assert sec == 60 -def test_read_all_rows(sample_file): +@sample_file_mmap +def test_read_all_rows(iasisndl2_file): """Test reading all rows from data.""" from satpy.readers.epsnative_reader import assemble_descriptor from satpy.readers.iasi_l2_eps import read_all_rows @@ -77,10 +96,9 @@ def test_read_all_rows(sample_file): descriptor = assemble_descriptor("IASISND02") mdr_class_offset = 271719118 # last row - with open(sample_file, "rb") as epsfile_obj: - data_before_errors_section, algorithms_data, errors_data = read_all_rows( - epsfile_obj, descriptor, mdr_class_offset, sensing_start, sensing_stop - ) + data_before_errors_section, algorithms_data, errors_data = read_all_rows( + iasisndl2_file, descriptor, mdr_class_offset, sensing_start, sensing_stop + ) assert len(data_before_errors_section) == 1 assert data_before_errors_section[0]["INTEGRATED_CO"].min() == 7205 @@ -91,21 +109,81 @@ def test_read_all_rows(sample_file): assert len(errors_data) == 1 assert errors_data[0]["SURFACE_Z"].min() == 30 + np.testing.assert_array_equal( + data_before_errors_section[0]["ATMOSPHERIC_TEMPERATURE"][50, 50:55], + np.array([22859, 22844, 22848, 22816, 22812], dtype=np.uint16)) -def test_read_product_data(sample_file): + +@sample_file_mmap +def test_read_nerr_values(iasisndl2_file): + """Test reading the number of uncertainty values.""" + from satpy.readers.epsnative_reader import assemble_descriptor + from satpy.readers.iasi_l2_eps import read_nerr_values + + descriptor = assemble_descriptor("IASISND02")[("mdr", 1, 4)] + nerr = read_nerr_values(iasisndl2_file, descriptor, 271719118) + np.testing.assert_array_equal(nerr, [69]) + + +@sample_file_mmap +def test_read_values(iasisndl2_file): + """Test reading values from a mmap.""" + from satpy.readers.iasi_l2_eps import read_values + + row = pd.Series( + np.array( + ['NERR', 'Number of error data records for current scan line', + 0.0, np.nan, 1, '1', 1, 'u-byte', 1, 1.0, 207747.0], + dtype=object), + index=pd.Index( + ['FIELD', 'DESCRIPTION', 'SF', 'UNITS', 'DIM1', 'DIM2', 'DIM3', + 'TYPE', 'TYPE_SIZE', 'FIELD_SIZE', 'OFFSET'], + dtype=object), + name=52) + try: + iasisndl2_file.seek(271926865) + except AttributeError: + iasisndl2_file = iasisndl2_file[271926865:] + + vals = read_values(iasisndl2_file, row, False) + np.testing.assert_array_equal(vals, [69]) + + +@sample_file_mmap +def test_read_records_before_error_section(iasisndl2_file): + """Test reading records before the error section.""" + from satpy.readers.epsnative_reader import assemble_descriptor + from satpy.readers.iasi_l2_eps import read_records_before_error_section + + descriptor = assemble_descriptor("IASISND02")[("mdr", 1, 4)] + data = read_records_before_error_section( + iasisndl2_file, + descriptor[1:54], + 271719118) + np.testing.assert_array_equal( + data["SURFACE_TEMPERATURE"][48:55], + np.array([29209, 29276, 29386, 29220, 29266, 29568, 29302], + dtype=np.uint16)) + np.testing.assert_array_equal( + data["ATMOSPHERIC_TEMPERATURE"][20:25, 50], + np.array([23917, 23741, 23583, 23436, 23290], + dtype=np.uint16)) + + +@sample_file_file +def test_read_product_data(iasisndl2_file): """Test reading product data.""" from satpy.readers.epsnative_reader import assemble_descriptor from satpy.readers.iasi_l2_eps import read_product_data descriptor = assemble_descriptor("IASISND02") mdr_class_offset = 260002007 # not quite the last row - with open(sample_file, "rb") as epsfile_obj: - (sdbs, sad, sed) = read_product_data( - epsfile_obj, - descriptor, - mdr_class_offset, - datetime.datetime(2019, 6, 5, 0, 23, 52), - datetime.datetime(2019, 6, 5, 2, 8, 56)) + (sdbs, sad, sed) = read_product_data( + iasisndl2_file, + descriptor, + mdr_class_offset, + datetime.datetime(2019, 6, 5, 0, 23, 52), + datetime.datetime(2019, 6, 5, 2, 8, 56)) assert sdbs["SURFACE_TEMPERATURE"].shape == (37, 120) assert sdbs["ATMOSPHERIC_TEMPERATURE"].shape == (37, 101, 120) np.testing.assert_array_equal( @@ -117,24 +195,22 @@ def test_read_product_data(sample_file): assert isinstance(sdbs["ATMOSPHERIC_TEMPERATURE"], dask.array.Array) -def test_load(sample_file, tmp_path): +@sample_file_str +def test_load(iasisndl2_file, tmp_path): """Test read at a scene and write again.""" from satpy import Scene with dask.config.set(scheduler=CustomScheduler(max_computes=0)): - sc = Scene(filenames=[sample_file], reader=["iasi_l2_eps2"]) + sc = Scene(filenames=[iasisndl2_file], reader=["iasi_l2_eps2"]) sc.load(["surface_temperature"]) assert sc["surface_temperature"].dims == ("y", "x") np.testing.assert_allclose( - sc["surface_temperature"][0, 104:110], - np.array([262.55, 261.83, 265.85, 265.81, 256.15, 255.23])) - np.testing.assert_allclose( - sc["surface_temperature"][1, 104:110], - np.array([261.9, 266.88, 268.47, 266.76, 256.63, 261.32])) + sc["surface_temperature"][0, 100:104], + np.array([270.39, 270.21, 269.65, 269.66])) np.testing.assert_allclose( - sc["surface_temperature"][2, 108:112], - np.array([262.33, 265.42, 268.32, 266.21])) + sc["surface_temperature"][1, 106:110], + np.array([269.91, 270.07, 270.02, 270.41])) np.testing.assert_allclose( - sc["surface_temperature"][10, 83:89], - np.array([271.11, 280.07, 278.65, 278.84, 278.37, 281.54])) + sc["surface_temperature"][30, 60:66], + np.array([282.27, 283.18, 285.67, 282.98, 282.81, 282.9])) - assert isinstance(sc["surface_temperature"], dask.array.Array) + assert isinstance(sc["surface_temperature"].data, dask.array.Array) From 345b6cc79ae1bfa0f5f05d45d113da80ac64d0a1 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Wed, 12 Jul 2023 15:38:35 +0200 Subject: [PATCH 16/23] Remove dead code --- satpy/readers/iasi_l2_eps.py | 11 +++-------- satpy/tests/conftest.py | 26 -------------------------- 2 files changed, 3 insertions(+), 34 deletions(-) diff --git a/satpy/readers/iasi_l2_eps.py b/satpy/readers/iasi_l2_eps.py index e21abec2d5..2a87844077 100644 --- a/satpy/readers/iasi_l2_eps.py +++ b/satpy/readers/iasi_l2_eps.py @@ -180,7 +180,7 @@ def read_giadr(input_product, descriptor, **kwargs): return class_data -def set_values_in_mdr_descriptor(epsfile_obj, mdr_descriptor, mdr_class_offset, field_name): +def set_values_in_mdr_descriptor(epsfile_mmap, mdr_descriptor, mdr_class_offset, field_name): """FIXME DOC. :param _io.BinaryIO epsfile_obj: @@ -193,9 +193,8 @@ def set_values_in_mdr_descriptor(epsfile_obj, mdr_descriptor, mdr_class_offset, row_idx = row.index[0] row = row.squeeze() # read values -# epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) - epsfile_obj = epsfile_obj[mdr_class_offset + int(row["OFFSET"]):] - value = read_values(epsfile_obj, row).astype(int)[0] + epsfile_mmap_subset = epsfile_mmap[mdr_class_offset + int(row["OFFSET"]):] + value = read_values(epsfile_mmap_subset, row).astype(int)[0] # set the read values in MDR-descriptor mdr_descriptor.loc[mdr_descriptor["DIM2"] == field_name.lower(), "DIM2"] = value return value, row_idx @@ -265,8 +264,6 @@ def read_errors(epsfile_obj, mdr_descriptor, mdr_class_offset, max_nerr): for field in fields_to_read: row = mdr_descriptor.loc[mdr_descriptor["FIELD"] == field].squeeze() epsfile_mmap_subset = epsfile_obj[mdr_class_offset + int(row["OFFSET"]):] -# epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) -# values = read_values(epsfile_obj, row, reshape=True) values = read_values(epsfile_mmap_subset, row, reshape=True) if field != "SURFACE_Z": if row["DIM2"] == 0: @@ -311,9 +308,7 @@ def read_algorithm_sections( section = {} for _, row in section_rows.iterrows(): epsfile_mmap_subset = epsfile_obj[mdr_class_offset + int(row["OFFSET"]):] -# epsfile_obj.seek(mdr_class_offset + int(row["OFFSET"])) values = read_values(epsfile_mmap_subset, row, reshape=True) -# values = read_values(epsfile_obj, row, reshape=True) if row["FIELD"] in vars_must_be_extend_to_fov: values = values.reshape(-1, 1) if values.ndim == 1 else values values = np.pad( diff --git a/satpy/tests/conftest.py b/satpy/tests/conftest.py index b2d061a833..8bcbea2093 100644 --- a/satpy/tests/conftest.py +++ b/satpy/tests/conftest.py @@ -21,12 +21,8 @@ """ import os -import pathlib -import shutil -import numpy as np import pytest -import requests import satpy @@ -59,25 +55,3 @@ def include_test_etc(): """Tell Satpy to use the config 'etc' directory from the tests directory.""" with satpy.config.set(config_path=[TEST_ETC_DIR]): yield TEST_ETC_DIR - - -_url_sample_file = ( - "https://go.dwd-nextcloud.de/index.php/s/z87KfL72b9dM5xm/download/" - "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") - - -@pytest.fixture(scope="session") -def sample_iasisndl2(tmp_path_factory, request): - """Obtain sample file.""" - fn = pathlib.Path("/media/nas/x21308/IASI/IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat") - if not fn.exists(): - fn = tmp_path_factory.mktemp("data") / "IASI_SND_02_M01_20190605002352Z_20190605020856Z_N_O_20190605011702Z.nat" - data = requests.get(_url_sample_file, stream=True) - with fn.open(mode="wb") as fp: - shutil.copyfileobj(data.raw, fp) - if request.param == "string": - return fn - if request.param == "file": - return open(fn, mode="rb") - if request.param == "mmap": - return np.memmap(fn, mode="r", offset=0) From e0a53ea1c1cc79d6d0ef8e29d9b93040e9c761e0 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Wed, 12 Jul 2023 15:56:52 +0200 Subject: [PATCH 17/23] Change sample data for test_eps_native_reader Remove more dead code and change the sample data used for test_eps_native_reader. Change the expected results accordilgny. --- .../reader_tests/test_epsnative_reader.py | 134 +++++------------- 1 file changed, 33 insertions(+), 101 deletions(-) diff --git a/satpy/tests/reader_tests/test_epsnative_reader.py b/satpy/tests/reader_tests/test_epsnative_reader.py index 2e32f28d71..84d722ec43 100644 --- a/satpy/tests/reader_tests/test_epsnative_reader.py +++ b/satpy/tests/reader_tests/test_epsnative_reader.py @@ -72,8 +72,10 @@ ), ] -SAMPLE_PATH = os.path.join(os.path.dirname(__file__), "sample_data") -SAMPLE_DATA = os.path.join(SAMPLE_PATH, "AVHRRL1.nat") +sample_file_str = pytest.mark.parametrize( + "iasisndl2_file", + ["string"], + indirect=["iasisndl2_file"]) def test_get_class_tuple(): @@ -85,29 +87,15 @@ def test_get_class_tuple(): assert epsnative_reader.get_class_tuple(class_string) == ("class", "dummy_dummy2") -def test_csv_list_by_product(): - """FIXME DOC.""" - csv_list = epsnative_reader.csv_list_by_product("AVHRRL1") - exp_csv = ["giadr_1_3.csv", "giadr_2_2.csv", "mdr_2_5.csv", "mphr_0_2.csv"] - - assert len(csv_list) == 4 - for csv_path in csv_list: - assert os.path.isfile(csv_path) - assert os.path.basename(csv_path) in exp_csv - - def test_assemble_descriptor(): """FIXME DOC.""" - descriptor = epsnative_reader.assemble_descriptor("AVHRRL1") - exp_keys = {("giadr", 1, 3), ("giadr", 2, 2), ("mdr", 2, 5), ("mphr", 0, 2)} + descriptor = epsnative_reader.assemble_descriptor("IASISNDL2") + exp_keys = {('mphr', 0, 2), ('giadr', 1, 4), ('mdr', 1, 4)} assert set(descriptor.keys()) == exp_keys -@pytest.mark.parametrize( - "iasisndl2_file", - ["string", "file", "mmap"], - indirect=["iasisndl2_file"]) +@sample_file_str def test_grh_reader(iasisndl2_file): """FIXME DOC.""" grh = epsnative_reader.grh_reader(iasisndl2_file) @@ -120,82 +108,47 @@ def test_grh_reader(iasisndl2_file): assert grh[5] == datetime.datetime(2019, 6, 5, 2, 8, 56, 225000) -def test_find_mphr_csv(): +@sample_file_str +def test_mphr_reader(iasisndl2_file): """FIXME DOC.""" - mphr_csv = epsnative_reader.find_mphr_csv() - - assert os.path.isfile(mphr_csv) - - -def test_mphr_reader(): - """FIXME DOC.""" - mphr_content = epsnative_reader.mphr_reader(SAMPLE_DATA) + mphr_content = epsnative_reader.mphr_reader(iasisndl2_file) assert len(mphr_content) == 72 - assert mphr_content["INSTRUMENT_ID"] == "AVHR" - assert mphr_content["ORBIT_START"] == 58431 - assert mphr_content["ORBIT_END"] == 58432 + assert mphr_content["INSTRUMENT_ID"] == "IASI" + assert mphr_content["ORBIT_START"] == 34826 + assert mphr_content["ORBIT_END"] == 34827 -def test_first_class_occurrence(): +@sample_file_str +def test_first_class_occurrence(iasisndl2_file): """FIXME DOC.""" class_name = "mphr" - grh, offset = epsnative_reader.first_class_occurrence(SAMPLE_DATA, class_name) + grh, offset = epsnative_reader.first_class_occurrence(iasisndl2_file, class_name) assert offset == 0 class_name = "giadr" - grh, offset = epsnative_reader.first_class_occurrence(SAMPLE_DATA, class_name) - assert offset == 4587 + grh, offset = epsnative_reader.first_class_occurrence(iasisndl2_file, class_name) + assert offset == 3361 assert grh[0] == class_name assert grh[1] == 1 - assert grh[2] == 3 + assert grh[2] == 4 class_name = "mdr" - grh, offset = epsnative_reader.first_class_occurrence(SAMPLE_DATA, class_name) - assert offset == 5077 + grh, offset = epsnative_reader.first_class_occurrence(iasisndl2_file, class_name) + assert offset == 4818 assert grh[0] == class_name - assert grh[1] == 2 + assert grh[1] == 1 assert grh[2] == 4 -def test_read_ipr_sequence(): +@sample_file_str +def test_read_ipr_sequence(iasisndl2_file): """FIXME DOC.""" - ipr_sequence = epsnative_reader.read_ipr_sequence(SAMPLE_DATA) + ipr_sequence = epsnative_reader.read_ipr_sequence(iasisndl2_file) - assert len(ipr_sequence) == 11 - assert ipr_sequence[-1]["class"] == ("mdr", 2) - assert ipr_sequence[-1]["offset"] == 5077 - - -def test_read_grh_of_target_class(): - """FIXME DOC.""" - target_offset = 0 - with open(SAMPLE_DATA, "rb") as eps_fileobj: - target_grh = epsnative_reader.read_grh_of_target_class(eps_fileobj, target_offset) - assert target_grh[0] == "mphr" - assert target_grh[1] == 0 - assert target_grh[2] == 2 - assert target_grh[3] == 3307 - - target_offset = 5077 - with open(SAMPLE_DATA, "rb") as eps_fileobj: - target_grh = epsnative_reader.read_grh_of_target_class(eps_fileobj, target_offset) - assert target_grh[0] == "mdr" - assert target_grh[1] == 2 - assert target_grh[2] == 4 - assert target_grh[3] == 26660 - - -def test_add_info_about_class(): - """FIXME DOC.""" - mphr = epsnative_reader.mphr_reader(SAMPLE_DATA) - current_class = {"class": "", "offset": 5077} - with open(SAMPLE_DATA, "rb") as eps_fileobj: - current_class = epsnative_reader.add_info_about_class(eps_fileobj, current_class, {}, mphr) - - assert current_class["class_id"] == ("mdr", 2, 4) - assert current_class["class_size"] == 26660 - assert current_class["nr_records"] == 10.0 + assert len(ipr_sequence) == 2 + assert ipr_sequence[-1]["class"] == ("mdr", 1) + assert ipr_sequence[-1]["offset"] == 4818 def test_reckon_dtype(): @@ -209,30 +162,9 @@ def test_reckon_dtype(): def test_bands_to_records_reader(): """FIXME DOC.""" - records = epsnative_reader.bands_to_records_reader("AVHRRL1") - - assert len(records) == 7 - assert records["channel_5"]["record_name"] == "SCENE_RADIANCES" - assert records["channel_5"]["band_position"] == 4 - assert records["channel_5"]["convention"] == "BIL" - - -def test_add_record_info_to_band(): - """FIXME DOC.""" - record_info = epsnative_reader.add_record_info_to_band("AVHRRL1") - - assert len(record_info) == 7 - assert record_info[0]["band_id"] == "channel_1" - assert record_info[0]["shape"] == [2048, 5] - assert record_info[-1]["band_id"] == "longitude" - assert record_info[-1]["shape"] == [2, 103] - - -def test_create_toc(): - """FIXME DOC.""" - toc = epsnative_reader.create_toc(SAMPLE_DATA) + records = epsnative_reader.bands_to_records_reader("IASISNDL2") - assert len(toc) == 11 - assert toc[0]["class_id"] == ("geadr", 1, 1) - assert toc[-1]["class_id"] == ("mdr", 2, 4) - assert toc[-1]["y_offset"] == 0 + assert len(records) == 74 + assert records["atmospheric_temperature"]["record_name"] == "ATMOSPHERIC_TEMPERATURE" + assert records["solar_azimuth"]["metadata"]["units"] == "degrees" + assert records["surface_emissivity"]["metadata"]["scale_factor"] == 4 From 3bdb5248e2d854d33bd95295421f95fc09b2292c Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Wed, 12 Jul 2023 16:17:20 +0200 Subject: [PATCH 18/23] swap readers --- satpy/etc/readers/iasi_l2_eps.yaml | 2 +- satpy/tests/reader_tests/test_iasisnd02.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/satpy/etc/readers/iasi_l2_eps.yaml b/satpy/etc/readers/iasi_l2_eps.yaml index 315770a36f..160e060c70 100644 --- a/satpy/etc/readers/iasi_l2_eps.yaml +++ b/satpy/etc/readers/iasi_l2_eps.yaml @@ -15,7 +15,7 @@ reader: file_types: iasi_l2_eps: - file_reader: !!python/name:satpy.readers.eps_l2.EPSIASIL2FileHandler + file_reader: !!python/name:satpy.readers.iasi_l2_eps.EPSIASIL2FileHandler # see https://www.eumetsat.int/media/40048 §7 for file pattern file_patterns: - IASI_SND_02_{platform_short_name}_{start_time:%Y%m%d%H%M%SZ}_{end_time:%Y%m%d%H%M%SZ}_{processing_mode}_{disposition_mode}_{creation_time:%Y%m%d%H%M%SZ} diff --git a/satpy/tests/reader_tests/test_iasisnd02.py b/satpy/tests/reader_tests/test_iasisnd02.py index c11159653c..39d9fe04cc 100644 --- a/satpy/tests/reader_tests/test_iasisnd02.py +++ b/satpy/tests/reader_tests/test_iasisnd02.py @@ -200,7 +200,7 @@ def test_load(iasisndl2_file, tmp_path): """Test read at a scene and write again.""" from satpy import Scene with dask.config.set(scheduler=CustomScheduler(max_computes=0)): - sc = Scene(filenames=[iasisndl2_file], reader=["iasi_l2_eps2"]) + sc = Scene(filenames=[iasisndl2_file], reader=["iasi_l2_eps"]) sc.load(["surface_temperature"]) assert sc["surface_temperature"].dims == ("y", "x") np.testing.assert_allclose( From 369ea3381b626b2b7453b4f73e8b8a4579028068 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Fri, 21 Jul 2023 11:16:34 +0200 Subject: [PATCH 19/23] Add geolocation to IASISND02 reader --- satpy/readers/iasi_l2_eps.py | 11 +++++++---- satpy/tests/reader_tests/test_iasisnd02.py | 14 +++++++++++++- 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/satpy/readers/iasi_l2_eps.py b/satpy/readers/iasi_l2_eps.py index 2a87844077..d45ffb53c0 100644 --- a/satpy/readers/iasi_l2_eps.py +++ b/satpy/readers/iasi_l2_eps.py @@ -61,9 +61,11 @@ def get_dataset(self, dataid, dataset_info): """Get dataset.""" if self._nc is None: self._nc = self._get_netcdf_dataset() - da = self._nc[dataid["name"]] - da = da * da.attrs.pop("scale_factor", 1) - return da + data = self._nc[dataid["name"]] + with xr.set_options(keep_attrs=True): + data = xr.where(data != np.iinfo(data.dtype).max, data, np.nan) + data = data * data.attrs.pop("scale_factor", 1) + return data def _get_netcdf_dataset(self): """Get full NetCDF dataset.""" @@ -84,7 +86,8 @@ def available_datasets(self, configured_datasets): """Get available datasets.""" # FIXME: do this without converting/reading the file — maybe hardcode # still? - common = {"file_type": "iasi_l2_eps", "resolution": 12000} + common = {"file_type": "iasi_l2_eps", "resolution": 12000, + "coordinates": ["lon", "lat"]} if self._nc is None: self._nc = self._get_netcdf_dataset() for var in self._nc.data_vars: diff --git a/satpy/tests/reader_tests/test_iasisnd02.py b/satpy/tests/reader_tests/test_iasisnd02.py index 39d9fe04cc..845ca65a65 100644 --- a/satpy/tests/reader_tests/test_iasisnd02.py +++ b/satpy/tests/reader_tests/test_iasisnd02.py @@ -41,6 +41,7 @@ import dask import numpy as np import pandas as pd +import pyresample import pytest from ..utils import CustomScheduler @@ -201,7 +202,8 @@ def test_load(iasisndl2_file, tmp_path): from satpy import Scene with dask.config.set(scheduler=CustomScheduler(max_computes=0)): sc = Scene(filenames=[iasisndl2_file], reader=["iasi_l2_eps"]) - sc.load(["surface_temperature"]) + sc.load(["surface_temperature", "atmospheric_temperature", + "atmospheric_water_vapour"]) assert sc["surface_temperature"].dims == ("y", "x") np.testing.assert_allclose( sc["surface_temperature"][0, 100:104], @@ -212,5 +214,15 @@ def test_load(iasisndl2_file, tmp_path): np.testing.assert_allclose( sc["surface_temperature"][30, 60:66], np.array([282.27, 283.18, 285.67, 282.98, 282.81, 282.9])) + np.testing.assert_array_equal( + sc["surface_temperature"][30, :5], + [np.nan]*5) + np.testing.assert_array_equal( + sc["atmospheric_water_vapour"][0, 0, :5], + [np.nan]*5) assert isinstance(sc["surface_temperature"].data, dask.array.Array) + assert isinstance(sc["surface_temperature"].attrs["area"], + pyresample.SwathDefinition) + assert sc["surface_temperature"].attrs["standard_name"] == "surface_temperature" + assert sc["surface_temperature"].attrs["units"] == "K" From f9b400aeddbccb07902d3185d034c6b0497ce2fe Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Fri, 21 Jul 2023 11:42:35 +0200 Subject: [PATCH 20/23] Change name of private method --- satpy/readers/iasi_l2_eps.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/satpy/readers/iasi_l2_eps.py b/satpy/readers/iasi_l2_eps.py index d45ffb53c0..b9c0891ea7 100644 --- a/satpy/readers/iasi_l2_eps.py +++ b/satpy/readers/iasi_l2_eps.py @@ -52,23 +52,22 @@ class EPSIASIL2FileHandler(BaseFileHandler): """ # TODO: - # - make dask-friendly / lazy - # - only read variables that are requested + # - make it faster _nc = None def get_dataset(self, dataid, dataset_info): """Get dataset.""" if self._nc is None: - self._nc = self._get_netcdf_dataset() + self._nc = self._get_xarray_dataset() data = self._nc[dataid["name"]] with xr.set_options(keep_attrs=True): data = xr.where(data != np.iinfo(data.dtype).max, data, np.nan) data = data * data.attrs.pop("scale_factor", 1) return data - def _get_netcdf_dataset(self): - """Get full NetCDF dataset.""" + def _get_xarray_dataset(self): + """Get full xarray dataset.""" input_product = self.filename descriptor = epsnative_reader.assemble_descriptor("IASISND02") ipr_sequence = epsnative_reader.read_ipr_sequence(input_product) @@ -84,12 +83,10 @@ def _get_netcdf_dataset(self): def available_datasets(self, configured_datasets): """Get available datasets.""" - # FIXME: do this without converting/reading the file — maybe hardcode - # still? common = {"file_type": "iasi_l2_eps", "resolution": 12000, "coordinates": ["lon", "lat"]} if self._nc is None: - self._nc = self._get_netcdf_dataset() + self._nc = self._get_xarray_dataset() for var in self._nc.data_vars: yield (True, {"name": var} | common | self._nc[var].attrs) From 7fce13ae8abdd2237d0ece94727138f9bb266b13 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Wed, 26 Jul 2023 11:50:46 +0200 Subject: [PATCH 21/23] Don't add coordinates where we shouldn't In the dynamic available datasets for the iasil2snd reader, don't add coordinates to variables that don't have corresponding dimensions. --- satpy/readers/iasi_l2_eps.py | 14 +++++++++++--- satpy/tests/reader_tests/test_iasisnd02.py | 3 ++- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/satpy/readers/iasi_l2_eps.py b/satpy/readers/iasi_l2_eps.py index b9c0891ea7..4b654aca70 100644 --- a/satpy/readers/iasi_l2_eps.py +++ b/satpy/readers/iasi_l2_eps.py @@ -28,6 +28,7 @@ import collections import datetime import itertools +import logging import dask.array as da import numpy as np @@ -36,6 +37,8 @@ from . import epsnative_reader from .file_handlers import BaseFileHandler +logger = logging.getLogger(__name__) + class EPSIASIL2FileHandler(BaseFileHandler): """EPS level 2 reader for IASI data. @@ -68,6 +71,7 @@ def get_dataset(self, dataid, dataset_info): def _get_xarray_dataset(self): """Get full xarray dataset.""" + logger.debug(f"Pre-loading {self.filename!s}") input_product = self.filename descriptor = epsnative_reader.assemble_descriptor("IASISND02") ipr_sequence = epsnative_reader.read_ipr_sequence(input_product) @@ -83,12 +87,16 @@ def _get_xarray_dataset(self): def available_datasets(self, configured_datasets): """Get available datasets.""" - common = {"file_type": "iasi_l2_eps", "resolution": 12000, - "coordinates": ["lon", "lat"]} + common = {"file_type": "iasi_l2_eps", "resolution": 12000} + coords = {"y": "lon", "x": "lat"} if self._nc is None: self._nc = self._get_xarray_dataset() for var in self._nc.data_vars: - yield (True, {"name": var} | common | self._nc[var].attrs) + extra = {} + coor = [coords[x] for x in self._nc[var].dims if x in coords] + if coor: + extra["coordinates"] = coor + yield (True, {"name": var} | common | extra | self._nc[var].attrs) missing_values = { diff --git a/satpy/tests/reader_tests/test_iasisnd02.py b/satpy/tests/reader_tests/test_iasisnd02.py index 845ca65a65..b83634d71d 100644 --- a/satpy/tests/reader_tests/test_iasisnd02.py +++ b/satpy/tests/reader_tests/test_iasisnd02.py @@ -203,7 +203,7 @@ def test_load(iasisndl2_file, tmp_path): with dask.config.set(scheduler=CustomScheduler(max_computes=0)): sc = Scene(filenames=[iasisndl2_file], reader=["iasi_l2_eps"]) sc.load(["surface_temperature", "atmospheric_temperature", - "atmospheric_water_vapour"]) + "atmospheric_water_vapour", "pressure_levels_temp"]) assert sc["surface_temperature"].dims == ("y", "x") np.testing.assert_allclose( sc["surface_temperature"][0, 100:104], @@ -226,3 +226,4 @@ def test_load(iasisndl2_file, tmp_path): pyresample.SwathDefinition) assert sc["surface_temperature"].attrs["standard_name"] == "surface_temperature" assert sc["surface_temperature"].attrs["units"] == "K" + assert "area" not in sc["pressure_levels_temp"].attrs From 802ceee354aa4c3c5933352fcf89ba2666952a08 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Wed, 23 Aug 2023 15:44:32 +0200 Subject: [PATCH 22/23] Add safeguards to prevent infinite loops Add safeguards in the IASI native reader to avoid going into infinite loops --- satpy/readers/iasi_l2_eps.py | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/satpy/readers/iasi_l2_eps.py b/satpy/readers/iasi_l2_eps.py index 4b654aca70..f4e4892a5a 100644 --- a/satpy/readers/iasi_l2_eps.py +++ b/satpy/readers/iasi_l2_eps.py @@ -202,6 +202,9 @@ def set_values_in_mdr_descriptor(epsfile_mmap, mdr_descriptor, mdr_class_offset, row = row.squeeze() # read values epsfile_mmap_subset = epsfile_mmap[mdr_class_offset + int(row["OFFSET"]):] + if epsfile_mmap_subset.size == 0: + raise ValueError(f"Could not read MDR for {row['FIELD']:s}. Found none " + "or with size zero. Corrupt file?") value = read_values(epsfile_mmap_subset, row).astype(int)[0] # set the read values in MDR-descriptor mdr_descriptor.loc[mdr_descriptor["DIM2"] == field_name.lower(), "DIM2"] = value @@ -493,12 +496,19 @@ def read_nerr_values(epsfile_mmap, mdr_descriptor, mdr_class_offset): while True: grh = epsnative_reader.grh_reader(epsfile_mmap_subset) if grh: + recsize = grh[3] if grh[0:3] == ("mdr", 1, 4): new_offset = mdr_class_offset + int(nerr_row["OFFSET"]) epsfile_mmap_subset = epsfile_mmap[new_offset:] + if epsfile_mmap_subset.size == 0: + raise ValueError(f'Unable to read {nerr_row["FIELD"]:s}. Corrupt file?') nerr_values.append(read_values( epsfile_mmap_subset, nerr_row, reshape=False)[0]) - mdr_class_offset += grh[3] + else: + # happens for IASI_SND_02_M01_20230724181453Z_20230724195652Z_N_O_20230724191952Z + # at least, and could cause infinite loop if not careful + logger.debug(f"Found unexpected or invalid record, skipping {recsize:d} bytes") + mdr_class_offset += recsize epsfile_mmap_subset = epsfile_mmap[mdr_class_offset:] else: break @@ -527,6 +537,7 @@ def read_all_rows(epsfile_mmap, descriptor, mdr_class_offset, sensing_start, sen while True: grh = epsnative_reader.grh_reader(epsfile_mmap_subset) if grh: + recsize = grh[3] if grh[0:3] == ("mdr", 1, 4): overlap = reckon_overlap( sensing_start, sensing_stop, grh[-2], grh[-1] @@ -550,11 +561,15 @@ def read_all_rows(epsfile_mmap, descriptor, mdr_class_offset, sensing_start, sen ) algorithms_data.append(algorithm_data) errors_data.append(errors) - mdr_class_offset += grh[3] - else: - mdr_class_offset += grh[3] + mdr_class_offset += recsize else: - algorithms_data.append("dummy_mdr") + logger.debug(f"Found something else of length {recsize:d}") + # this is probably salvageable, but there are too many places in + # the code where putting a dummy value has implications, so for + # now we skip the entire file + # happens for IASI_SND_02_M01_20230724181453Z_20230724195652Z_N_O_20230724191952Z + # at least, and could cause infinite loop if not careful + raise ValueError("Potential data corruption") epsfile_mmap_subset = epsfile_mmap[mdr_class_offset:] else: break From 54c6826725cc4ab1a93b016931ccb74c39935146 Mon Sep 17 00:00:00 2001 From: Gerrit Holl Date: Wed, 23 Aug 2023 18:26:26 +0200 Subject: [PATCH 23/23] Remove/restore unworking implementation attempts Remove implementation attempts that didn't work, restoring those files. --- satpy/etc/eps_iasil2_9.0.xml | 909 --------------------- satpy/readers/eps_base.py | 209 ----- satpy/readers/eps_l1b.py | 198 ++++- satpy/readers/eps_l2.py | 219 ----- satpy/readers/xmlformat.py | 96 +-- satpy/tests/reader_tests/test_eps_l1b.py | 6 +- satpy/tests/reader_tests/test_xmlformat.py | 78 -- 7 files changed, 227 insertions(+), 1488 deletions(-) delete mode 100644 satpy/etc/eps_iasil2_9.0.xml delete mode 100644 satpy/readers/eps_base.py delete mode 100644 satpy/readers/eps_l2.py delete mode 100644 satpy/tests/reader_tests/test_xmlformat.py diff --git a/satpy/etc/eps_iasil2_9.0.xml b/satpy/etc/eps_iasil2_9.0.xml deleted file mode 100644 index f5219776b5..0000000000 --- a/satpy/etc/eps_iasil2_9.0.xml +++ /dev/null @@ -1,909 +0,0 @@ - - - - - - - - 11 - 0 - current - - - 9 - 0 - - - EPS IASI Level 2 Format - - This IASI 2 description was generated using the IASI PFS Excel document Issue 6 Revision 3 - - - IASI_???_02_*Z* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Geolocation - Geolocation Coverage (Latitude, Longitude) - mdr[].EARTH_LOCATION[][0] - mdr[].EARTH_LOCATION[][1] - - - - PRESSURE - - - Atmospheric Temperature - Atmospheric Temperature - mdr[].ATMOSPHERIC_TEMPERATURE[][$PRESSURE] - Geolocation - - - - PRESSURE - - - Water Vapour - Atmospheric Water Vapour - mdr[].ATMOSPHERIC_WATER_VAPOUR[][$PRESSURE] - Geolocation - - - - PRESSURE - - - Ozone - Atmospheric Ozone - mdr[].ATMOSPHERIC_OZONE[][$PRESSURE] - Geolocation - - - Integrated Ozone - Integrated Ozone - mdr[].INTEGRATED_OZONE[] - Geolocation - - - Surface Temperature 1 - Surface Temperature 1 - mdr[].SURFACE_TEMPERATURE[][0] - Geolocation - - - Surface Temperature 2 - Surface Temperature 2 - mdr[].SURFACE_TEMPERATURE[][1] - Geolocation - - - Integrated N20 - Integrated N20 - mdr[].INEGRATED_N2O[] - Geolocation - - - Integrated CO - Integrated CO - mdr[].INTEGRATED_CO[] - Geolocation - - - Integrated CH4 - Integrated CH4 - mdr[].INTEGRATED_CH4[] - Geolocation - - - Integrated CO2 - Integrated CO2 - mdr[].INTEGRATED_CO2[] - Geolocation - - - - WAVELENGTH - - - Surface Emissivity - Surface Emissivity - mdr[].SURFACE_EMISSIVITY[][$WAVELENGTH] - Geolocation - - - Number of cloud formations - Number of cloud formations - mdr[].NUMBER_CLOUD_FORMATIONS[] - Geolocation - - - - CLOUD - - - Fractional cloud cover - Fractional cloud cover - mdr[].FRACTIONAL_CLOUD_COVER[][$CLOUD] - Geolocation - - - - CLOUD - - - Cloud top temperature - Cloud top temperature - mdr[].CLOUD_TOP_TEMPERATURE[][$CLOUD] - Geolocation - - - - CLOUD - - - Cloud top pressure - Cloud top pressure - mdr[].CLOUD_TOP_PRESSURE[][$CLOUD] - Geolocation - - - - CLOUD - - - Cloud phase - Cloud phase (0=no cloud, 1=liquid, 2=ice, 3=mixed) - mdr[].CLOUD_PHASE[][$CLOUD] - Geolocation - - - FLG_AVHAVL - Completness flag for AVHRR products - mdr[].FLG_AVHAVL[] - Geolocation - - - FLG_CHNSEL - Channel selection flag - mdr[].FLG_CHNSEL[] - Geolocation - - - FLG_CLDAVH - Cloud formations anaylsed in AVHRR flag - mdr[].FLG_CLDAVH[] - Geolocation - - - FLG_CLDPHA - Cloud phase flag - mdr[].FLG_CLDPHA[] - Geolocation - - - FLG_DAYNIT - Flag discriminating between day and night - mdr[].FLG_DAYNIT[] - Geolocation - - - FLG_IASICLR - IASI-l2 IFOV clear, partly cloudy or cloudy flag - mdr[].FLG_IASICLR[] - Geolocation - - - FLG_ITCONV - Convergence of the iterative retrieval flag - mdr[].FLG_ITCONV[] - Geolocation - - - FLG_ITRBOU - Validation of iterated state vector flag - mdr[].FLG_ITRBOU[] - Geolocation - - - FLG_LANSEA - Flag specifying surface type - mdr[].FLG_LANSEA[] - Geolocation - - - FLG_NUMIT - Number of iterations used for retrieval - mdr[].FLG_NUMIT[] - Geolocation - - - FLG_NWPBAD - Validation flag of NWP forecast - mdr[].FLG_NWPBAD[] - Geolocation - - - FLG_QUAL - Quality and completness of the retrieval - mdr[].FLG_QUAL[] - Geolocation - - - FLG_RESID - Retrieval residual acceptance flag - mdr[].FLG_RESID[] - Geolocation - - - FLG_SATMAN - Flag indicating satellite manoeuvre - mdr[].FLG_SATMAN[] - Geolocation - - - FLG_SELBAC - Selection of background state flag - mdr[].FLG_SELBAC[] - Geolocation - - - FLG_SUNGLNT - Flag identifying of sun glint - mdr[].FLG_SUNGLNT[] - Geolocation - - - FLG_SUPADI - Flag indicating super-adiabatic conditions in final retrieval - mdr[].FLG_SUPADI[] - Geolocation - - - FLG_SUPSAT - Flag indicating super-saturation in final retrieval - mdr[].FLG_SUPSAT[] - Geolocation - - - FLG_THICIR - Thin cirrus cloud test - mdr[].FLG_THICIR[] - Geolocation - - - FLG_THICOR - Thin cirrus has been corrected for - mdr[].FLG_THICOR[] - Geolocation - - - FLG_VARCLR - Cloud clearing by variational analysis - mdr[].FLG_VARCLR[] - Geolocation - - eps-product - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/satpy/readers/eps_base.py b/satpy/readers/eps_base.py deleted file mode 100644 index 6ffcaaefdb..0000000000 --- a/satpy/readers/eps_base.py +++ /dev/null @@ -1,209 +0,0 @@ -# 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 . - -"""Reader for eps level 1b data. Uses xml files as a format description.""" - -import logging - -import dask.array as da -import numpy as np -import xarray as xr -from pyresample.geometry import SwathDefinition - -from satpy._config import get_config_path -from satpy.readers.file_handlers import BaseFileHandler -from satpy.readers.xmlformat import XMLFormat -from satpy.utils import get_legacy_chunk_size - -logger = logging.getLogger(__name__) - -CHUNK_SIZE = get_legacy_chunk_size() - -C1 = 1.191062e-05 # mW/(m2*sr*cm-4) -C2 = 1.4387863 # K/cm-1 - - -record_class = ["Reserved", "mphr", "sphr", - "ipr", "geadr", "giadr", - "veadr", "viadr", "mdr"] - - -def read_records(filename, format_definition): - """Read ``filename`` without scaling it afterwards. - - Args: - filename (str): Filename to read. - format_definition (str): XML file containing the format definition. - Must be findable in the Satpy config path. - """ - format_fn = get_config_path(format_definition) - form = XMLFormat(format_fn) - - grh_dtype = np.dtype([("record_class", "|i1"), - ("INSTRUMENT_GROUP", "|i1"), - ("RECORD_SUBCLASS", "|i1"), - ("RECORD_SUBCLASS_VERSION", "|i1"), - ("RECORD_SIZE", ">u4"), - ("RECORD_START_TIME", "S6"), - ("RECORD_STOP_TIME", "S6")]) - - max_lines = np.floor((CHUNK_SIZE ** 2) / 2048) - - dtypes = [] - cnt = 0 - counts = [] - classes = [] - prev = None - with open(filename, "rb") as fdes: - while True: - grh = np.fromfile(fdes, grh_dtype, 1) - if grh.size == 0: - break - rec_class = record_class[int(grh["record_class"])] - sub_class = grh["RECORD_SUBCLASS"][0] - - expected_size = int(grh["RECORD_SIZE"]) - bare_size = expected_size - grh_dtype.itemsize - try: - the_type = form.dtype((rec_class, sub_class)) - # the_descr = grh_dtype.descr + the_type.descr - except KeyError: - the_type = np.dtype([('unknown', 'V%d' % bare_size)]) - the_descr = grh_dtype.descr + the_type.descr - the_type = np.dtype(the_descr) - if the_type.itemsize < expected_size: - padding = [('unknown%d' % cnt, 'V%d' % (expected_size - the_type.itemsize))] - cnt += 1 - the_descr += padding - new_dtype = np.dtype(the_descr) - key = (rec_class, sub_class) - if key == prev: - counts[-1] += 1 - else: - dtypes.append(new_dtype) - counts.append(1) - classes.append(key) - prev = key - fdes.seek(expected_size - grh_dtype.itemsize, 1) - - sections = {} - offset = 0 - for dtype, count, rec_class in zip(dtypes, counts, classes): - fdes.seek(offset) - if rec_class == ('mdr', 2): - record = da.from_array(np.memmap(fdes, mode='r', dtype=dtype, shape=count, offset=offset), - chunks=(max_lines,)) - else: - record = np.fromfile(fdes, dtype=dtype, count=count) - offset += dtype.itemsize * count - if rec_class in sections: - logger.debug('Multiple records for ', str(rec_class)) - sections[rec_class] = np.hstack((sections[rec_class], record)) - else: - sections[rec_class] = record - - return sections, form - - -def create_xarray(arr, dims=("y", "x"), attrs=None): - """Create xarray with correct dimensions.""" - res = arr - if attrs is None: - attrs = {} - res = xr.DataArray(res, dims=dims, attrs=attrs) - return res - - -class EPSBaseFileHandler(BaseFileHandler): - """Base class for EPS level 1b readers.""" - - spacecrafts = {"M01": "Metop-B", - "M02": "Metop-A", - "M03": "Metop-C", } - - xml_conf: str - mdr_subclass: int - - def __init__(self, filename, filename_info, filetype_info): - """Initialize FileHandler.""" - super().__init__(filename, filename_info, filetype_info) - - self.area = None - self._start_time = filename_info['start_time'] - self._end_time = filename_info['end_time'] - self.form = None - self.scanlines = None - self.sections = None - - def _read_all(self): - logger.debug("Reading %s", self.filename) - self.sections, self.form = read_records(self.filename, self.xml_conf) - self.scanlines = self['TOTAL_MDR'] - if self.scanlines != len(self.sections[('mdr', self.mdr_subclass)]): - logger.warning("Number of declared records doesn't match number of scanlines in the file.") - self.scanlines = len(self.sections[('mdr', self.mdr_subclass)]) - - def __getitem__(self, key): - """Get value for given key.""" - for altkey in self.form.scales: - try: - try: - return self.sections[altkey][key] * self.form.scales[altkey][key] - except TypeError: - val = self.sections[altkey][key].item().decode().split("=")[1] - try: - return float(val) * self.form.scales[altkey][key].item() - except ValueError: - return val.strip() - except (KeyError, ValueError): - continue - raise KeyError("No matching value for " + str(key)) - - def keys(self): - """List of reader's keys.""" - keys = [] - for val in self.form.scales.values(): - keys += val.dtype.fields.keys() - return keys - - def get_lonlats(self): - """Get lonlats.""" - if self.area is None: - lons, lats = self.get_full_lonlats() - self.area = SwathDefinition(lons, lats) - self.area.name = '_'.join([self.platform_name, str(self.start_time), - str(self.end_time)]) - return self.area - - @property - def platform_name(self): - """Get platform name.""" - return self.spacecrafts[self["SPACECRAFT_ID"]] - - @property - def sensor_name(self): - """Get sensor name.""" - return self.sensors[self["INSTRUMENT_ID"]] - - @property - def start_time(self): - """Get start time.""" - return self._start_time - - @property - def end_time(self): - """Get end time.""" - return self._end_time diff --git a/satpy/readers/eps_l1b.py b/satpy/readers/eps_l1b.py index 23a6071d7a..1cc098a612 100644 --- a/satpy/readers/eps_l1b.py +++ b/satpy/readers/eps_l1b.py @@ -1,4 +1,6 @@ -# Copyright (c) 2017-2023 Satpy developers +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2017-2020 Satpy developers # # This file is part of satpy. # @@ -14,7 +16,7 @@ # You should have received a copy of the GNU General Public License along with # satpy. If not, see . -"""Reader for EPS level 1b data. Uses xml files as a format description.""" +"""Reader for eps level 1b data. Uses xml files as a format description.""" import functools import logging @@ -23,9 +25,12 @@ import numpy as np import xarray as xr from dask.delayed import delayed +from pyresample.geometry import SwathDefinition from satpy._compat import cached_property -from satpy.readers.eps_base import EPSBaseFileHandler, XMLFormat, read_records, record_class # noqa +from satpy._config import get_config_path +from satpy.readers.file_handlers import BaseFileHandler +from satpy.readers.xmlformat import XMLFormat from satpy.utils import get_legacy_chunk_size logger = logging.getLogger(__name__) @@ -46,6 +51,82 @@ def radiance_to_refl(arr, solar_flux): return arr * np.pi * 100.0 / solar_flux +record_class = ["Reserved", "mphr", "sphr", + "ipr", "geadr", "giadr", + "veadr", "viadr", "mdr"] + + +def read_records(filename): + """Read *filename* without scaling it afterwards.""" + format_fn = get_config_path("eps_avhrrl1b_6.5.xml") + form = XMLFormat(format_fn) + + grh_dtype = np.dtype([("record_class", "|i1"), + ("INSTRUMENT_GROUP", "|i1"), + ("RECORD_SUBCLASS", "|i1"), + ("RECORD_SUBCLASS_VERSION", "|i1"), + ("RECORD_SIZE", ">u4"), + ("RECORD_START_TIME", "S6"), + ("RECORD_STOP_TIME", "S6")]) + + max_lines = np.floor((CHUNK_SIZE ** 2) / 2048) + + dtypes = [] + cnt = 0 + counts = [] + classes = [] + prev = None + with open(filename, "rb") as fdes: + while True: + grh = np.fromfile(fdes, grh_dtype, 1) + if grh.size == 0: + break + rec_class = record_class[int(grh["record_class"])] + sub_class = grh["RECORD_SUBCLASS"][0] + + expected_size = int(grh["RECORD_SIZE"]) + bare_size = expected_size - grh_dtype.itemsize + try: + the_type = form.dtype((rec_class, sub_class)) + # the_descr = grh_dtype.descr + the_type.descr + except KeyError: + the_type = np.dtype([('unknown', 'V%d' % bare_size)]) + the_descr = grh_dtype.descr + the_type.descr + the_type = np.dtype(the_descr) + if the_type.itemsize < expected_size: + padding = [('unknown%d' % cnt, 'V%d' % (expected_size - the_type.itemsize))] + cnt += 1 + the_descr += padding + new_dtype = np.dtype(the_descr) + key = (rec_class, sub_class) + if key == prev: + counts[-1] += 1 + else: + dtypes.append(new_dtype) + counts.append(1) + classes.append(key) + prev = key + fdes.seek(expected_size - grh_dtype.itemsize, 1) + + sections = {} + offset = 0 + for dtype, count, rec_class in zip(dtypes, counts, classes): + fdes.seek(offset) + if rec_class == ('mdr', 2): + record = da.from_array(np.memmap(fdes, mode='r', dtype=dtype, shape=count, offset=offset), + chunks=(max_lines,)) + else: + record = np.fromfile(fdes, dtype=dtype, count=count) + offset += dtype.itemsize * count + if rec_class in sections: + logger.debug('Multiple records for ', str(rec_class)) + sections[rec_class] = np.hstack((sections[rec_class], record)) + else: + sections[rec_class] = record + + return sections, form + + def create_xarray(arr): """Create xarray with correct dimensions.""" res = arr @@ -53,33 +134,69 @@ def create_xarray(arr): return res -class EPSAVHRRFile(EPSBaseFileHandler): - """EPS level 1b reader for AVHRR data.""" +class EPSAVHRRFile(BaseFileHandler): + """Eps level 1b reader for AVHRR data.""" + + spacecrafts = {"M01": "Metop-B", + "M02": "Metop-A", + "M03": "Metop-C", } sensors = {"AVHR": "avhrr-3"} units = {"reflectance": "%", "brightness_temperature": "K"} - xml_conf = "eps_avhrrl1b_6.5.xml" - mdr_subclass = 2 - def __init__(self, filename, filename_info, filetype_info): """Initialize FileHandler.""" - super().__init__(filename, filename_info, filetype_info) - + super(EPSAVHRRFile, self).__init__( + filename, filename_info, filetype_info) + + self.area = None + self._start_time = filename_info['start_time'] + self._end_time = filename_info['end_time'] + self.form = None + self.scanlines = None + self.pixels = None + self.sections = None self.get_full_angles = functools.lru_cache(maxsize=1)( self._get_full_angles_uncached ) self.get_full_lonlats = functools.lru_cache(maxsize=1)( self._get_full_lonlats_uncached ) - self.pixels = None def _read_all(self): - super()._read_all() + logger.debug("Reading %s", self.filename) + self.sections, self.form = read_records(self.filename) + self.scanlines = self['TOTAL_MDR'] + if self.scanlines != len(self.sections[('mdr', 2)]): + logger.warning("Number of declared records doesn't match number of scanlines in the file.") + self.scanlines = len(self.sections[('mdr', 2)]) self.pixels = self["EARTH_VIEWS_PER_SCANLINE"] + def __getitem__(self, key): + """Get value for given key.""" + for altkey in self.form.scales: + try: + try: + return self.sections[altkey][key] * self.form.scales[altkey][key] + except TypeError: + val = self.sections[altkey][key].item().decode().split("=")[1] + try: + return float(val) * self.form.scales[altkey][key].item() + except ValueError: + return val.strip() + except (KeyError, ValueError): + continue + raise KeyError("No matching value for " + str(key)) + + def keys(self): + """List of reader's keys.""" + keys = [] + for val in self.form.scales.values(): + keys += val.dtype.fields.keys() + return keys + def _get_full_lonlats_uncached(self): """Get the interpolated longitudes and latitudes.""" raw_lats = np.hstack((self["EARTH_LOCATION_FIRST"][:, [0]], @@ -165,19 +282,6 @@ def get_bounding_box(self): self["EARTH_LOCATION_FIRST"][-1, [1]]]) return lons.ravel(), lats.ravel() - def _get_angle_dataarray(self, key): - """Get an angle dataarray.""" - sun_azi, sun_zen, sat_azi, sat_zen = self.get_full_angles() - if key['name'] == 'solar_zenith_angle': - return create_xarray(sun_zen) - if key['name'] == 'solar_azimuth_angle': - return create_xarray(sun_azi) - if key['name'] == 'satellite_zenith_angle': - return create_xarray(sat_zen) - if key['name'] == 'satellite_azimuth_angle': - return create_xarray(sat_azi) - raise ValueError(f"Unknown angle data-array: {key['name']:s}") - def get_dataset(self, key, info): """Get calibrated channel data.""" if self.sections is None: @@ -207,6 +311,19 @@ def get_dataset(self, key, info): dataset.attrs.update(key.to_dict()) return dataset + def _get_angle_dataarray(self, key): + """Get an angle dataarray.""" + sun_azi, sun_zen, sat_azi, sat_zen = self.get_full_angles() + if key['name'] == 'solar_zenith_angle': + dataset = create_xarray(sun_zen) + elif key['name'] == 'solar_azimuth_angle': + dataset = create_xarray(sun_azi) + if key['name'] == 'satellite_zenith_angle': + dataset = create_xarray(sat_zen) + elif key['name'] == 'satellite_azimuth_angle': + dataset = create_xarray(sat_azi) + return dataset + @cached_property def three_a_mask(self): """Mask for 3A.""" @@ -250,3 +367,34 @@ def _get_calibrated_dataarray(self, key): if mask is not None: dataset = dataset.where(~mask) return dataset + + def get_lonlats(self): + """Get lonlats.""" + if self.area is None: + lons, lats = self.get_full_lonlats() + self.area = SwathDefinition(lons, lats) + self.area.name = '_'.join([self.platform_name, str(self.start_time), + str(self.end_time)]) + return self.area + + @property + def platform_name(self): + """Get platform name.""" + return self.spacecrafts[self["SPACECRAFT_ID"]] + + @property + def sensor_name(self): + """Get sensor name.""" + return self.sensors[self["INSTRUMENT_ID"]] + + @property + def start_time(self): + """Get start time.""" + # return datetime.strptime(self["SENSING_START"], "%Y%m%d%H%M%SZ") + return self._start_time + + @property + def end_time(self): + """Get end time.""" + # return datetime.strptime(self["SENSING_END"], "%Y%m%d%H%M%SZ") + return self._end_time diff --git a/satpy/readers/eps_l2.py b/satpy/readers/eps_l2.py deleted file mode 100644 index 7a61c9ecb2..0000000000 --- a/satpy/readers/eps_l2.py +++ /dev/null @@ -1,219 +0,0 @@ -# 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 . - -"""Reader for EPS level 2 data. Uses xml files as a format description.""" - -import logging - -import numpy as np - -from satpy.readers.eps_base import EPSBaseFileHandler, create_xarray - -from . import epsnative_reader - -logger = logging.getLogger(__name__) - -var_2nd_dim = { - "co_cp_air": "co_nbr", - "co_cp_co_a": "co_nbr", - "co_x_co": "co_nbr", - "co_h_eigenvalues": "co_nbr", - "co_h_eigenvectors": "co_nbr", - "ozone_error": "nerr", - "temperature_error": "nerr", - "water_vapour_error": "nerr", -} -var_3rd_dim = { - "fg_atmospheric_temperature": "nlt", - "fg_atmospheric_water_vapour": "nlq", - "fg_atmospheric_ozone": "nlo", - "atmospheric_temperature": "nlt", - "atmospheric_water_vapour": "nlq", - "atmospheric_ozone": "nlo", -} - -NR_FOV = 120 - - -class EPSIASIL2FileHandler(EPSBaseFileHandler): - """EPS level 2 reader for IASI data. - - Reader for the IASI Level 2 combined sounding products in native format. - - Overview of the data including links to the product user guide, product format - specification, validation reports, and other documents, can be found at the - EUMETSAT Data Services at https://data.eumetsat.int/product/EO:EUM:DAT:METOP:IASSND02 - - """ - - xml_conf = "eps_iasil2_9.0.xml" - mdr_subclass = 1 - - def __init__(self, filename, filename_info, filetype_info): - """Initialise handler.""" - super().__init__(filename, filename_info, filetype_info) - self.descr = epsnative_reader.assemble_descriptor("IASISND02") - - def get_dataset(self, key, info): - """Get calibrated data.""" - dsname = key["name"].upper() - vals = self[dsname] - dims = get_var_dimensions(dsname, vals, self.dimensions) - # FIXME: get variable attributes - return create_xarray(vals, dims=list(dims.keys())) - - def available_datasets(self, configured_datasets=None): - """Get available datasets.""" - if self.sections is None: - self._read_all() - common = {"file_type": "iasi_l2_eps", "resolution": 12000} - for key in self.keys(): - yield (True, {"name": key.lower()} | common) - - def _read_all(self): - """Read all variables, lazily.""" - super()._read_all() - self.giadr = read_giadr(self.filename, self.descr) - self.dimensions = assemble_dimensions( - self.giadr, - self.scanlines, - self["NUM_ERROR_DATA"].max()) # Correct? - - -def read_giadr(input_product, descriptor, **kwargs): - """FIXME DOC. - - :param input_product: - :param descriptor: - :param kwargs: - :return: - """ - grh_size = 20 # bytes - ipr_sequence = epsnative_reader.read_ipr_sequence(input_product) - giadr_class = [cl for cl in ipr_sequence if "giadr" in cl["class"]][0] - with open(input_product, "rb") as eps_fileobj: - eps_fileobj.seek(giadr_class["offset"] + grh_size) - class_data = {} - giadr_descriptor = descriptor.get(("giadr", 1, 4)) - for _, row in giadr_descriptor.iterrows(): - # empty row or general header - if np.isnan(row["FIELD_SIZE"]) or row["FIELD"] == "RECORD_HEADER": - continue - else: - class_data[row["FIELD"].strip()] = { - "description": row["DESCRIPTION"], - "units": row["UNITS"], - "values": read_values(eps_fileobj, row, **kwargs), - } - return class_data - - -def read_values(eps_fileobj, row, reshape=False): - """FIXME DOC. - - :param _io.BinaryIO eps_fileobj: - :param pandas.core.series.Series row: - :param bool reshape: - :return numpy.ndarray: - """ - dtype = epsnative_reader.reckon_dtype(row["TYPE"].strip()) - shape = [int(row[k]) for k in row.keys() if k.lower().startswith("dim")] - shape = [k for k in shape if k > 1] - count = int(np.product(shape)) - values = np.fromfile(eps_fileobj, dtype, count) - # if np.isfinite(row['SF']): - # values = values / 10 ** row['SF'] - if row["TYPE"].startswith("v-"): - values = values["f1"] * 10.0 ** -values["f0"] - if row["TYPE"].strip().lower() == "short_cds_time": - values = [epsnative_reader.scds_to_datetime(*v) for v in values] - values = values[0] if len(values) == 1 else values - if reshape: - try: - values = values.reshape(tuple(reversed(shape))).T - except AttributeError: - values = values - return values - - -def assemble_dimensions(giadr, nr_rows, max_nerr): - """See Parameter Table (8.6) in EUM/OPS-EPS/MAN/04/0033. - - :param dict giadr: - :param int nr_rows: - :param int max_nerr: - :return dict: - """ - dimensions = { - "across_track": NR_FOV, - "along_track": nr_rows, - "cloud_formations": 3, - "co_nbr": NR_FOV, # pad to FOV number; real value for CO_NBR is provided in MDR - "nerr": max_nerr, # max value of NERR values provided in MDR - "nerro": giadr["NUM_OZONE_PCS"]["values"][0] - * (giadr["NUM_OZONE_PCS"]["values"][0] + 1) - / 2, - "nerrt": giadr["NUM_TEMPERATURE_PCS"]["values"][0] - * (giadr["NUM_TEMPERATURE_PCS"]["values"][0] + 1) - / 2, - "nerrw": giadr["NUM_WATER_VAPOUR_PCS"]["values"][0] - * (giadr["NUM_WATER_VAPOUR_PCS"]["values"][0] + 1) - / 2, - "neva_co": round(giadr["FORLI_NUM_LAYERS_CO"]["values"][0] / 2), - "neve_co": round(giadr["FORLI_NUM_LAYERS_CO"]["values"][0] / 2) - * giadr["FORLI_NUM_LAYERS_CO"]["values"][0], - "new": giadr["NUM_SURFACE_EMISSIVITY_WAVELENGTHS"]["values"][0], - "nl_co": giadr["FORLI_NUM_LAYERS_CO"]["values"][0], - "nl_hno3": giadr["FORLI_NUM_LAYERS_HNO3"]["values"][0], - "nl_o3": giadr["FORLI_NUM_LAYERS_O3"]["values"][0], - "nl_so2": giadr["BRESCIA_NUM_ALTITUDES_SO2"]["values"][0], - "nlo": giadr["NUM_PRESSURE_LEVELS_OZONE"]["values"][0], - "nlq": giadr["NUM_PRESSURE_LEVELS_HUMIDITY"]["values"][0], - "nlt": giadr["NUM_PRESSURE_LEVELS_TEMP"]["values"][0], - "npco": giadr["NUM_OZONE_PCS"]["values"][0], - "npct": giadr["NUM_TEMPERATURE_PCS"]["values"][0], - "npcw": giadr["NUM_WATER_VAPOUR_PCS"]["values"][0], - } - return dimensions - - -def get_var_dimensions(var_name, values, dimensions): - """FIXME DOC. - - :param str var_name: - :param numpy.ndarray values: - :param dict dimensions: - :return dict: - """ - # nlt, nlq, nlo are all equals to 101, so the automatic detection fails - dims_keys = list(dimensions.keys()) - sizes = list(dimensions.values()) - dims = {} - for idx, k in enumerate(values.shape): - if idx == 0: - dims["y"] = dimensions["along_track"] - elif idx == 1: - if var_name in var_2nd_dim: - dims[var_2nd_dim[var_name]] = k - else: - dims["x"] = k - else: - if var_name in var_3rd_dim: - dims[var_3rd_dim[var_name]] = k - else: - dim_name = dims_keys[sizes.index(k)] - dims[dim_name] = k - return dims diff --git a/satpy/readers/xmlformat.py b/satpy/readers/xmlformat.py index 4b14bbd7a7..0c46a3595e 100644 --- a/satpy/readers/xmlformat.py +++ b/satpy/readers/xmlformat.py @@ -19,26 +19,25 @@ from __future__ import annotations -from xml.etree.ElementTree import ElementTree # nosec +from xml.etree.ElementTree import ElementTree import numpy as np +VARIABLES: dict[str, str] = {} + TYPEC = {"boolean": ">i1", "integer2": ">i2", "integer4": ">i4", - "vinteger4": ">i4", - "uinteger1": ">u1", "uinteger2": ">u2", - "vuinteger2": ">u2", - "uinteger4": ">u4", - "vuinteger4": ">u4"} + "uinteger4": ">u4", } -def process_delimiter(elt, variables, ascii=False): +def process_delimiter(elt, ascii=False): """Process a 'delimiter' tag.""" + del elt, ascii -def process_field(elt, variables, ascii=False): +def process_field(elt, ascii=False): """Process a 'field' tag.""" # NOTE: if there is a variable defined in this field and it is different # from the default, we could change the value and restart. @@ -66,26 +65,28 @@ def process_field(elt, variables, ascii=False): return ((elt.get("name"), current_type, scale)) -def process_array(elt, variables, ascii=False): +def process_array(elt, ascii=False): """Process an 'array' tag.""" + del ascii chld = list(elt) if len(chld) > 1: raise ValueError() chld = chld[0] try: - name, current_type, scale = CASES[chld.tag](chld, variables, ascii) + name, current_type, scale = CASES[chld.tag](chld) size = None except ValueError: - name, current_type, size, scale = CASES[chld.tag](chld, variables, ascii) + name, current_type, size, scale = CASES[chld.tag](chld) + del name myname = elt.get("name") or elt.get("label") if elt.get("length").startswith("$"): - dimname = elt.get("length")[1:] - length = int(variables[dimname]) + length = int(VARIABLES[elt.get("length")[1:]]) else: length = int(elt.get("length")) if size is not None: return (myname, current_type, (length, ) + size, scale) - return (myname, current_type, (length, ), scale) + else: + return (myname, current_type, (length, ), scale) CASES = {"delimiter": process_delimiter, @@ -138,6 +139,36 @@ def to_scales(val): return scales +def parse_format(xml_file): + """Parse the xml file to create types, scaling factor types, and scales.""" + tree = ElementTree() + tree.parse(xml_file) + for param in tree.find("parameters"): + VARIABLES[param.get("name")] = param.get("value") + + types_scales = {} + + for prod in tree.find("product"): + ascii = (prod.tag in ["mphr", "sphr"]) + res = [] + for i in prod: + lres = CASES[i.tag](i, ascii) + if lres is not None: + res.append(lres) + types_scales[(prod.tag, int(prod.get("subclass")))] = res + + types = {} + stypes = {} + scales = {} + + for key, val in types_scales.items(): + types[key] = to_dtype(val) + stypes[key] = to_scaled_dtype(val) + scales[key] = to_scales(val) + + return types, stypes, scales + + def _apply_scales(array, scales, dtype): """Apply scales to the array.""" new_array = np.empty(array.shape, dtype) @@ -152,14 +183,12 @@ def _apply_scales(array, scales, dtype): return new_array -class XMLFormat: +class XMLFormat(object): """XMLFormat object.""" def __init__(self, filename): """Init the format reader.""" - self.dims = {} - self.variables = {} - self.types, self.stypes, self.scales = self.parse_format(filename) + self.types, self.stypes, self.scales = parse_format(filename) self.translator = {} @@ -174,31 +203,6 @@ def apply_scales(self, array): """Apply scales to *array*.""" return _apply_scales(array, *self.translator[array.dtype]) - def parse_format(self, xml_file): - """Parse the xml file to create types, scaling factor types, and scales.""" - tree = ElementTree() - tree.parse(xml_file) - for param in tree.find("parameters"): - self.variables[param.get("name")] = param.get("value") - - types_scales = {} - - for prod in tree.find("product"): - ascii = (prod.tag in ["mphr", "sphr"]) - res = [] - for i in prod: - lres = CASES[i.tag](i, self.variables, ascii) - if lres is not None: - res.append(lres) - types_scales[(prod.tag, int(prod.get("subclass")))] = res - - types = {} - stypes = {} - scales = {} - - for key, val in types_scales.items(): - types[key] = to_dtype(val) - stypes[key] = to_scaled_dtype(val) - scales[key] = to_scales(val) - - return types, stypes, scales + +if __name__ == '__main__': + pass diff --git a/satpy/tests/reader_tests/test_eps_l1b.py b/satpy/tests/reader_tests/test_eps_l1b.py index 4e3435ad23..3d085d5106 100644 --- a/satpy/tests/reader_tests/test_eps_l1b.py +++ b/satpy/tests/reader_tests/test_eps_l1b.py @@ -1,4 +1,6 @@ -# Copyright (c) 2019, 2022-2023 Satpy developers +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2019, 2022 Satpy developers # # This file is part of satpy. # @@ -13,7 +15,7 @@ # # You should have received a copy of the GNU General Public License along with # satpy. If not, see . -"""Test the EPS l1b format reader.""" +"""Test the eps l1b format.""" import os from contextlib import suppress diff --git a/satpy/tests/reader_tests/test_xmlformat.py b/satpy/tests/reader_tests/test_xmlformat.py deleted file mode 100644 index 1b22089ea5..0000000000 --- a/satpy/tests/reader_tests/test_xmlformat.py +++ /dev/null @@ -1,78 +0,0 @@ -# Copyright (c) 2018-2020 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 . -"""Tests for reading a format from an XML file.""" - -from xml.etree.ElementTree import Element - -import numpy as np -import pytest - - -@pytest.mark.parametrize("form", ["eps_avhrrl1b_6.5.xml", "eps_iasil2_9.0.xml"]) -def test_parse_format(form): - """Test parsing the XML format.""" - from satpy._config import get_config_path - from satpy.readers.xmlformat import XMLFormat - filename = get_config_path(form) - XMLFormat(filename) - - -def test_process_array(monkeypatch): - """Test processing an array tag.""" - from satpy.readers.xmlformat import process_array - elt = Element( - "array", - {"name": "SCENE_RADIANCES", - "length": "5", - "labels": "channel1,channel2,channel3,channel4,channel5"}) - elt2 = Element( - "array", - {"length": "$NE", "label": "FOV"}) - elt3 = Element( - "field", - {"type": "integer2", - "scaling-factor": "10^2,10^2,10^4,10^2,10^2"}) - elt.append(elt2) - elt2.append(elt3) - variables = {"NE": 10} -# dims = {} - (name, tp, shp, scl) = process_array(elt, variables, False) - assert name == "SCENE_RADIANCES" - assert tp == ">i2" - assert shp == (5, 10) - np.testing.assert_allclose( - scl, - np.array([0.01, 0.01, 0.0001, 0.01, 0.01])) -# assert dims == {"FOV": "NE"} - - elt = Element( - "array", - {"name": "ATMOSPHERIC_WATER_VAPOUR", - "length": "120"}) - elt2 = Element( - "array", - {"length": "$NLQ"}) - elt3 = Element( - "field", - {"type": "uinteger4", - "scaling-factor": "10^7", - "units": "kg/kg"}) - elt.append(elt2) - elt2.append(elt3) - variables = {"NLQ": 125} -# dims = {} - (name, tp, shp, scl) = process_array(elt, variables, False) -# assert dims == {"ATMOSPHERIC_WATER_VAPOUR": "NLQ"}