Merge pull request #886 from adybbroe/feature-old-nwcsaf-geo-format

Reader for NWCSAF/MSG 2013 format

satpy/composites/__init__.py

@@ -931,6 +931,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
 

satpy/etc/readers/nwcsaf-msg2013-hdf5.yaml

@@ -0,0 +1,190 @@
+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
+    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_effective_cloudiness:
+    name: ctth_effective_cloudiness
+    sensor: seviri
+    resolution: 3000
+    file_type: h5_nwcsaf_ctth
+    file_key: CTTH_EFFECT
+
+  ctth_effective_cloudiness_pal:
+    name: ctth_eff_pal
+    sensor: seviri
+    resolution: 3000
+    file_type: h5_nwcsaf_ctth
+    file_key: 04-PALETTE
+
+  ctth_quality:
+    name: ctth_quality
+    sensor: seviri
+    resolution: 3000
+    file_type: h5_nwcsaf_ctth
+    file_key: CTTH_QUALITY
+

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

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:
@@ -59,6 +60,20 @@ def _collect_attrs(self, name, attrs):
                 self.file_content[fc_key] = np2str(value)
             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 isinstance(hf[name].attrs[key], h5py.h5r.Reference):
+                ref_name = h5py.h5r.get_name(hf[name].attrs[key], hf.id)
+                return hf[ref_name][()]
 
     def collect_metadata(self, name, obj):
         if isinstance(obj, h5py.Dataset):

satpy/readers/nwcsaf_msg2013_hdf5.py

@@ -0,0 +1,151 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# Copyright (c) 2019 Pytroll
+
+# Author(s):
+
+#   Adam.Dybbroe <[email protected]>
+
+# 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 <http://www.gnu.org/licenses/>.
+
+"""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
+from datetime import datetime
+import numpy as np
+from satpy.readers.hdf5_utils import HDF5FileHandler
+from pyresample.geometry import AreaDefinition
+import h5py
+
+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]
+
+        nodata = None
+        if 'SCALING_FACTOR' in data.attrs and 'OFFSET' in data.attrs:
+            dtype = np.dtype(data.data)
+            if dataset_id.name in ['ctth_alti']:
+                data.attrs['valid_range'] = (0, 27000)
+                data.attrs['_FillValue'] = np.nan
+
+            if dataset_id.name in ['ctth_alti', 'ctth_pres', 'ctth_tempe', 'ctth_effective_cloudiness']:
+                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))
+
+            attrs = data.attrs
+            scaled_data = (data * data.attrs['SCALING_FACTOR'] + data.attrs['OFFSET']).astype(dtype)
+            if nodata:
+                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()):
+            val = data.attrs[key]
+            if isinstance(val, h5py.h5r.Reference):
+                del data.attrs[key]
+
+        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')
+
+
+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