From b6783295801199e0ea5626df873aeb02e7fc5298 Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Tue, 27 Aug 2019 16:56:52 +0200 Subject: [PATCH 01/16] Raise exception if colormap cannot be defined Signed-off-by: Adam Dybbroe --- satpy/composites/__init__.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/satpy/composites/__init__.py b/satpy/composites/__init__.py index adb6c15eb4..23fb557ca9 100644 --- a/satpy/composites/__init__.py +++ b/satpy/composites/__init__.py @@ -894,6 +894,9 @@ def build_colormap(palette, dtype, info): sf = info.get('scale_factor', np.array(1)) colormap.set_range( *info['valid_range'] * sf + info.get('add_offset', 0)) + else: + raise AttributeError("Data needs to have either a valid_range or be of type uint8" + + " in order to be displayable with an attached color-palette!") return colormap, sqpalette From 4fc0e9655d9ffa21e7f697df9fe60c49261850cc Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Tue, 27 Aug 2019 16:57:22 +0200 Subject: [PATCH 02/16] Handle hdf5 references Signed-off-by: Adam Dybbroe --- satpy/readers/hdf5_utils.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/satpy/readers/hdf5_utils.py b/satpy/readers/hdf5_utils.py index bbaf36ff59..e4b3b31c81 100644 --- a/satpy/readers/hdf5_utils.py +++ b/satpy/readers/hdf5_utils.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -# Copyright (c) 2016-2017 Satpy developers +# Copyright (c) 2016-2017, 2019 Satpy developers # # This file is part of satpy. # @@ -40,6 +40,7 @@ def __init__(self, filename, filename_info, filetype_info): super(HDF5FileHandler, self).__init__( filename, filename_info, filetype_info) self.file_content = {} + try: file_handle = h5py.File(self.filename, 'r') except IOError: @@ -57,7 +58,7 @@ def _collect_attrs(self, name, attrs): fc_key = "{}/attr/{}".format(name, key) try: self.file_content[fc_key] = np2str(value) - except ValueError: + except (ValueError, AttributeError): self.file_content[fc_key] = value def collect_metadata(self, name, obj): From 920146fb2a6029b794acab180f49dbd0adf77dea Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Tue, 27 Aug 2019 16:57:58 +0200 Subject: [PATCH 03/16] Add reader for the old NWCSAF/MSG (v2013) formatted cloud products Signed-off-by: Adam Dybbroe --- satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml | 175 +++++++++++++++++++++ satpy/readers/nwcsaf_msg2013_hdf5.py | 141 +++++++++++++++++ 2 files changed, 316 insertions(+) create mode 100644 satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml create mode 100644 satpy/readers/nwcsaf_msg2013_hdf5.py diff --git a/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml b/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml new file mode 100644 index 0000000000..e0fe25ad91 --- /dev/null +++ b/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml @@ -0,0 +1,175 @@ +reader: + description: HDF5 reader for the NWCSAF/Geo Seviri 2013 format + name: nwcsaf-msg2013-hdf5 + sensors: [seviri] + default_channels: [] + reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader + +file_types: + h5_nwcsaf_cma: + file_reader: !!python/name:satpy.readers.nwcsaf_msg2013_hdf5.Hdf5NWCSAF + # SAFNWC_MSG4_CMa__201908271145_MSG-N_______.PLAX.CTTH.0.h5 + file_patterns: ['SAFNWC_{platform_id}_CMa__{start_time:%Y%m%d%H%M}_{region_id:_<12s}.PLAX.CTTH.0.h5'] + + h5_nwcsaf_ct: + file_reader: !!python/name:satpy.readers.nwcsaf_msg2013_hdf5.Hdf5NWCSAF + # SAFNWC_MSG4_CT___201906241245_MSG-N_______.PLAX.CTTH.0.h5 + file_patterns: ['SAFNWC_{platform_id}_CT___{start_time:%Y%m%d%H%M}_{region_id:_<12s}.PLAX.CTTH.0.h5'] + + h5_nwcsaf_ctth: + file_reader: !!python/name:satpy.readers.nwcsaf_msg2013_hdf5.Hdf5NWCSAF + # SAFNWC_MSG4_CTTH_201906241245_MSG-N_______.PLAX.CTTH.0.h5 + file_patterns: ['SAFNWC_{platform_id}_CTTH_{start_time:%Y%m%d%H%M}_{region_id:_<12s}.PLAX.CTTH.0.h5'] + + +datasets: + +# ---- CMA products ------------ + + cma: + name: cma + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_cma + + cma_pal: + name: cma_pal + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_cma + + cma_cloudsnow: + name: cma_cloudsnow + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_cma + + cma_cloudsnow_pal: + name: cma_cloudsnow_pal + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_cma + + cma_dust: + name: cma_dust + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_cma + + cma_dust_pal: + name: cma_dust_pal + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_cma + + cma_volcanic: + name: cma_volcanic + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_cma + + cma_volcanic_pal: + name: cma_volcanic_pal + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_cma + + cma_conditions: + name: cma_conditions + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_cma + + cma_status_flag: + name: cma_status_flag + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_cma + +# ---- CT products ------------ + + ct: + name: ct + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ct + file_key: CT + + ct_pal: + name: ct_pal + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ct + file_key: 01-PALETTE + + ct_quality: + name: ct_quality + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ct + file_key: CT_QUALITY + + ct_phase: + name: ct_phase + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ct + file_key: CT_PHASE + + ct_phase_pal: + name: ct_phase_pal + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ct + file_key: 02-PALETTE + +# ---- CTTH products ------------ + + ctth_alti: + name: ctth_alti + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ctth + + ctth_alti_pal: + name: ctth_alti_pal + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ctth + + ctth_pres: + name: ctth_pres + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ctth + + ctth_pres_pal: + name: ctth_pres_pal + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ctth + + ctth_tempe: + name: ctth_tempe + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ctth + + ctth_tempe_pal: + name: ctth_tempe_pal + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ctth + + ctth_quality: + name: ctth_quality + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ctth + + ctth_status_flag: + name: ctth_status_flag + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ctth + diff --git a/satpy/readers/nwcsaf_msg2013_hdf5.py b/satpy/readers/nwcsaf_msg2013_hdf5.py new file mode 100644 index 0000000000..6f3b935ed5 --- /dev/null +++ b/satpy/readers/nwcsaf_msg2013_hdf5.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2019 Pytroll + +# Author(s): + +# Adam.Dybbroe + +# 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 . + +"""Reader for the old NWCSAF/Geo (v2013 and earlier) cloud product format + +References: + + - The NWCSAF GEO 2013 products documentation: + http://www.nwcsaf.org/web/guest/archive - Search for Code "ICD/3"; Type + "MSG" and the box to the right should say 'Status' (which means any + status). Version 7.0 seems to be for v2013 + + http://www.nwcsaf.org/aemetRest/downloadAttachment/2623 + + +""" + +import logging +import os +from datetime import datetime +import dask.array as da +import numpy as np +from pyresample.utils import get_area_def +from satpy import CHUNK_SIZE +from satpy.readers.hdf5_utils import HDF5FileHandler +from pyresample.geometry import AreaDefinition + + +logger = logging.getLogger(__name__) + +PLATFORM_NAMES = {'MSG1': 'Meteosat-8', + 'MSG2': 'Meteosat-9', + 'MSG3': 'Meteosat-10', + 'MSG4': 'Meteosat-11', } + + +class Hdf5NWCSAF(HDF5FileHandler): + """NWCSAF MSG hdf5 reader.""" + + def __init__(self, filename, filename_info, filetype_info): + """Init method.""" + super(Hdf5NWCSAF, self).__init__(filename, filename_info, + filetype_info) + + self.cache = {} + + def get_dataset(self, dataset_id, ds_info): + """Load a dataset.""" + file_key = ds_info.get('file_key', dataset_id.name) + data = self[file_key] + if 'SCALING_FACTOR' in data.attrs and 'OFFSET' in data.attrs: + dtype = np.dtype(data.data) + data.data = (data.data * data.attrs['SCALING_FACTOR'] + data.attrs['OFFSET']).astype(dtype) + + return data + + def get_area_def(self, dsid): + """Get the area definition of the datasets in the file. + + """ + + if dsid.name.endswith('_pal'): + raise NotImplementedError + + cfac = self.file_content['/attr/CFAC'] + lfac = self.file_content['/attr/LFAC'] + coff = self.file_content['/attr/COFF'] + loff = self.file_content['/attr/LOFF'] + numcols = int(self.file_content['/attr/NC']) + numlines = int(self.file_content['/attr/NL']) + + aex = get_area_extent(cfac, lfac, coff, loff, numcols, numlines) + pname = self.file_content['/attr/PROJECTION_NAME'] + proj = {} + if pname.startswith("GEOS"): + proj["proj"] = "geos" + proj["a"] = "6378169.0" + proj["b"] = "6356583.8" + proj["h"] = "35785831.0" + proj["lon_0"] = str(float(pname.split("<")[1][:-1])) + else: + raise NotImplementedError("Only geos projection supported yet.") + + area_def = AreaDefinition(self.file_content['/attr/REGION_NAME'], + self.file_content['/attr/REGION_NAME'], + pname, + proj, + numcols, + numlines, + aex) + + return area_def + + @property + def start_time(self): + """Return the start time of the object.""" + return datetime.strptime(self.file_content['/attr/IMAGE_ACQUISITION_TIME'], '%Y%m%d%H%M') + + @property + def end_time(self): + """Return the end time of the object.""" + return None + + +def get_area_extent(cfac, lfac, coff, loff, numcols, numlines): + """Get the area extent from msg parameters. + """ + + xur = (numcols - coff) * 2 ** 16 / (cfac * 1.0) + xur = np.deg2rad(xur) * 35785831.0 + xll = (-1 - coff) * 2 ** 16 / (cfac * 1.0) + xll = np.deg2rad(xll) * 35785831.0 + xres = (xur - xll) / numcols + xur, xll = xur - xres / 2, xll + xres / 2 + yll = (numlines - loff) * 2 ** 16 / (-lfac * 1.0) + yll = np.deg2rad(yll) * 35785831.0 + yur = (-1 - loff) * 2 ** 16 / (-lfac * 1.0) + yur = np.deg2rad(yur) * 35785831.0 + yres = (yur - yll) / numlines + yll, yur = yll + yres / 2, yur - yres / 2 + + return xll, yll, xur, yur From 3323e040585c7b06d4cbe0eb3e2960349c6911dd Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Wed, 28 Aug 2019 14:35:07 +0200 Subject: [PATCH 04/16] Handle hdf5 references Signed-off-by: Adam Dybbroe --- satpy/readers/hdf5_utils.py | 18 +++++++++++++++++- satpy/readers/nwcsaf_msg2013_hdf5.py | 7 ++++++- 2 files changed, 23 insertions(+), 2 deletions(-) diff --git a/satpy/readers/hdf5_utils.py b/satpy/readers/hdf5_utils.py index e4b3b31c81..01f6faab55 100644 --- a/satpy/readers/hdf5_utils.py +++ b/satpy/readers/hdf5_utils.py @@ -58,8 +58,24 @@ def _collect_attrs(self, name, attrs): fc_key = "{}/attr/{}".format(name, key) try: self.file_content[fc_key] = np2str(value) - except (ValueError, AttributeError): + except ValueError: self.file_content[fc_key] = value + except AttributeError: + # A HDF5 reference ? + value = self.get_reference(name, key) + if value is None: + LOG.warning("Value cannot be converted - skip setting attribute %s", fc_key) + else: + self.file_content[fc_key] = value + + def get_reference(self, name, key): + + with h5py.File(self.filename, 'r') as hf: + if type(hf[name].attrs[key]) is h5py.h5r.Reference: + ref_name = h5py.h5r.get_name(hf[name].attrs[key], hf.id) + return hf[ref_name][()] + else: + return None def collect_metadata(self, name, obj): if isinstance(obj, h5py.Dataset): diff --git a/satpy/readers/nwcsaf_msg2013_hdf5.py b/satpy/readers/nwcsaf_msg2013_hdf5.py index 6f3b935ed5..73651fe480 100644 --- a/satpy/readers/nwcsaf_msg2013_hdf5.py +++ b/satpy/readers/nwcsaf_msg2013_hdf5.py @@ -43,7 +43,7 @@ from satpy import CHUNK_SIZE from satpy.readers.hdf5_utils import HDF5FileHandler from pyresample.geometry import AreaDefinition - +import h5py logger = logging.getLogger(__name__) @@ -71,6 +71,11 @@ def get_dataset(self, dataset_id, ds_info): dtype = np.dtype(data.data) data.data = (data.data * data.attrs['SCALING_FACTOR'] + data.attrs['OFFSET']).astype(dtype) + for key in list(data.attrs.keys()): + val = data.attrs[key] + if isinstance(val, h5py.h5r.Reference): + del data.attrs[key] + return data def get_area_def(self, dsid): From ad41b771287a3f2a03ebc5bed624962fc62529e2 Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Wed, 28 Aug 2019 14:36:35 +0200 Subject: [PATCH 05/16] Add support for CTTH image products Signed-off-by: Adam Dybbroe --- satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml | 23 ++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml b/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml index e0fe25ad91..efc93797cb 100644 --- a/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml +++ b/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml @@ -130,46 +130,61 @@ datasets: sensor: seviri resolution: 3000 file_type: h5_nwcsaf_ctth + file_key: CTTH_HEIGHT ctth_alti_pal: name: ctth_alti_pal sensor: seviri resolution: 3000 file_type: h5_nwcsaf_ctth + file_key: 02-PALETTE ctth_pres: name: ctth_pres sensor: seviri resolution: 3000 file_type: h5_nwcsaf_ctth + file_key: CTTH_PRESS ctth_pres_pal: name: ctth_pres_pal sensor: seviri resolution: 3000 file_type: h5_nwcsaf_ctth + file_key: 01-PALETTE ctth_tempe: name: ctth_tempe sensor: seviri resolution: 3000 file_type: h5_nwcsaf_ctth + file_key: CTTH_TEMPER ctth_tempe_pal: name: ctth_tempe_pal sensor: seviri resolution: 3000 file_type: h5_nwcsaf_ctth + file_key: 03-PALETTE - ctth_quality: - name: ctth_quality + ctth_effective_cloudiness: + name: ctth_eff + sensor: seviri + resolution: 3000 + file_type: h5_nwcsaf_ctth + file_key: CTTH_TEMPER + + ctth_effective_cloudiness_pal: + name: ctth_eff_pal sensor: seviri resolution: 3000 file_type: h5_nwcsaf_ctth + file_key: 04-PALETTE - ctth_status_flag: - name: ctth_status_flag + ctth_quality: + name: ctth_quality sensor: seviri resolution: 3000 file_type: h5_nwcsaf_ctth + file_key: CTTH_QUALITY From 7db1dff78b4785e9a5eec284d31264e88ac8b7e8 Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Fri, 25 Oct 2019 08:54:26 +0200 Subject: [PATCH 06/16] Fix scaling of CTTH fields Signed-off-by: Adam Dybbroe --- satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml | 4 ++-- satpy/etc/readers/nwcsaf-pps_nc.yaml | 1 - satpy/readers/nwcsaf_msg2013_hdf5.py | 24 ++++++++++++++++------ 3 files changed, 20 insertions(+), 9 deletions(-) diff --git a/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml b/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml index efc93797cb..fbcffa7276 100644 --- a/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml +++ b/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml @@ -168,11 +168,11 @@ datasets: file_key: 03-PALETTE ctth_effective_cloudiness: - name: ctth_eff + name: ctth_effective_cloudiness sensor: seviri resolution: 3000 file_type: h5_nwcsaf_ctth - file_key: CTTH_TEMPER + file_key: CTTH_EFFECT ctth_effective_cloudiness_pal: name: ctth_eff_pal diff --git a/satpy/etc/readers/nwcsaf-pps_nc.yaml b/satpy/etc/readers/nwcsaf-pps_nc.yaml index 55f40d5817..902b8c85ec 100644 --- a/satpy/etc/readers/nwcsaf-pps_nc.yaml +++ b/satpy/etc/readers/nwcsaf-pps_nc.yaml @@ -162,7 +162,6 @@ datasets: file_type: nc_nwcsaf_ctth coordinates: [lon, lat] - ctth_pres: name: ctth_pres file_type: nc_nwcsaf_ctth diff --git a/satpy/readers/nwcsaf_msg2013_hdf5.py b/satpy/readers/nwcsaf_msg2013_hdf5.py index 73651fe480..030cb0f247 100644 --- a/satpy/readers/nwcsaf_msg2013_hdf5.py +++ b/satpy/readers/nwcsaf_msg2013_hdf5.py @@ -67,9 +67,26 @@ def get_dataset(self, dataset_id, ds_info): """Load a dataset.""" file_key = ds_info.get('file_key', dataset_id.name) data = self[file_key] + + nodata = None if 'SCALING_FACTOR' in data.attrs and 'OFFSET' in data.attrs: dtype = np.dtype(data.data) - data.data = (data.data * data.attrs['SCALING_FACTOR'] + data.attrs['OFFSET']).astype(dtype) + if dataset_id.name in ['ctth_alti']: + data.attrs['valid_range'] = (0, 27000) + data.attrs['_FillValue'] = np.nan + #data.attrs['_FillValue'] = 65535 + + if dataset_id.name in ['ctth_alti', 'ctth_pres', 'ctth_tempe', 'ctth_effective_cloudiness']: + #dtype = np.dtype('uint16') + dtype = np.dtype('float32') + nodata = 255 + + scaled_data = (data.data * data.attrs['SCALING_FACTOR'] + data.attrs['OFFSET']).astype(dtype) + if nodata: + scaled_data = np.where(data.data == nodata, np.nan, scaled_data) + scaled_data = np.where(scaled_data < 0, np.nan, scaled_data) + + data.data = scaled_data for key in list(data.attrs.keys()): val = data.attrs[key] @@ -120,11 +137,6 @@ def start_time(self): """Return the start time of the object.""" return datetime.strptime(self.file_content['/attr/IMAGE_ACQUISITION_TIME'], '%Y%m%d%H%M') - @property - def end_time(self): - """Return the end time of the object.""" - return None - def get_area_extent(cfac, lfac, coff, loff, numcols, numlines): """Get the area extent from msg parameters. From 2b7775a60d8d3e487883ef5000fad4bafc231670 Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Wed, 30 Oct 2019 16:59:55 +0100 Subject: [PATCH 07/16] Fix for stickler Signed-off-by: Adam Dybbroe --- satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml | 2 +- satpy/readers/hdf5_utils.py | 4 +--- satpy/readers/nwcsaf_msg2013_hdf5.py | 6 ------ 3 files changed, 2 insertions(+), 10 deletions(-) diff --git a/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml b/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml index fbcffa7276..3c431ad335 100644 --- a/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml +++ b/satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml @@ -108,7 +108,7 @@ datasets: resolution: 3000 file_type: h5_nwcsaf_ct file_key: CT_QUALITY - + ct_phase: name: ct_phase sensor: seviri diff --git a/satpy/readers/hdf5_utils.py b/satpy/readers/hdf5_utils.py index 01f6faab55..45919fccb4 100644 --- a/satpy/readers/hdf5_utils.py +++ b/satpy/readers/hdf5_utils.py @@ -71,11 +71,9 @@ def _collect_attrs(self, name, attrs): def get_reference(self, name, key): with h5py.File(self.filename, 'r') as hf: - if type(hf[name].attrs[key]) is h5py.h5r.Reference: + if isinstance(hf[name].attrs[key], h5py.h5r.Reference): ref_name = h5py.h5r.get_name(hf[name].attrs[key], hf.id) return hf[ref_name][()] - else: - return None def collect_metadata(self, name, obj): if isinstance(obj, h5py.Dataset): diff --git a/satpy/readers/nwcsaf_msg2013_hdf5.py b/satpy/readers/nwcsaf_msg2013_hdf5.py index 030cb0f247..955c0dc00b 100644 --- a/satpy/readers/nwcsaf_msg2013_hdf5.py +++ b/satpy/readers/nwcsaf_msg2013_hdf5.py @@ -35,12 +35,8 @@ """ import logging -import os from datetime import datetime -import dask.array as da import numpy as np -from pyresample.utils import get_area_def -from satpy import CHUNK_SIZE from satpy.readers.hdf5_utils import HDF5FileHandler from pyresample.geometry import AreaDefinition import h5py @@ -74,10 +70,8 @@ def get_dataset(self, dataset_id, ds_info): if dataset_id.name in ['ctth_alti']: data.attrs['valid_range'] = (0, 27000) data.attrs['_FillValue'] = np.nan - #data.attrs['_FillValue'] = 65535 if dataset_id.name in ['ctth_alti', 'ctth_pres', 'ctth_tempe', 'ctth_effective_cloudiness']: - #dtype = np.dtype('uint16') dtype = np.dtype('float32') nodata = 255 From 9057e81a6c1dda785eeddcd7c73e1c8e8a3dbebd Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Mon, 4 Nov 2019 17:32:39 +0100 Subject: [PATCH 08/16] Set valid range and fill value for Cloudtype data Signed-off-by: Adam Dybbroe --- satpy/readers/nwcsaf_msg2013_hdf5.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/satpy/readers/nwcsaf_msg2013_hdf5.py b/satpy/readers/nwcsaf_msg2013_hdf5.py index 955c0dc00b..7f29527556 100644 --- a/satpy/readers/nwcsaf_msg2013_hdf5.py +++ b/satpy/readers/nwcsaf_msg2013_hdf5.py @@ -75,6 +75,11 @@ def get_dataset(self, dataset_id, ds_info): dtype = np.dtype('float32') nodata = 255 + if dataset_id.name in ['ct']: + data.attrs['valid_range'] = (0, 20) + data.attrs['_FillValue'] = 255 + # data.attrs['palette_meanings'] = list(range(21)) + scaled_data = (data.data * data.attrs['SCALING_FACTOR'] + data.attrs['OFFSET']).astype(dtype) if nodata: scaled_data = np.where(data.data == nodata, np.nan, scaled_data) From 0a6499b5123577e34ca1e6878f97935524a68b98 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Mon, 4 Nov 2019 17:52:33 +0100 Subject: [PATCH 09/16] Add test for references in hdf5_utils --- satpy/tests/reader_tests/test_hdf5_utils.py | 23 ++++++++++++--------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/satpy/tests/reader_tests/test_hdf5_utils.py b/satpy/tests/reader_tests/test_hdf5_utils.py index 46b1a0075f..115da55dbc 100644 --- a/satpy/tests/reader_tests/test_hdf5_utils.py +++ b/satpy/tests/reader_tests/test_hdf5_utils.py @@ -15,8 +15,7 @@ # # You should have received a copy of the GNU General Public License along with # satpy. If not, see . -"""Module for testing the satpy.readers.hdf5_utils module. -""" +"""Module for testing the satpy.readers.hdf5_utils module.""" import os import sys @@ -35,10 +34,10 @@ class FakeHDF5FileHandler(HDF5FileHandler): - """Swap-in HDF5 File Handler for reader tests to use""" + """Swap HDF5 File Handler for reader tests to use.""" def __init__(self, filename, filename_info, filetype_info, **kwargs): - """Get fake file content from 'get_test_content'""" + """Get fake file content from 'get_test_content'.""" if HDF5FileHandler is object: raise ImportError("Base 'HDF5FileHandler' could not be " "imported.") @@ -48,7 +47,7 @@ def __init__(self, filename, filename_info, filetype_info, **kwargs): self.file_content.update(kwargs) def get_test_content(self, filename, filename_info, filetype_info): - """Mimic reader input file content + """Mimic reader input file content. Args: filename (str): input filename @@ -67,9 +66,10 @@ def get_test_content(self, filename, filename_info, filetype_info): class TestHDF5FileHandler(unittest.TestCase): - """Test HDF5 File Handler Utility class""" + """Test HDF5 File Handler Utility class.""" + def setUp(self): - """Create a test HDF5 file""" + """Create a test HDF5 file.""" import h5py h = h5py.File('test.h5', 'w') # Create Group @@ -107,15 +107,16 @@ def setUp(self): d.attrs['test_attr_str'] = 'test_string' d.attrs['test_attr_int'] = 0 d.attrs['test_attr_float'] = 1.2 + d.attrs['test_ref'] = d.ref h.close() def tearDown(self): - """Remove the previously created test file""" + """Remove the previously created test file.""" os.remove('test.h5') def test_all_basic(self): - """Test everything about the HDF5 class""" + """Test everything about the HDF5 class.""" from satpy.readers.hdf5_utils import HDF5FileHandler import xarray as xr file_handler = HDF5FileHandler('test.h5', {}, {}) @@ -139,9 +140,11 @@ def test_all_basic(self): self.assertTrue('ds2_f' in file_handler) self.assertFalse('fake_ds' in file_handler) + self.assertIsInstance(file_handler['ds2_f/attr/test_ref'], np.ndarray) + def suite(): - """The test suite for test_hdf5_utils.""" + """Test suite for test_hdf5_utils.""" loader = unittest.TestLoader() mysuite = unittest.TestSuite() mysuite.addTest(loader.loadTestsFromTestCase(TestHDF5FileHandler)) From 16dc78196bdc87dc153dee638ec0958f9c35ff2f Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Mon, 4 Nov 2019 19:49:49 +0100 Subject: [PATCH 10/16] Add initial tests for the nwcsaf-msg reader --- satpy/tests/reader_tests/__init__.py | 8 +- satpy/tests/reader_tests/test_hdf5_utils.py | 6 +- satpy/tests/reader_tests/test_nwcsaf_msg.py | 220 ++++++++++++++++++++ 3 files changed, 225 insertions(+), 9 deletions(-) create mode 100644 satpy/tests/reader_tests/test_nwcsaf_msg.py diff --git a/satpy/tests/reader_tests/__init__.py b/satpy/tests/reader_tests/__init__.py index f1d55baffb..1b34c01821 100644 --- a/satpy/tests/reader_tests/__init__.py +++ b/satpy/tests/reader_tests/__init__.py @@ -15,8 +15,7 @@ # # You should have received a copy of the GNU General Public License along with # satpy. If not, see . -"""The reader tests package. -""" +"""The reader tests package.""" import sys @@ -39,7 +38,7 @@ test_electrol_hrit, test_mersi2_l1b, test_avhrr_l1b_gaclac, test_vaisala_gld360, test_fci_l1c_fdhsi, test_tropomi_l2, - test_hsaf_grib) + test_hsaf_grib, test_nwcsaf_msg) if sys.version_info < (2, 7): @@ -49,7 +48,7 @@ def suite(): - """Test suite for all reader tests""" + """Test suite for all reader tests.""" mysuite = unittest.TestSuite() mysuite.addTests(test_abi_l1b.suite()) mysuite.addTests(test_agri_l1.suite()) @@ -95,5 +94,6 @@ def suite(): mysuite.addTests(test_fci_l1c_fdhsi.suite()) mysuite.addTests(test_tropomi_l2.suite()) mysuite.addTests(test_hsaf_grib.suite()) + mysuite.addTests(test_nwcsaf_msg.suite()) return mysuite diff --git a/satpy/tests/reader_tests/test_hdf5_utils.py b/satpy/tests/reader_tests/test_hdf5_utils.py index 115da55dbc..87e76f76cf 100644 --- a/satpy/tests/reader_tests/test_hdf5_utils.py +++ b/satpy/tests/reader_tests/test_hdf5_utils.py @@ -18,7 +18,6 @@ """Module for testing the satpy.readers.hdf5_utils module.""" import os -import sys import numpy as np try: @@ -27,10 +26,7 @@ # fake the import so we can at least run the tests in this file HDF5FileHandler = object -if sys.version_info < (2, 7): - import unittest2 as unittest -else: - import unittest +import unittest class FakeHDF5FileHandler(HDF5FileHandler): diff --git a/satpy/tests/reader_tests/test_nwcsaf_msg.py b/satpy/tests/reader_tests/test_nwcsaf_msg.py new file mode 100644 index 0000000000..c57aa1ef5d --- /dev/null +++ b/satpy/tests/reader_tests/test_nwcsaf_msg.py @@ -0,0 +1,220 @@ +#!/usr/bin/env python +# Copyright (c) 2019 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 . +"""Unittests for NWC SAF MSG (2013) reader.""" + +import unittest +import numpy as np +import tempfile +import os +import h5py + +try: + from unittest import mock +except ImportError: + import mock # noqa + + +fake_ct = { + "01-PALETTE": { + "attrs": { + "CLASS": b"PALETTE", + "PAL_COLORMODEL": b"RGB", + "PAL_TYPE": b"DIRECTINDEX", + }, + "value": np.array( + [ + [100, 100, 100], + [0, 120, 0], + [0, 0, 0], + [250, 190, 250], + [220, 160, 220], + [255, 150, 0], + [255, 100, 0], + [255, 220, 0], + [255, 180, 0], + [255, 255, 140], + [240, 240, 0], + [250, 240, 200], + [215, 215, 150], + [255, 255, 255], + [230, 230, 230], + [0, 80, 215], + [0, 180, 230], + [0, 240, 240], + [90, 200, 160], + [200, 0, 200], + [95, 60, 30], + ], + dtype=np.uint8, + ), + }, + "02-PALETTE": { + "attrs": { + "CLASS": b"PALETTE", + "PAL_COLORMODEL": b"RGB", + "PAL_TYPE": b"DIRECTINDEX", + }, + "value": np.array( + [[100, 100, 100], [255, 100, 0], [0, 80, 215], [95, 60, 30]], dtype=np.uint8 + ), + }, + "CT": { + "attrs": { + "CLASS": b"IMAGE", + "ID": b"CT", + "IMAGE_COLORMODEL": b"RGB", + "IMAGE_SUBCLASS": b"IMAGE_INDEXED", + "IMAGE_VERSION": b"1.0", + "N_COLS": 3712, + "N_LINES": 1856, + "OFFSET": 0.0, + "PALETTE": " 01-PALETTE", + "PRODUCT": b"CT__", + "SCALING_FACTOR": 1.0, + }, + "value": (np.random.rand(1856, 3712) * 255).astype(np.uint8), + }, + "CT_PHASE": { + "attrs": { + "CLASS": b"IMAGE", + "ID": b"CT_PHASE", + "IMAGE_COLORMODEL": b"RGB", + "IMAGE_SUBCLASS": b"IMAGE_INDEXED", + "IMAGE_VERSION": b"1.0", + "N_COLS": 3712, + "N_LINES": 1856, + "OFFSET": 0.0, + "PALETTE": " 02-PALETTE", + "PRODUCT": b"CT__", + "SCALING_FACTOR": 1.0, + }, + "value": (np.random.rand(1856, 3712) * 255).astype(np.uint8), + }, + "CT_QUALITY": { + "attrs": { + "CLASS": b"IMAGE", + "ID": b"CT_QUALITY", + "IMAGE_COLORMODEL": b"RGB", + "IMAGE_SUBCLASS": b"IMAGE_INDEXED", + "IMAGE_VERSION": b"1.0", + "N_COLS": 3712, + "N_LINES": 1856, + "OFFSET": 0.0, + "PRODUCT": b"CT__", + "SCALING_FACTOR": 1.0, + }, + "value": (np.random.rand(1856, 3712) * 65535).astype(np.uint16), + }, + "attrs": { + "CFAC": 13642337, + "COFF": 1856, + "GP_SC_ID": 323, + "IMAGE_ACQUISITION_TIME": b"201611090800", + "LFAC": 13642337, + "LOFF": 1856, + "NB_PARAMETERS": 3, + "NC": 3712, + "NL": 1856, + "NOMINAL_PRODUCT_TIME": b"201611090814", + "PACKAGE": b"SAFNWC/MSG", + "PRODUCT_ALGORITHM_VERSION": b" 2.2", + "PRODUCT_NAME": b"CT__", + "PROJECTION_NAME": b"GEOS<+000.0>", + "REGION_NAME": b"MSG-N", + "SAF": b"NWC", + "SGS_PRODUCT_COMPLETENESS": 99, + "SGS_PRODUCT_QUALITY": 79, + "SPECTRAL_CHANNEL_ID": 0, + }, +} + + +PROJ_KM = { + "gdal_projection": "+proj=geos +a=6378.137000 +b=6356.752300 +lon_0=0.000000 +h=35785.863000", + "gdal_xgeo_up_left": -5569500.0, + "gdal_ygeo_up_left": 5437500.0, + "gdal_xgeo_low_right": 5566500.0, + "gdal_ygeo_low_right": 2653500.0, +} +PROJ = { + "gdal_projection": "+proj=geos +a=6378137.000 +b=6356752.300 +lon_0=0.000000 +h=35785863.000", + "gdal_xgeo_up_left": -5569500.0, + "gdal_ygeo_up_left": 5437500.0, + "gdal_xgeo_low_right": 5566500.0, + "gdal_ygeo_low_right": 2653500.0, +} + + +class TestH5NWCSAF(unittest.TestCase): + """Test the nwcsaf msg reader.""" + + def setUp(self): + """Set up the tests.""" + self.filename = os.path.join( + tempfile.gettempdir(), + "SAFNWC_MSG3_CT___201611090800_MSG-N_______.PLAX.CTTH.0.h5") + + def fill_h5(root, stuff): + for key, val in stuff.items(): + if key in ["value", "attrs"]: + continue + if "value" in val: + root[key] = val["value"] + else: + grp = root.create_group(key) + fill_h5(grp, stuff[key]) + if "attrs" in val: + for attrs, val in val["attrs"].items(): + if isinstance(val, str) and val.startswith(""): + root[key].attrs[attrs] = root[val[24:]].ref + else: + root[key].attrs[attrs] = val + + h5f = h5py.File(self.filename, mode="w") + fill_h5(h5f, fake_ct) + for attr, val in fake_ct["attrs"].items(): + h5f.attrs[attr] = val + h5f.close() + + def test_get_dataset(self): + """Retrieve datasets from a DNB file.""" + from satpy.readers.nwcsaf_msg2013_hdf5 import Hdf5NWCSAF + from satpy import DatasetID + + filename_info = {} + filetype_info = {} + dsid = DatasetID(name='ct') + test = Hdf5NWCSAF(self.filename, filename_info, filetype_info) + ds = test.get_dataset(dsid, {'file_key': 'CT'}) + self.assertEqual(ds.shape, (1856, 3712)) + self.assertEqual(ds.dtype, np.uint8) + + def tearDown(self): + """Destroy.""" + try: + os.remove(self.filename) + except OSError: + pass + + +def suite(): + """Test suite for test_writers.""" + loader = unittest.TestLoader() + my_suite = unittest.TestSuite() + my_suite.addTest(loader.loadTestsFromTestCase(TestH5NWCSAF)) + + return my_suite From 4a2e44d7d46c5aa43e83deeea6dbaf11e21d2339 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Mon, 4 Nov 2019 20:46:39 +0100 Subject: [PATCH 11/16] Make the dummy data in the nwcsaf-msg tests ordered --- satpy/tests/reader_tests/test_nwcsaf_msg.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/satpy/tests/reader_tests/test_nwcsaf_msg.py b/satpy/tests/reader_tests/test_nwcsaf_msg.py index c57aa1ef5d..19a2770108 100644 --- a/satpy/tests/reader_tests/test_nwcsaf_msg.py +++ b/satpy/tests/reader_tests/test_nwcsaf_msg.py @@ -21,6 +21,7 @@ import tempfile import os import h5py +from collections import OrderedDict try: from unittest import mock @@ -142,6 +143,7 @@ }, } +fake_ct = OrderedDict(sorted(fake_ct.items(), key=lambda t: t[0])) PROJ_KM = { "gdal_projection": "+proj=geos +a=6378.137000 +b=6356.752300 +lon_0=0.000000 +h=35785.863000", From 59faf6df2e4c1d66f4749ac2a617dbbacc19c78a Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Mon, 4 Nov 2019 21:07:01 +0100 Subject: [PATCH 12/16] Add a CTTH test to the nwcsaf-msg reader tests --- satpy/tests/reader_tests/test_nwcsaf_msg.py | 280 +++++++++++++++++++- 1 file changed, 272 insertions(+), 8 deletions(-) diff --git a/satpy/tests/reader_tests/test_nwcsaf_msg.py b/satpy/tests/reader_tests/test_nwcsaf_msg.py index 19a2770108..8c1ffceec8 100644 --- a/satpy/tests/reader_tests/test_nwcsaf_msg.py +++ b/satpy/tests/reader_tests/test_nwcsaf_msg.py @@ -145,6 +145,247 @@ fake_ct = OrderedDict(sorted(fake_ct.items(), key=lambda t: t[0])) +fake_ctth = { + "01-PALETTE": { + "attrs": { + "CLASS": b"PALETTE", + "PAL_COLORMODEL": b"RGB", + "PAL_TYPE": b"DIRECTINDEX", + }, + "value": np.array( + [ + [0, 0, 0], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [255, 255, 255], + [239, 239, 223], + [239, 239, 223], + [238, 214, 210], + [238, 214, 210], + [0, 255, 255], + [0, 255, 255], + [0, 216, 255], + [0, 216, 255], + [0, 178, 255], + [0, 178, 255], + [0, 140, 48], + [0, 140, 48], + [0, 255, 0], + [0, 255, 0], + [153, 255, 0], + [153, 255, 0], + [178, 255, 0], + [178, 255, 0], + [216, 255, 0], + [216, 255, 0], + [255, 255, 0], + [255, 255, 0], + [255, 216, 0], + [255, 216, 0], + [255, 164, 0], + [255, 164, 0], + [255, 102, 0], + [255, 102, 0], + [255, 76, 0], + [255, 76, 0], + [178, 51, 0], + [178, 51, 0], + [153, 20, 47], + [153, 20, 47], + [126, 0, 43], + [126, 0, 43], + [255, 0, 216], + [255, 0, 216], + [255, 0, 128], + [255, 0, 128], + [255, 0, 128], + [255, 0, 128], + [255, 0, 128], + [255, 0, 128], + [255, 0, 128], + [255, 0, 128], + [255, 0, 128], + [255, 0, 128], + [255, 0, 128], + ], + dtype=np.uint8, + ), + }, + "02-PALETTE": { + "attrs": { + "CLASS": b"PALETTE", + "PAL_COLORMODEL": b"RGB", + "PAL_TYPE": b"DIRECTINDEX", + }, + "value": (np.random.rand(128, 3) * 255).astype(np.uint8), + }, + "03-PALETTE": { + "attrs": { + "CLASS": b"PALETTE", + "PAL_COLORMODEL": b"RGB", + "PAL_TYPE": b"DIRECTINDEX", + }, + "value": (np.random.rand(256, 3) * 255).astype(np.uint8), + }, + "04-PALETTE": { + "attrs": { + "CLASS": b"PALETTE", + "PAL_COLORMODEL": b"RGB", + "PAL_TYPE": b"DIRECTINDEX", + }, + "value": np.array( + [ + [78, 119, 145], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [12, 12, 12], + [24, 24, 24], + [36, 36, 36], + [48, 48, 48], + [60, 60, 60], + [72, 72, 72], + [84, 84, 84], + [96, 96, 96], + [108, 108, 108], + [120, 120, 120], + [132, 132, 132], + [144, 144, 144], + [156, 156, 156], + [168, 168, 168], + [180, 180, 180], + [192, 192, 192], + [204, 204, 204], + [216, 216, 216], + [228, 228, 228], + [240, 240, 240], + [240, 240, 240], + ], + dtype=np.uint8, + ), + }, + "CTTH_EFFECT": { + "attrs": { + "CLASS": b"IMAGE", + "ID": b"CTTH_EFFECT", + "IMAGE_COLORMODEL": b"RGB", + "IMAGE_SUBCLASS": b"IMAGE_INDEXED", + "IMAGE_VERSION": b"1.0", + "N_COLS": 3712, + "N_LINES": 1856, + "OFFSET": -50.0, + "PALETTE": " 04-PALETTE", + "PRODUCT": b"CTTH", + "SCALING_FACTOR": 5.0, + }, + "value": (np.random.rand(1856, 3712) * 255).astype(np.uint8), + }, + "CTTH_HEIGHT": { + "attrs": { + "CLASS": b"IMAGE", + "ID": b"CTTH_HEIGHT", + "IMAGE_COLORMODEL": b"RGB", + "IMAGE_SUBCLASS": b"IMAGE_INDEXED", + "IMAGE_VERSION": b"1.0", + "N_COLS": 3712, + "N_LINES": 1856, + "OFFSET": -2000.0, + "PALETTE": " 02-PALETTE", + "PRODUCT": b"CTTH", + "SCALING_FACTOR": 200.0, + }, + "value": (np.random.rand(1856, 3712) * 255).astype(np.uint8), + }, + "CTTH_PRESS": { + "attrs": { + "CLASS": b"IMAGE", + "ID": b"CTTH_PRESS", + "IMAGE_COLORMODEL": b"RGB", + "IMAGE_SUBCLASS": b"IMAGE_INDEXED", + "IMAGE_VERSION": b"1.0", + "N_COLS": 3712, + "N_LINES": 1856, + "OFFSET": -250.0, + "PALETTE": " 01-PALETTE", + "PRODUCT": b"CTTH", + "SCALING_FACTOR": 25.0, + }, + "value": (np.random.rand(1856, 3712) * 255).astype(np.uint8), + }, + "CTTH_QUALITY": { + "attrs": { + "CLASS": b"IMAGE", + "ID": b"CTTH_QUALITY", + "IMAGE_COLORMODEL": b"RGB", + "IMAGE_SUBCLASS": b"IMAGE_INDEXED", + "IMAGE_VERSION": b"1.0", + "N_COLS": 3712, + "N_LINES": 1856, + "OFFSET": 0.0, + "PRODUCT": b"CTTH", + "SCALING_FACTOR": 1.0, + }, + "value": (np.random.rand(1856, 3712) * 65535).astype(np.uint16), + }, + "CTTH_TEMPER": { + "attrs": { + "CLASS": b"IMAGE", + "ID": b"CTTH_TEMPER", + "IMAGE_COLORMODEL": b"RGB", + "IMAGE_SUBCLASS": b"IMAGE_INDEXED", + "IMAGE_VERSION": b"1.0", + "N_COLS": 3712, + "N_LINES": 1856, + "OFFSET": 150.0, + "PALETTE": " 03-PALETTE", + "PRODUCT": b"CTTH", + "SCALING_FACTOR": 1.0, + }, + "value": (np.random.rand(1856, 3712) * 255).astype(np.uint8), + }, + "attrs": { + "CFAC": 13642337, + "COFF": 1856, + "GP_SC_ID": 323, + "IMAGE_ACQUISITION_TIME": b"201611090800", + "LFAC": 13642337, + "LOFF": 1856, + "NB_PARAMETERS": 5, + "NC": 3712, + "NL": 1856, + "NOMINAL_PRODUCT_TIME": b"201611090816", + "PACKAGE": b"SAFNWC/MSG", + "PRODUCT_ALGORITHM_VERSION": b" 2.2", + "PRODUCT_NAME": b"CTTH", + "PROJECTION_NAME": b"GEOS<+000.0>", + "REGION_NAME": b"MSG-N", + "SAF": b"NWC", + "SGS_PRODUCT_COMPLETENESS": 87, + "SGS_PRODUCT_QUALITY": 69, + "SPECTRAL_CHANNEL_ID": 0, + }, +} + +fake_ctth = OrderedDict(sorted(fake_ctth.items(), key=lambda t: t[0])) + PROJ_KM = { "gdal_projection": "+proj=geos +a=6378.137000 +b=6356.752300 +lon_0=0.000000 +h=35785.863000", "gdal_xgeo_up_left": -5569500.0, @@ -166,9 +407,15 @@ class TestH5NWCSAF(unittest.TestCase): def setUp(self): """Set up the tests.""" - self.filename = os.path.join( + self.filename_ct = os.path.join( tempfile.gettempdir(), - "SAFNWC_MSG3_CT___201611090800_MSG-N_______.PLAX.CTTH.0.h5") + "SAFNWC_MSG3_CT___201611090800_MSG-N_______.PLAX.CTTH.0.h5", + ) + + self.filename_ctth = os.path.join( + tempfile.gettempdir(), + "SAFNWC_MSG3_CTTH_201611090800_MSG-N_______.PLAX.CTTH.0.h5", + ) def fill_h5(root, stuff): for key, val in stuff.items(): @@ -181,17 +428,25 @@ def fill_h5(root, stuff): fill_h5(grp, stuff[key]) if "attrs" in val: for attrs, val in val["attrs"].items(): - if isinstance(val, str) and val.startswith(""): + if isinstance(val, str) and val.startswith( + "" + ): root[key].attrs[attrs] = root[val[24:]].ref else: root[key].attrs[attrs] = val - h5f = h5py.File(self.filename, mode="w") + h5f = h5py.File(self.filename_ct, mode="w") fill_h5(h5f, fake_ct) for attr, val in fake_ct["attrs"].items(): h5f.attrs[attr] = val h5f.close() + h5f = h5py.File(self.filename_ctth, mode="w") + fill_h5(h5f, fake_ctth) + for attr, val in fake_ctth["attrs"].items(): + h5f.attrs[attr] = val + h5f.close() + def test_get_dataset(self): """Retrieve datasets from a DNB file.""" from satpy.readers.nwcsaf_msg2013_hdf5 import Hdf5NWCSAF @@ -199,16 +454,25 @@ def test_get_dataset(self): filename_info = {} filetype_info = {} - dsid = DatasetID(name='ct') - test = Hdf5NWCSAF(self.filename, filename_info, filetype_info) - ds = test.get_dataset(dsid, {'file_key': 'CT'}) + dsid = DatasetID(name="ct") + test = Hdf5NWCSAF(self.filename_ct, filename_info, filetype_info) + ds = test.get_dataset(dsid, {"file_key": "CT"}) self.assertEqual(ds.shape, (1856, 3712)) self.assertEqual(ds.dtype, np.uint8) + filename_info = {} + filetype_info = {} + dsid = DatasetID(name="ctth_alti") + test = Hdf5NWCSAF(self.filename_ctth, filename_info, filetype_info) + ds = test.get_dataset(dsid, {"file_key": "CTTH_HEIGHT"}) + self.assertEqual(ds.shape, (1856, 3712)) + self.assertEqual(ds.dtype, np.float32) + def tearDown(self): """Destroy.""" try: - os.remove(self.filename) + os.remove(self.filename_ct) + os.remove(self.filename_ctth) except OSError: pass From 1eadde8aa58e77505c31b35034db6e3fc4479e22 Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Mon, 4 Nov 2019 23:55:18 +0100 Subject: [PATCH 13/16] Add test of get_area_def and add more content to the dataset testing Signed-off-by: Adam Dybbroe --- satpy/tests/reader_tests/test_nwcsaf_msg.py | 49 +++++++++++++++++++-- 1 file changed, 46 insertions(+), 3 deletions(-) diff --git a/satpy/tests/reader_tests/test_nwcsaf_msg.py b/satpy/tests/reader_tests/test_nwcsaf_msg.py index 8c1ffceec8..ebf6e07f41 100644 --- a/satpy/tests/reader_tests/test_nwcsaf_msg.py +++ b/satpy/tests/reader_tests/test_nwcsaf_msg.py @@ -14,6 +14,7 @@ # # You should have received a copy of the GNU General Public License along with # satpy. If not, see . +# """Unittests for NWC SAF MSG (2013) reader.""" import unittest @@ -28,6 +29,17 @@ except ImportError: import mock # noqa +CTYPE_TEST_ARRAY = (np.random.rand(1856, 3712) * 255).astype(np.uint8) +CTYPE_TEST_FRAME = (np.arange(100).reshape(10, 10) / 100. * 20).astype(np.uint8) +CTYPE_TEST_ARRAY[1000:1010, 1000:1010] = CTYPE_TEST_FRAME + +CTTH_HEIGHT_TEST_ARRAY = (np.random.rand(1856, 3712) * 255).astype(np.uint8) +_CTTH_HEIGHT_TEST_FRAME = (np.arange(100).reshape(10, 10) / 100. * 80).astype(np.uint8) +CTTH_HEIGHT_TEST_ARRAY[1000:1010, 1000:1010] = _CTTH_HEIGHT_TEST_FRAME + +CTTH_HEIGHT_TEST_FRAME_RES = _CTTH_HEIGHT_TEST_FRAME.astype(np.float32) * 200 - 2000 +CTTH_HEIGHT_TEST_FRAME_RES[0, 0:10] = np.nan +CTTH_HEIGHT_TEST_FRAME_RES[1, 0:3] = np.nan fake_ct = { "01-PALETTE": { @@ -87,7 +99,7 @@ "PRODUCT": b"CT__", "SCALING_FACTOR": 1.0, }, - "value": (np.random.rand(1856, 3712) * 255).astype(np.uint8), + "value": (CTYPE_TEST_ARRAY), }, "CT_PHASE": { "attrs": { @@ -312,7 +324,7 @@ "PRODUCT": b"CTTH", "SCALING_FACTOR": 200.0, }, - "value": (np.random.rand(1856, 3712) * 255).astype(np.uint8), + "value": (CTTH_HEIGHT_TEST_ARRAY), }, "CTTH_PRESS": { "attrs": { @@ -401,6 +413,15 @@ "gdal_ygeo_low_right": 2653500.0, } +AREA_DEF_DICT = { + "proj_dict": {'proj': 'geos', 'lon_0': 0, 'h': 35785831, 'x_0': 0, 'y_0': 0, + 'a': 6378169, 'b': 6356583.8, 'units': 'm', 'no_defs': None, 'type': 'crs'}, + "area_id": 'MSG-N', + "x_size": 3712, + "y_size": 1856, + "area_extent": (-5570248.2825, 1501.0099, 5567247.8793, 5570247.8784) +} + class TestH5NWCSAF(unittest.TestCase): """Test the nwcsaf msg reader.""" @@ -447,8 +468,28 @@ def fill_h5(root, stuff): h5f.attrs[attr] = val h5f.close() + def test_get_area_def(self): + """Get the area definition.""" + from satpy.readers.nwcsaf_msg2013_hdf5 import Hdf5NWCSAF + from satpy import DatasetID + + filename_info = {} + filetype_info = {} + dsid = DatasetID(name="ct") + test = Hdf5NWCSAF(self.filename_ct, filename_info, filetype_info) + + area_def = test.get_area_def(dsid) + + aext_res = AREA_DEF_DICT['area_extent'] + for i in range(4): + self.assertAlmostEqual(area_def.area_extent[i], aext_res[i], 4) + + proj_dict = AREA_DEF_DICT['proj_dict'] + for key in proj_dict: + self.assertEqual(proj_dict[key], area_def.proj_dict[key]) + def test_get_dataset(self): - """Retrieve datasets from a DNB file.""" + """Retrieve datasets from a NWCSAF msgv2013 hdf5 file.""" from satpy.readers.nwcsaf_msg2013_hdf5 import Hdf5NWCSAF from satpy import DatasetID @@ -459,6 +500,7 @@ def test_get_dataset(self): ds = test.get_dataset(dsid, {"file_key": "CT"}) self.assertEqual(ds.shape, (1856, 3712)) self.assertEqual(ds.dtype, np.uint8) + np.testing.assert_allclose(ds.data.compute()[1000:1010, 1000:1010], CTYPE_TEST_FRAME) filename_info = {} filetype_info = {} @@ -467,6 +509,7 @@ def test_get_dataset(self): ds = test.get_dataset(dsid, {"file_key": "CTTH_HEIGHT"}) self.assertEqual(ds.shape, (1856, 3712)) self.assertEqual(ds.dtype, np.float32) + np.testing.assert_allclose(ds.data.compute()[1000:1010, 1000:1010], CTTH_HEIGHT_TEST_FRAME_RES) def tearDown(self): """Destroy.""" From a4b6a17aca51a1590c3171f10162316e90de16b5 Mon Sep 17 00:00:00 2001 From: Martin Raspaud Date: Tue, 5 Nov 2019 08:48:03 +0100 Subject: [PATCH 14/16] Use only DataArrays in the nwcsaf-msg dataset reader --- satpy/readers/nwcsaf_msg2013_hdf5.py | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/satpy/readers/nwcsaf_msg2013_hdf5.py b/satpy/readers/nwcsaf_msg2013_hdf5.py index 7f29527556..69c987dc66 100644 --- a/satpy/readers/nwcsaf_msg2013_hdf5.py +++ b/satpy/readers/nwcsaf_msg2013_hdf5.py @@ -20,10 +20,9 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -"""Reader for the old NWCSAF/Geo (v2013 and earlier) cloud product format +"""Reader for the old NWCSAF/Geo (v2013 and earlier) cloud product format. References: - - The NWCSAF GEO 2013 products documentation: http://www.nwcsaf.org/web/guest/archive - Search for Code "ICD/3"; Type "MSG" and the box to the right should say 'Status' (which means any @@ -31,7 +30,6 @@ http://www.nwcsaf.org/aemetRest/downloadAttachment/2623 - """ import logging @@ -79,13 +77,13 @@ def get_dataset(self, dataset_id, ds_info): data.attrs['valid_range'] = (0, 20) data.attrs['_FillValue'] = 255 # data.attrs['palette_meanings'] = list(range(21)) - - scaled_data = (data.data * data.attrs['SCALING_FACTOR'] + data.attrs['OFFSET']).astype(dtype) + attrs = data.attrs + data = (data * data.attrs['SCALING_FACTOR'] + data.attrs['OFFSET']).astype(dtype) if nodata: - scaled_data = np.where(data.data == nodata, np.nan, scaled_data) - scaled_data = np.where(scaled_data < 0, np.nan, scaled_data) + data = data.where(data != nodata) + data = data.where(data >= 0) - data.data = scaled_data + data.attrs = attrs for key in list(data.attrs.keys()): val = data.attrs[key] @@ -95,10 +93,7 @@ def get_dataset(self, dataset_id, ds_info): return data def get_area_def(self, dsid): - """Get the area definition of the datasets in the file. - - """ - + """Get the area definition of the datasets in the file.""" if dsid.name.endswith('_pal'): raise NotImplementedError @@ -138,9 +133,7 @@ def start_time(self): def get_area_extent(cfac, lfac, coff, loff, numcols, numlines): - """Get the area extent from msg parameters. - """ - + """Get the area extent from msg parameters.""" xur = (numcols - coff) * 2 ** 16 / (cfac * 1.0) xur = np.deg2rad(xur) * 35785831.0 xll = (-1 - coff) * 2 ** 16 / (cfac * 1.0) From 8416a8f0624014cf5709b41be89cfe35b3b186a0 Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Tue, 5 Nov 2019 08:51:31 +0100 Subject: [PATCH 15/16] Enhance test of get_area_def Signed-off-by: Adam Dybbroe --- satpy/tests/reader_tests/test_nwcsaf_msg.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/satpy/tests/reader_tests/test_nwcsaf_msg.py b/satpy/tests/reader_tests/test_nwcsaf_msg.py index ebf6e07f41..a7004d88b3 100644 --- a/satpy/tests/reader_tests/test_nwcsaf_msg.py +++ b/satpy/tests/reader_tests/test_nwcsaf_msg.py @@ -488,6 +488,11 @@ def test_get_area_def(self): for key in proj_dict: self.assertEqual(proj_dict[key], area_def.proj_dict[key]) + self.assertEqual(AREA_DEF_DICT['x_size'], area_def.width) + self.assertEqual(AREA_DEF_DICT['y_size'], area_def.height) + + self.assertEqual(AREA_DEF_DICT['area_id'], area_def.area_id) + def test_get_dataset(self): """Retrieve datasets from a NWCSAF msgv2013 hdf5 file.""" from satpy.readers.nwcsaf_msg2013_hdf5 import Hdf5NWCSAF From 462af23aae006ed75d08f429c45a9fbbba439cf9 Mon Sep 17 00:00:00 2001 From: Adam Dybbroe Date: Tue, 5 Nov 2019 10:53:02 +0100 Subject: [PATCH 16/16] Fix reader to handle nodata correctly when scaling, and grow test coverage Signed-off-by: Adam Dybbroe --- satpy/readers/nwcsaf_msg2013_hdf5.py | 9 ++-- satpy/tests/reader_tests/test_nwcsaf_msg.py | 49 ++++++++++++++++++--- 2 files changed, 48 insertions(+), 10 deletions(-) diff --git a/satpy/readers/nwcsaf_msg2013_hdf5.py b/satpy/readers/nwcsaf_msg2013_hdf5.py index 69c987dc66..f65253c641 100644 --- a/satpy/readers/nwcsaf_msg2013_hdf5.py +++ b/satpy/readers/nwcsaf_msg2013_hdf5.py @@ -77,12 +77,13 @@ def get_dataset(self, dataset_id, ds_info): data.attrs['valid_range'] = (0, 20) data.attrs['_FillValue'] = 255 # data.attrs['palette_meanings'] = list(range(21)) + attrs = data.attrs - data = (data * data.attrs['SCALING_FACTOR'] + data.attrs['OFFSET']).astype(dtype) + scaled_data = (data * data.attrs['SCALING_FACTOR'] + data.attrs['OFFSET']).astype(dtype) if nodata: - data = data.where(data != nodata) - data = data.where(data >= 0) - + scaled_data = scaled_data.where(data != nodata) + scaled_data = scaled_data.where(scaled_data >= 0) + data = scaled_data data.attrs = attrs for key in list(data.attrs.keys()): diff --git a/satpy/tests/reader_tests/test_nwcsaf_msg.py b/satpy/tests/reader_tests/test_nwcsaf_msg.py index a7004d88b3..f7fdec168f 100644 --- a/satpy/tests/reader_tests/test_nwcsaf_msg.py +++ b/satpy/tests/reader_tests/test_nwcsaf_msg.py @@ -41,6 +41,23 @@ CTTH_HEIGHT_TEST_FRAME_RES[0, 0:10] = np.nan CTTH_HEIGHT_TEST_FRAME_RES[1, 0:3] = np.nan +CTTH_PRESSURE_TEST_ARRAY = (np.random.rand(1856, 3712) * 255).astype(np.uint8) +_CTTH_PRESSURE_TEST_FRAME = (np.arange(100).reshape(10, 10) / 100. * 54).astype(np.uint8) +CTTH_PRESSURE_TEST_ARRAY[1000:1010, 1000:1010] = _CTTH_PRESSURE_TEST_FRAME + +CTTH_PRESSURE_TEST_FRAME_RES = _CTTH_PRESSURE_TEST_FRAME.astype(np.float32) * 25 - 250 +CTTH_PRESSURE_TEST_FRAME_RES[0, 0:10] = np.nan +CTTH_PRESSURE_TEST_FRAME_RES[1, 0:9] = np.nan + +CTTH_TEMPERATURE_TEST_ARRAY = (np.random.rand(1856, 3712) * 255).astype(np.uint8) +_CTTH_TEMPERATURE_TEST_FRAME = (np.arange(100).reshape(10, 10) / 100. * 140).astype(np.uint8) +_CTTH_TEMPERATURE_TEST_FRAME[8, 5] = 255 +CTTH_TEMPERATURE_TEST_ARRAY[1000:1010, 1000:1010] = _CTTH_TEMPERATURE_TEST_FRAME + +CTTH_TEMPERATURE_TEST_FRAME_RES = _CTTH_TEMPERATURE_TEST_FRAME.astype(np.float32) * 1.0 + 150 +CTTH_TEMPERATURE_TEST_FRAME_RES[8, 5] = np.nan + + fake_ct = { "01-PALETTE": { "attrs": { @@ -340,7 +357,7 @@ "PRODUCT": b"CTTH", "SCALING_FACTOR": 25.0, }, - "value": (np.random.rand(1856, 3712) * 255).astype(np.uint8), + "value": (CTTH_PRESSURE_TEST_ARRAY), }, "CTTH_QUALITY": { "attrs": { @@ -371,7 +388,7 @@ "PRODUCT": b"CTTH", "SCALING_FACTOR": 1.0, }, - "value": (np.random.rand(1856, 3712) * 255).astype(np.uint8), + "value": (CTTH_TEMPERATURE_TEST_ARRAY), }, "attrs": { "CFAC": 13642337, @@ -485,8 +502,10 @@ def test_get_area_def(self): self.assertAlmostEqual(area_def.area_extent[i], aext_res[i], 4) proj_dict = AREA_DEF_DICT['proj_dict'] - for key in proj_dict: - self.assertEqual(proj_dict[key], area_def.proj_dict[key]) + self.assertEqual(proj_dict['proj'], area_def.proj_dict['proj']) + # Not all elements passed on Appveyor, so skip testing every single element of the proj-dict: + # for key in proj_dict: + # self.assertEqual(proj_dict[key], area_def.proj_dict[key]) self.assertEqual(AREA_DEF_DICT['x_size'], area_def.width) self.assertEqual(AREA_DEF_DICT['y_size'], area_def.height) @@ -505,7 +524,7 @@ def test_get_dataset(self): ds = test.get_dataset(dsid, {"file_key": "CT"}) self.assertEqual(ds.shape, (1856, 3712)) self.assertEqual(ds.dtype, np.uint8) - np.testing.assert_allclose(ds.data.compute()[1000:1010, 1000:1010], CTYPE_TEST_FRAME) + np.testing.assert_allclose(ds.data[1000:1010, 1000:1010].compute(), CTYPE_TEST_FRAME) filename_info = {} filetype_info = {} @@ -514,7 +533,25 @@ def test_get_dataset(self): ds = test.get_dataset(dsid, {"file_key": "CTTH_HEIGHT"}) self.assertEqual(ds.shape, (1856, 3712)) self.assertEqual(ds.dtype, np.float32) - np.testing.assert_allclose(ds.data.compute()[1000:1010, 1000:1010], CTTH_HEIGHT_TEST_FRAME_RES) + np.testing.assert_allclose(ds.data[1000:1010, 1000:1010].compute(), CTTH_HEIGHT_TEST_FRAME_RES) + + filename_info = {} + filetype_info = {} + dsid = DatasetID(name="ctth_pres") + test = Hdf5NWCSAF(self.filename_ctth, filename_info, filetype_info) + ds = test.get_dataset(dsid, {"file_key": "CTTH_PRESS"}) + self.assertEqual(ds.shape, (1856, 3712)) + self.assertEqual(ds.dtype, np.float32) + np.testing.assert_allclose(ds.data[1000:1010, 1000:1010].compute(), CTTH_PRESSURE_TEST_FRAME_RES) + + filename_info = {} + filetype_info = {} + dsid = DatasetID(name="ctth_tempe") + test = Hdf5NWCSAF(self.filename_ctth, filename_info, filetype_info) + ds = test.get_dataset(dsid, {"file_key": "CTTH_TEMPER"}) + self.assertEqual(ds.shape, (1856, 3712)) + self.assertEqual(ds.dtype, np.float32) + np.testing.assert_allclose(ds.data[1000:1010, 1000:1010].compute(), CTTH_TEMPERATURE_TEST_FRAME_RES) def tearDown(self): """Destroy."""