From 4164bb5ab5c3d234d068b016f01827a386544bb0 Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Tue, 21 Aug 2018 12:07:44 +0300 Subject: [PATCH 01/15] Convert IASI L2 reader to xarray/dask --- satpy/readers/iasi_l2.py | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/satpy/readers/iasi_l2.py b/satpy/readers/iasi_l2.py index 8c36afddc5..a76fdd9b06 100644 --- a/satpy/readers/iasi_l2.py +++ b/satpy/readers/iasi_l2.py @@ -24,11 +24,14 @@ import h5py import numpy as np +import xarray as xr +import xarray.ufuncs as xu +import dask.array as da import datetime as dt import logging -from satpy.dataset import Dataset from satpy.readers.file_handlers import BaseFileHandler +from satpy import CHUNK_SIZE # Scan timing values taken from # http://oiswww.eumetsat.org/WEBOPS/eps-pg/IASI-L1/IASIL1-PG-4ProdOverview.htm @@ -108,7 +111,7 @@ def get_dataset(self, key, info): m_data = read_dataset(fid, key, info) else: m_data = read_geo(fid, key, info) - m_data.info.update(info) + m_data.attrs.update(info) return m_data @@ -116,18 +119,25 @@ def get_dataset(self, key, info): def read_dataset(fid, key, info): """Read dataset""" dsid = DSET_NAMES[key.name] - data = fid["/PWLR/" + dsid].value + dset = fid["/PWLR/" + dsid] + if dset.ndim == 3: + dims = ['y', 'x', 'level'] + else: + dims = ['y', 'x'] + data = xr.DataArray(da.from_array(dset.value, chunks=CHUNK_SIZE), + name=key.name, dims=dims).astype(np.float32) try: - unit = fid["/PWLR/" + dsid].attrs['units'] - long_name = fid["/PWLR/" + dsid].attrs['long_name'] + units = dset.attrs['units'] + long_name = dset.attrs['long_name'] except KeyError: - unit = '' + units = '' long_name = '' + data.attrs['units'] = units + data.attrs['long_name'] = long_name - data = np.ma.masked_where(data > 1e30, data) + data = xr.where(data > 1e30, np.nan, data) - return Dataset(data, copy=False, long_name=long_name, - **info) + return data def read_geo(fid, key, info): @@ -136,13 +146,15 @@ def read_geo(fid, key, info): if "time" in key.name: days = fid["/L1C/" + dsid["day"]].value msecs = fid["/L1C/" + dsid["msec"]].value - unit = "" + units = "" data = _form_datetimes(days, msecs) else: data = fid["/L1C/" + dsid].value - unit = fid["/L1C/" + dsid].attrs['units'] + units = fid["/L1C/" + dsid].attrs['units'] - data = Dataset(data, copy=False, **info) + data = xr.DataArray(da.from_array(data, chunks=CHUNK_SIZE), + name=key.name, dims=['y', 'x']).astype(np.float32) + data.attrs['uints'] = units return data @@ -156,7 +168,7 @@ def _form_datetimes(days, msecs): day = int(days[i]) msec = msecs[i] scanline_datetimes = [] - for j in range(VALUES_PER_SCAN_LINE / 4): + for j in range(int(VALUES_PER_SCAN_LINE / 4)): usec = 1000 * (j * VIEW_TIME_ADJUSTMENT + msec) for k in range(4): delta = (dt.timedelta(days=day, microseconds=usec)) From ee22a07e157d5832aae19c6c8f7f16f3ef60572d Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Tue, 21 Aug 2018 12:27:24 +0300 Subject: [PATCH 02/15] Use seconds instead of datetime objects for the sensing time dataset --- satpy/readers/iasi_l2.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/satpy/readers/iasi_l2.py b/satpy/readers/iasi_l2.py index a76fdd9b06..9888db8189 100644 --- a/satpy/readers/iasi_l2.py +++ b/satpy/readers/iasi_l2.py @@ -143,25 +143,28 @@ def read_dataset(fid, key, info): def read_geo(fid, key, info): """Read geolocation and related datasets.""" dsid = GEO_NAMES[key.name] + add_epoch = False if "time" in key.name: days = fid["/L1C/" + dsid["day"]].value msecs = fid["/L1C/" + dsid["msec"]].value - units = "" data = _form_datetimes(days, msecs) + add_epoch = True else: data = fid["/L1C/" + dsid].value units = fid["/L1C/" + dsid].attrs['units'] data = xr.DataArray(da.from_array(data, chunks=CHUNK_SIZE), name=key.name, dims=['y', 'x']).astype(np.float32) - data.attrs['uints'] = units + # The units are taken from the reader configuration YAML file + # data.attrs['units'] = units + if add_epoch: + data.attrs['scan_time_epoch'] = EPOCH return data def _form_datetimes(days, msecs): - """Form datetimes from days and milliseconds relative to EPOCH for - each of IASI scans.""" + """Calculate seconds since EPOCH from days and milliseconds for each of IASI scan.""" all_datetimes = [] for i in range(days.size): @@ -172,7 +175,7 @@ def _form_datetimes(days, msecs): usec = 1000 * (j * VIEW_TIME_ADJUSTMENT + msec) for k in range(4): delta = (dt.timedelta(days=day, microseconds=usec)) - scanline_datetimes.append(EPOCH + delta) + scanline_datetimes.append(delta.total_seconds()) all_datetimes.append(np.array(scanline_datetimes)) return np.array(all_datetimes) From 4e60c1810c335c98f334082c5e311e2c17998fe0 Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Tue, 21 Aug 2018 13:10:53 +0300 Subject: [PATCH 03/15] Add rest of the datasets, handle attributes properly --- satpy/etc/readers/iasi_l2.yaml | 42 ++++++++++++++++++++++++++++++++++ satpy/readers/iasi_l2.py | 21 ++++++++++------- 2 files changed, 55 insertions(+), 8 deletions(-) diff --git a/satpy/etc/readers/iasi_l2.yaml b/satpy/etc/readers/iasi_l2.yaml index e86458b03f..b409262c37 100644 --- a/satpy/etc/readers/iasi_l2.yaml +++ b/satpy/etc/readers/iasi_l2.yaml @@ -33,6 +33,13 @@ datasets: resolution: 12000 coordinates: [longitude, latitude] + ozone_total_column: + name: ozone_total_column + file_type: iasi_l2_hdf5 + units: "kg/m^2" + resolution: 12000 + coordinates: [longitude, latitude] + pressure: name: pressure file_type: iasi_l2_hdf5 @@ -68,6 +75,41 @@ datasets: resolution: 12000 coordinates: [longitude, latitude] + water_total_column: + name: water_total_column + file_type: iasi_l2_hdf5 + units: "mm" + resolution: 12000 + coordinates: [longitude, latitude] + + surface_skin_temperature: + name: surface_skin_temperature + file_type: iasi_l2_hdf5 + units: "K" + resolution: 12000 + coordinates: [longitude, latitude] + + surface_skin_temperature_quality: + name: surface_skin_temperature_quality + file_type: iasi_l2_hdf5 + units: "" + resolution: 12000 + coordinates: [longitude, latitude] + + emissivity: + name: emissivity + file_type: iasi_l2_hdf5 + units: "1" + resolution: 12000 + coordinates: [longitude, latitude] + + emissivity_quality: + name: emissivity_quality + file_type: iasi_l2_hdf5 + units: "" + resolution: 12000 + coordinates: [longitude, latitude] + water_mixing_ratio_quality: name: water_mixing_ratio_quality file_type: iasi_l2_hdf5 diff --git a/satpy/readers/iasi_l2.py b/satpy/readers/iasi_l2.py index 9888db8189..887255369b 100644 --- a/satpy/readers/iasi_l2.py +++ b/satpy/readers/iasi_l2.py @@ -58,7 +58,14 @@ 'temperature': 'T', 'temperature_quality': 'QT', 'water_mixing_ratio': 'W', - 'water_mixing_ratio_quality': 'QW'} + 'water_mixing_ratio_quality': 'QW', + 'water_total_column': 'WC', + 'ozone_total_column': 'OC', + 'surface_skin_temperature': 'Ts', + 'surface_skin_temperature_quality': 'QTs', + 'emissivity': 'E', + 'emissivity_quality': 'EQ' +} GEO_NAMES = {'latitude': 'Latitude', 'longitude': 'Longitude', @@ -126,16 +133,14 @@ def read_dataset(fid, key, info): dims = ['y', 'x'] data = xr.DataArray(da.from_array(dset.value, chunks=CHUNK_SIZE), name=key.name, dims=dims).astype(np.float32) + data = xr.where(data > 1e30, np.nan, data) + try: - units = dset.attrs['units'] - long_name = dset.attrs['long_name'] + dset_attrs = dict(dset.attrs) except KeyError: - units = '' - long_name = '' - data.attrs['units'] = units - data.attrs['long_name'] = long_name + dset_attrs = {} - data = xr.where(data > 1e30, np.nan, data) + data.attrs.update(dset_attrs) return data From b0df3796836239c38e10f8b63bbc1b0fc03b9956 Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Tue, 21 Aug 2018 13:12:26 +0300 Subject: [PATCH 04/15] Fix typo in emissivity quality dataset name --- satpy/readers/iasi_l2.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/iasi_l2.py b/satpy/readers/iasi_l2.py index 887255369b..6f6fd6a754 100644 --- a/satpy/readers/iasi_l2.py +++ b/satpy/readers/iasi_l2.py @@ -64,7 +64,7 @@ 'surface_skin_temperature': 'Ts', 'surface_skin_temperature_quality': 'QTs', 'emissivity': 'E', - 'emissivity_quality': 'EQ' + 'emissivity_quality': 'QE' } GEO_NAMES = {'latitude': 'Latitude', From 5868d3281dfe44e65e2b1c5f3eafd674b379c3d5 Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Tue, 21 Aug 2018 13:13:31 +0300 Subject: [PATCH 05/15] Rename 'scan_time_epoch' to 'sensing_time_epoch' for consistency --- satpy/readers/iasi_l2.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/readers/iasi_l2.py b/satpy/readers/iasi_l2.py index 6f6fd6a754..19a6cde4b9 100644 --- a/satpy/readers/iasi_l2.py +++ b/satpy/readers/iasi_l2.py @@ -163,7 +163,7 @@ def read_geo(fid, key, info): # The units are taken from the reader configuration YAML file # data.attrs['units'] = units if add_epoch: - data.attrs['scan_time_epoch'] = EPOCH + data.attrs['sensing_time_epoch'] = EPOCH return data From 975239d1bf1c9ec1a8fde0ed0805f52411351c38 Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Tue, 21 Aug 2018 13:18:40 +0300 Subject: [PATCH 06/15] Remove unnecessary imports, arguments and comments --- satpy/readers/iasi_l2.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/satpy/readers/iasi_l2.py b/satpy/readers/iasi_l2.py index 19a6cde4b9..3c5c4125bd 100644 --- a/satpy/readers/iasi_l2.py +++ b/satpy/readers/iasi_l2.py @@ -25,7 +25,6 @@ import h5py import numpy as np import xarray as xr -import xarray.ufuncs as xu import dask.array as da import datetime as dt import logging @@ -115,15 +114,15 @@ def get_dataset(self, key, info): with h5py.File(self.filename, 'r') as fid: LOGGER.debug('Reading %s.', key.name) if key.name in DSET_NAMES: - m_data = read_dataset(fid, key, info) + m_data = read_dataset(fid, key) else: - m_data = read_geo(fid, key, info) + m_data = read_geo(fid, key) m_data.attrs.update(info) return m_data -def read_dataset(fid, key, info): +def read_dataset(fid, key): """Read dataset""" dsid = DSET_NAMES[key.name] dset = fid["/PWLR/" + dsid] @@ -145,7 +144,7 @@ def read_dataset(fid, key, info): return data -def read_geo(fid, key, info): +def read_geo(fid, key): """Read geolocation and related datasets.""" dsid = GEO_NAMES[key.name] add_epoch = False @@ -156,12 +155,9 @@ def read_geo(fid, key, info): add_epoch = True else: data = fid["/L1C/" + dsid].value - units = fid["/L1C/" + dsid].attrs['units'] data = xr.DataArray(da.from_array(data, chunks=CHUNK_SIZE), name=key.name, dims=['y', 'x']).astype(np.float32) - # The units are taken from the reader configuration YAML file - # data.attrs['units'] = units if add_epoch: data.attrs['sensing_time_epoch'] = EPOCH From cc15f2a7beafd34ee5353d90c7427e6924cd7093 Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Tue, 21 Aug 2018 13:24:06 +0300 Subject: [PATCH 07/15] Fix stickler complaint --- satpy/readers/iasi_l2.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/satpy/readers/iasi_l2.py b/satpy/readers/iasi_l2.py index 3c5c4125bd..4f2172b99e 100644 --- a/satpy/readers/iasi_l2.py +++ b/satpy/readers/iasi_l2.py @@ -63,8 +63,7 @@ 'surface_skin_temperature': 'Ts', 'surface_skin_temperature_quality': 'QTs', 'emissivity': 'E', - 'emissivity_quality': 'QE' -} + 'emissivity_quality': 'QE'} GEO_NAMES = {'latitude': 'Latitude', 'longitude': 'Longitude', From 5496f043d1f409789c6f9a2ec1f48695c323d06b Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Wed, 22 Aug 2018 11:19:06 +0300 Subject: [PATCH 08/15] Add unittests for IASI L2 reader --- satpy/tests/reader_tests/test_iasi_l2.py | 241 +++++++++++++++++++++++ 1 file changed, 241 insertions(+) create mode 100644 satpy/tests/reader_tests/test_iasi_l2.py diff --git a/satpy/tests/reader_tests/test_iasi_l2.py b/satpy/tests/reader_tests/test_iasi_l2.py new file mode 100644 index 0000000000..17b753dd99 --- /dev/null +++ b/satpy/tests/reader_tests/test_iasi_l2.py @@ -0,0 +1,241 @@ +#!/usr/bin/python +# Copyright (c) 2018. +# +# +# Author(s): +# Panu Lahtinen +# +# This program 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. +# +# This program 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 this program. If not, see . +# + +"""Unit tests for IASI L2 reader""" + +import os +import unittest + +import numpy as np + +SCAN_WIDTH = 120 +NUM_LEVELS = 138 +NUM_SCANLINES = 1 +FNAME = "W_XX-EUMETSAT-kan,iasi,metopb+kan_C_EUMS_20170920103559_IASI_PW3_02_M01_20170920102217Z_20170920102912Z.hdf" +# Structure for the test data, to be written to HDF5 file +TEST_DATA = { + # Not implemented in the reader + 'Amsu': { + 'FLG_AMSUBAD': {'data': np.zeros((NUM_SCANLINES, 30), dtype=np.uint8), + 'attrs': {} + } + }, +# Not implemented in the reader + 'INFO': { + 'OmC': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'long_name': "Cloud signal. Predicted average window channel 'Obs minus Calc", + 'units': 'K'} + }, + 'mdist': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {} + } + }, + 'L1C': { + 'Latitude': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'units': 'degrees_north'} + }, + 'Longitude': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'units': 'degrees_north'} + }, + 'SatAzimuth': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'units': 'degrees'} + }, + 'SatZenith': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'units': 'degrees'} + }, + 'SensingTime_day': {'data': np.array([6472], dtype=np.uint16), + 'attrs': {} + }, + 'SensingTime_msec': {'data': np.array([37337532], dtype=np.uint32), + 'attrs': {} + }, + 'SunAzimuth': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'units': 'degrees'} + }, + 'SunZenith': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'units': 'degrees'} + }, + }, + # Not implemented in the reader + 'Maps': { + 'Height': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'units': 'm'} + }, + 'HeightStd': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'units': 'm'} + }, + }, + # Not implemented in the reader + 'Mhs': { + 'FLG_MHSBAD': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.uint8), + 'attrs': {} + } + }, + 'PWLR': { + 'E': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH, 10), dtype=np.float32), + 'attrs': {'emissivity_wavenumbers': np.array([699.3, 826.4, + 925.9, 1075.2, + 1204.8, 1315.7, + 1724.1, 2000.0, + 2325.5, 2702.7], + dtype=np.float32)} + }, + 'O': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH, NUM_LEVELS), dtype=np.float32), + 'attrs': {'long_name': 'Ozone mixing ratio vertical profile', + 'units': 'kg/kg'} + }, + 'OC': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {} + }, + 'P': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH, NUM_LEVELS), dtype=np.float32), + 'attrs': {'long_name': 'Atmospheric pressures at which the vertical profiles are given. Last value is the surface pressure', + 'units': 'hpa'} + }, + 'QE': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {} + }, + 'QO': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {} + }, + 'QO': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {} + }, + 'QT': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {} + }, + 'QTs': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {} + }, + 'QW': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {} + }, + 'T': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH, NUM_LEVELS), dtype=np.float32), + 'attrs': {'long_name': 'Temperature vertical profile', + 'units': 'K'} + }, + 'Ts': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'long_name': 'Surface skin temperature', + 'units': 'K'} + }, + 'W': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH, NUM_LEVELS), dtype=np.float32), + 'attrs': {'long_name': 'Water vapour mixing ratio vertical profile', + 'units': 'kg/kg'} + }, + 'WC': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'attrs': {'long_name': 'Water vapour total columnar amount', + 'units': 'mm'} + }, + } +} + + +def save_test_data(path): + """Save the test to the indicated directory""" + import h5py + with h5py.File(os.path.join(path, FNAME), 'w') as fid: + # Create groups + for grp in TEST_DATA: + fid.create_group(grp) + # Write datasets + for dset in TEST_DATA[grp]: + fid[grp][dset] = TEST_DATA[grp][dset]['data'] + # Write dataset attributes + for attr in TEST_DATA[grp][dset]['attrs']: + fid[grp][dset].attrs[attr] = \ + TEST_DATA[grp][dset]['attrs'][attr] + + +class TestIasiL2(unittest.TestCase): + """Test IASI L2 reader.""" + + def setUp(self): + """Create temporary data to test on.""" + import tempfile + from datetime import datetime + + self.base_dir = tempfile.mkdtemp() + save_test_data(self.base_dir) + + def tearDown(self): + """Remove the temporary directory created for a test""" + try: + import shutil + shutil.rmtree(self.base_dir, ignore_errors=True) + except OSError: + pass + + def test_scene(self): + """Test scene creation""" + from satpy import Scene + fname = os.path.join(self.base_dir, FNAME) + scn = Scene(reader='iasi_l2', filenames=[fname,]) + self.assertTrue('start_time' in scn.attrs) + self.assertTrue('end_time' in scn.attrs) + self.assertTrue('sensor' in scn.attrs) + self.assertTrue('iasi' in scn.attrs['sensor']) + + def test_load_available_datasets(self): + """Test that all datasets are available""" + from satpy import Scene + fname = os.path.join(self.base_dir, FNAME) + scn = Scene(reader='iasi_l2', filenames=[fname,]) + scn.load(scn.available_dataset_names()) + + def test_load_pressure(self): + """Test loading pressure data""" + from satpy import Scene + fname = os.path.join(self.base_dir, FNAME) + scn = Scene(reader='iasi_l2', filenames=[fname,]) + scn.load(['pressure']) + pres = scn['pressure'].compute() + self.assertTrue(np.all(pres == 0.0)) + self.assertEqual(pres.x.size, SCAN_WIDTH) + self.assertEqual(pres.y.size, NUM_SCANLINES) + self.assertEqual(pres.level.size, NUM_LEVELS) + self.assertEqual(pres.attrs['start_time'], scn.attrs['start_time']) + self.assertEqual(pres.attrs['end_time'], scn.attrs['end_time']) + self.assertTrue('long_name' in pres.attrs) + self.assertTrue('units' in pres.attrs) + + def test_load_emissivity(self): + """Test loading emissivity data""" + from satpy import Scene + fname = os.path.join(self.base_dir, FNAME) + scn = Scene(reader='iasi_l2', filenames=[fname,]) + scn.load(['emissivity']) + emis = scn['emissivity'].compute() + self.assertTrue(np.all(emis == 0.0)) + self.assertEqual(emis.x.size, SCAN_WIDTH) + self.assertEqual(emis.y.size, NUM_SCANLINES) + self.assertTrue('emissivity_wavenumbers' in emis.attrs) + + +def suite(): + """The test suite for test_scene. + """ + loader = unittest.TestLoader() + mysuite = unittest.TestSuite() + mysuite.addTest(loader.loadTestsFromTestCase(TestIasiL2)) + return mysuite + + +if __name__ == '__main__': + unittest.main() From a45a2d3775545ac34120125620d364797bc5327c Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Wed, 22 Aug 2018 15:02:05 +0300 Subject: [PATCH 09/15] Use 64-bit floats for sensing times --- satpy/readers/iasi_l2.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/satpy/readers/iasi_l2.py b/satpy/readers/iasi_l2.py index 4f2172b99e..a4e6100a24 100644 --- a/satpy/readers/iasi_l2.py +++ b/satpy/readers/iasi_l2.py @@ -152,11 +152,13 @@ def read_geo(fid, key): msecs = fid["/L1C/" + dsid["msec"]].value data = _form_datetimes(days, msecs) add_epoch = True + dtype = np.float64 else: data = fid["/L1C/" + dsid].value - + dtype = np.float32 data = xr.DataArray(da.from_array(data, chunks=CHUNK_SIZE), - name=key.name, dims=['y', 'x']).astype(np.float32) + name=key.name, dims=['y', 'x']).astype(dtype) + if add_epoch: data.attrs['sensing_time_epoch'] = EPOCH @@ -173,9 +175,9 @@ def _form_datetimes(days, msecs): scanline_datetimes = [] for j in range(int(VALUES_PER_SCAN_LINE / 4)): usec = 1000 * (j * VIEW_TIME_ADJUSTMENT + msec) + delta = (dt.timedelta(days=day, microseconds=usec)) for k in range(4): - delta = (dt.timedelta(days=day, microseconds=usec)) scanline_datetimes.append(delta.total_seconds()) - all_datetimes.append(np.array(scanline_datetimes)) + all_datetimes.append(scanline_datetimes) - return np.array(all_datetimes) + return np.array(all_datetimes, dtype=np.float64) From 908a6b1c1dba94e43a875952b63eee0756c953bc Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Wed, 22 Aug 2018 15:03:52 +0300 Subject: [PATCH 10/15] Add unit test for sensing times --- satpy/tests/reader_tests/test_iasi_l2.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/satpy/tests/reader_tests/test_iasi_l2.py b/satpy/tests/reader_tests/test_iasi_l2.py index 17b753dd99..636bd952c0 100644 --- a/satpy/tests/reader_tests/test_iasi_l2.py +++ b/satpy/tests/reader_tests/test_iasi_l2.py @@ -227,6 +227,19 @@ def test_load_emissivity(self): self.assertEqual(emis.y.size, NUM_SCANLINES) self.assertTrue('emissivity_wavenumbers' in emis.attrs) + def test_load_sensing_times(self): + """Test loading sensing times""" + from satpy import Scene + fname = os.path.join(self.base_dir, FNAME) + scn = Scene(reader='iasi_l2', filenames=[fname,]) + scn.load(['sensing_time']) + times = scn['sensing_time'].compute() + # Times should be equal in blocks of four, but not beyond, so + # there should be SCAN_WIDTH/4 different values + for i in range(int(SCAN_WIDTH / 4)): + self.assertTrue(np.unique(times[0, i*4:i*4+4]).size == 1) + self.assertTrue(np.unique(times[0, :]).size == SCAN_WIDTH / 4) + def suite(): """The test suite for test_scene. From fa1c1714ca09460ae8be4577ca373c702f804e12 Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Thu, 23 Aug 2018 10:51:30 +0300 Subject: [PATCH 11/15] Add tests for bare reader without Scene() and module functions --- satpy/tests/reader_tests/test_iasi_l2.py | 129 +++++++++++++++++++---- 1 file changed, 110 insertions(+), 19 deletions(-) diff --git a/satpy/tests/reader_tests/test_iasi_l2.py b/satpy/tests/reader_tests/test_iasi_l2.py index 636bd952c0..46642af6ee 100644 --- a/satpy/tests/reader_tests/test_iasi_l2.py +++ b/satpy/tests/reader_tests/test_iasi_l2.py @@ -169,10 +169,23 @@ class TestIasiL2(unittest.TestCase): def setUp(self): """Create temporary data to test on.""" import tempfile - from datetime import datetime + import datetime as dt + from satpy.readers.iasi_l2 import IASIL2HDF5 self.base_dir = tempfile.mkdtemp() save_test_data(self.base_dir) + self.fname = os.path.join(self.base_dir, FNAME) + self.fname_info = {'start_time': dt.datetime(2017, 9, 20, 10, 22, 17), + 'end_time': dt.datetime(2017, 9, 20, 10, 29, 12), + 'processing_time': dt.datetime(2017, 9, 20, 10, 35, 59), + 'processing_location': 'kan', + 'long_platform_id': 'metopb', + 'instrument': 'iasi', + 'platform_id': 'M01'} + self.ftype_info = {'file_reader': IASIL2HDF5, + 'file_patterns': ['{fname}.hdf'], + 'file_type': 'iasi_l2_hdf5'} + self.reader = IASIL2HDF5(self.fname, self.fname_info, self.ftype_info) def tearDown(self): """Remove the temporary directory created for a test""" @@ -192,53 +205,131 @@ def test_scene(self): self.assertTrue('sensor' in scn.attrs) self.assertTrue('iasi' in scn.attrs['sensor']) - def test_load_available_datasets(self): + def test_scene_load_available_datasets(self): """Test that all datasets are available""" from satpy import Scene fname = os.path.join(self.base_dir, FNAME) scn = Scene(reader='iasi_l2', filenames=[fname,]) scn.load(scn.available_dataset_names()) - def test_load_pressure(self): + def test_scene_load_pressure(self): """Test loading pressure data""" from satpy import Scene fname = os.path.join(self.base_dir, FNAME) scn = Scene(reader='iasi_l2', filenames=[fname,]) scn.load(['pressure']) pres = scn['pressure'].compute() - self.assertTrue(np.all(pres == 0.0)) - self.assertEqual(pres.x.size, SCAN_WIDTH) - self.assertEqual(pres.y.size, NUM_SCANLINES) - self.assertEqual(pres.level.size, NUM_LEVELS) - self.assertEqual(pres.attrs['start_time'], scn.attrs['start_time']) - self.assertEqual(pres.attrs['end_time'], scn.attrs['end_time']) - self.assertTrue('long_name' in pres.attrs) - self.assertTrue('units' in pres.attrs) + self.check_pressure(pres, scn.attrs) - def test_load_emissivity(self): + def test_scene_load_emissivity(self): """Test loading emissivity data""" from satpy import Scene fname = os.path.join(self.base_dir, FNAME) scn = Scene(reader='iasi_l2', filenames=[fname,]) scn.load(['emissivity']) emis = scn['emissivity'].compute() - self.assertTrue(np.all(emis == 0.0)) - self.assertEqual(emis.x.size, SCAN_WIDTH) - self.assertEqual(emis.y.size, NUM_SCANLINES) - self.assertTrue('emissivity_wavenumbers' in emis.attrs) + self.check_emissivity(emis) - def test_load_sensing_times(self): + def test_scene_load_sensing_times(self): """Test loading sensing times""" from satpy import Scene fname = os.path.join(self.base_dir, FNAME) scn = Scene(reader='iasi_l2', filenames=[fname,]) scn.load(['sensing_time']) times = scn['sensing_time'].compute() + self.check_sensing_times(times) + + def test_init(self): + """Test reader initialization""" + self.assertEqual(self.reader.filename, self.fname) + self.assertEqual(self.reader.finfo, self.fname_info) + self.assertTrue(self.reader.lons is None) + self.assertTrue(self.reader.lats is None) + self.assertEqual(self.reader.mda['platform_name'], 'Metop-B') + self.assertEqual(self.reader.mda['sensor'], 'iasi') + + def test_time_properties(self): + """Test time properties""" + import datetime as dt + self.assertTrue(isinstance(self.reader.start_time, dt.datetime)) + self.assertTrue(isinstance(self.reader.end_time, dt.datetime)) + + def test_get_dataset(self): + """Test get_dataset() for different datasets""" + from satpy import DatasetID + info = {'eggs': 'spam'} + key = DatasetID(name='pressure') + data = self.reader.get_dataset(key, info).compute() + self.check_pressure(data) + self.assertTrue('eggs' in data.attrs) + self.assertEqual(data.attrs['eggs'], 'spam') + key = DatasetID(name='emissivity') + data = self.reader.get_dataset(key, info).compute() + self.check_emissivity(data) + key = DatasetID(name='sensing_time') + data = self.reader.get_dataset(key, info).compute() + self.assertEqual(data.shape, (NUM_SCANLINES, SCAN_WIDTH)) + + def check_pressure(self, pres, attrs=None): + """Helper method for testing reading pressure dataset""" + self.assertTrue(np.all(pres == 0.0)) + self.assertEqual(pres.x.size, SCAN_WIDTH) + self.assertEqual(pres.y.size, NUM_SCANLINES) + self.assertEqual(pres.level.size, NUM_LEVELS) + if attrs: + self.assertEqual(pres.attrs['start_time'], attrs['start_time']) + self.assertEqual(pres.attrs['end_time'], attrs['end_time']) + self.assertTrue('long_name' in pres.attrs) + self.assertTrue('units' in pres.attrs) + + def check_emissivity(self, emis): + """Helper method for testing reading emissivity dataset.""" + self.assertTrue(np.all(emis == 0.0)) + self.assertEqual(emis.x.size, SCAN_WIDTH) + self.assertEqual(emis.y.size, NUM_SCANLINES) + self.assertTrue('emissivity_wavenumbers' in emis.attrs) + + def check_sensing_times(self, times): + """Helper method for testing reading sensing times""" # Times should be equal in blocks of four, but not beyond, so # there should be SCAN_WIDTH/4 different values for i in range(int(SCAN_WIDTH / 4)): - self.assertTrue(np.unique(times[0, i*4:i*4+4]).size == 1) - self.assertTrue(np.unique(times[0, :]).size == SCAN_WIDTH / 4) + self.assertEqual(np.unique(times[0, i*4:i*4+4]).size, 1) + self.assertEqual(np.unique(times[0, :]).size, SCAN_WIDTH / 4) + + def test_read_dataset(self): + """Test read_dataset() function""" + import h5py + from satpy.readers.iasi_l2 import read_dataset + from satpy import DatasetID + with h5py.File(self.fname, 'r') as fid: + key = DatasetID(name='pressure') + data = read_dataset(fid, key).compute() + self.check_pressure(data) + key = DatasetID(name='emissivity') + data = read_dataset(fid, key).compute() + self.check_emissivity(data) + + def test_read_geo(self): + """Test read_geo() function""" + import h5py + from satpy.readers.iasi_l2 import read_geo + from satpy import DatasetID + with h5py.File(self.fname, 'r') as fid: + key = DatasetID(name='sensing_time') + data = read_geo(fid, key).compute() + self.assertEqual(data.shape, (NUM_SCANLINES, SCAN_WIDTH)) + key = DatasetID(name='latitude') + data = read_geo(fid, key).compute() + self.assertEqual(data.shape, (NUM_SCANLINES, SCAN_WIDTH)) + + def test_form_datetimes(self): + """Test _form_datetimes() function""" + from satpy.readers.iasi_l2 import _form_datetimes + days = TEST_DATA['L1C']['SensingTime_day']['data'] + msecs = TEST_DATA['L1C']['SensingTime_msec']['data'] + times = _form_datetimes(days, msecs) + self.check_sensing_times(times) def suite(): From 4da3ad8888e5b3d7bb3eb1599c720dd78c1fd44a Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Thu, 23 Aug 2018 12:11:15 +0300 Subject: [PATCH 12/15] Remove unnecessary try/except --- satpy/readers/iasi_l2.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/satpy/readers/iasi_l2.py b/satpy/readers/iasi_l2.py index a4e6100a24..6569e08424 100644 --- a/satpy/readers/iasi_l2.py +++ b/satpy/readers/iasi_l2.py @@ -133,11 +133,7 @@ def read_dataset(fid, key): name=key.name, dims=dims).astype(np.float32) data = xr.where(data > 1e30, np.nan, data) - try: - dset_attrs = dict(dset.attrs) - except KeyError: - dset_attrs = {} - + dset_attrs = dict(dset.attrs) data.attrs.update(dset_attrs) return data From 1c1d7dcf54b9d331c4d4f048085907a2371efc35 Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Thu, 23 Aug 2018 12:11:52 +0300 Subject: [PATCH 13/15] Test reading dataset with empty attributes --- satpy/tests/reader_tests/test_iasi_l2.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/satpy/tests/reader_tests/test_iasi_l2.py b/satpy/tests/reader_tests/test_iasi_l2.py index 46642af6ee..e704cb1194 100644 --- a/satpy/tests/reader_tests/test_iasi_l2.py +++ b/satpy/tests/reader_tests/test_iasi_l2.py @@ -309,6 +309,10 @@ def test_read_dataset(self): key = DatasetID(name='emissivity') data = read_dataset(fid, key).compute() self.check_emissivity(data) + # This dataset doesn't have any attributes + key = DatasetID(name='ozone_total_column') + data = read_dataset(fid, key).compute() + self.assertEqual(len(data.attrs), 0) def test_read_geo(self): """Test read_geo() function""" From 300fd9d37463bf764300a0ec8ab39243b1299594 Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Thu, 23 Aug 2018 12:12:22 +0300 Subject: [PATCH 14/15] Add test_iasi_l2 suite to reader test suite --- satpy/tests/reader_tests/__init__.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/satpy/tests/reader_tests/__init__.py b/satpy/tests/reader_tests/__init__.py index d1575c8548..9b7814a78a 100644 --- a/satpy/tests/reader_tests/__init__.py +++ b/satpy/tests/reader_tests/__init__.py @@ -32,7 +32,8 @@ test_acspo, test_amsr2_l1b, test_omps_edr, test_nucaps, test_geocat, test_seviri_calibration, test_clavrx, - test_grib, test_hrit_goes, test_ahi_hsd) + test_grib, test_hrit_goes, test_ahi_hsd, + test_iasi_l2) if sys.version_info < (2, 7): import unittest2 as unittest @@ -63,5 +64,6 @@ def suite(): mysuite.addTests(test_grib.suite()) mysuite.addTests(test_hrit_goes.suite()) mysuite.addTests(test_ahi_hsd.suite()) + mysuite.addTests(test_iasi_l2.suite()) return mysuite From c20fddfd82062de24219baafa4d246f61cd9b18f Mon Sep 17 00:00:00 2001 From: Panu Lahtinen Date: Fri, 24 Aug 2018 12:57:29 +0300 Subject: [PATCH 15/15] Fix duplicate dataset name --- satpy/tests/reader_tests/test_iasi_l2.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/satpy/tests/reader_tests/test_iasi_l2.py b/satpy/tests/reader_tests/test_iasi_l2.py index e704cb1194..128f1355e8 100644 --- a/satpy/tests/reader_tests/test_iasi_l2.py +++ b/satpy/tests/reader_tests/test_iasi_l2.py @@ -115,7 +115,7 @@ 'QO': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), 'attrs': {} }, - 'QO': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), + 'QP': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32), 'attrs': {} }, 'QT': {'data': np.zeros((NUM_SCANLINES, SCAN_WIDTH), dtype=np.float32),