Source code for satpy.readers.seviri_l1b_icare

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# 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 <>.
r"""Interface to SEVIRI L1B data from ICARE (Lille).


The ``seviri_l1b_icare`` reader reads MSG-SEVIRI L1.5 image data in HDF format
that has been produced by the ICARE Data and Services Center
Data can be accessed via:

Each SEVIRI timeslot comes as 12 HDF files, one per band. Only those bands
that are of interest need to be passed to the reader. Others can be ignored.
Filenames follow the format: GEO_L1B-MSG1_YYYY-MM-DDTHH-MM-SS_G_CHANN_VX-XX.hdf
YYYY, MM, DD, HH, MM, SS specify the timeslot starting time.
CHANN is the channel (i.e: HRV, IR016, WV073, etc)
VX-XX is the processing version number

Here is an example how to read the data in satpy:

.. code-block:: python

    from satpy import Scene
    import glob

    filenames = glob.glob('data/*2019-03-01T12-00-00*.hdf')
    scn = Scene(filenames=filenames, reader='seviri_l1b_icare')
    scn.load(['VIS006', 'IR_108'])


.. code-block:: none

    <xarray.DataArray 'array-a1d52b7e19ec5a875e2f038df5b60d7e' (y: 3712, x: 3712)>
    dask.array<add, shape=(3712, 3712), dtype=float32, chunksize=(1024, 1024), chunktype=numpy.ndarray>
        crs      object +proj=geos +a=6378169.0 +b=6356583.8 +lon_0=0.0 +h=35785831.0 +units=m +type=crs
      * y        (y) float64 5.566e+06 5.563e+06 5.56e+06 ... -5.566e+06 -5.569e+06
      * x        (x) float64 -5.566e+06 -5.563e+06 -5.56e+06 ... 5.566e+06 5.569e+06
        start_time:           2004-12-29 12:15:00
        end_time:             2004-12-29 12:27:44
        area:                 Area ID: geosmsg\nDescription: MSG/SEVIRI low resol...
        name:                 IR_108
        resolution:           3000.403165817
        calibration:          brightness_temperature
        polarization:         None
        level:                None
        modifiers:            ()
        ancillary_variables:  []

from datetime import datetime

import numpy as np

from satpy.readers._geos_area import get_area_definition, get_area_extent
from satpy.readers.hdf4_utils import HDF4FileHandler

[docs] class SEVIRI_ICARE(HDF4FileHandler): """SEVIRI L1B handler for HDF4 files.""" def __init__(self, filename, filename_info, filetype_info): """Init the file handler.""" super(SEVIRI_ICARE, self).__init__(filename, filename_info, filetype_info) # These are VIS bands self.ref_bands = ['HRV', 'VIS006', 'VIS008', 'IR_016'] # And these are IR bands self.bt_bands = ['IR_039', 'IR_062', 'IR_073', 'IR_087', 'IR_097', 'IR_108', 'IR_120', 'IR_134', 'WV_062', 'WV_073'] @property def sensor_name(self): """Get the sensor name.""" # the sensor and platform names are stored together, eg: MSG1/SEVIRI attr = self['/attr/Sensors'] if isinstance(attr, np.ndarray): attr = str(attr.astype(str)).lower() else: attr = attr.lower() plat = attr[0:4] sens = attr[5:] # icare uses non-standard platform names if plat == 'msg1': plat = 'Meteosat-08' elif plat == 'msg2': plat = 'Meteosat-09' elif plat == 'msg3': plat = 'Meteosat-10' elif plat == 'msg4': plat = 'Meteosat-11' else: raise NameError("Unsupported satellite platform:"+plat) return [plat, sens] @property def satlon(self): """Get the satellite longitude.""" attr = self['/attr/Sub_Satellite_Longitude'] if isinstance(attr, np.ndarray): attr = float(attr.astype(str)) return attr @property def projlon(self): """Get the projection longitude.""" attr = self['/attr/Projection_Longitude'] if isinstance(attr, np.ndarray): attr = float(attr.astype(str)) return attr @property def projection(self): """Get the projection.""" attr = self['/attr/Geographic_Projection'] if isinstance(attr, np.ndarray): attr = str(attr.astype(str)) attr = attr.lower() if attr != 'geos': raise NotImplementedError("Only the GEOS projection is supported.\ This is:", attr) return attr @property def zone(self): """Get the zone.""" attr = self['/attr/Zone'] if isinstance(attr, np.ndarray): attr = str(attr.astype(str)).lower() return attr @property def res(self): """Get the resolution.""" attr = self['/attr/Nadir_Pixel_Size'] if isinstance(attr, np.ndarray): attr = str(attr.astype(str)).lower() return float(attr) @property def end_time(self): """Get the end time.""" attr = self['/attr/End_Acquisition_Date'] if isinstance(attr, np.ndarray): attr = str(attr.astype(str)) # In some versions milliseconds are present, sometimes not. try: endacq = datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ") except ValueError: endacq = datetime.strptime(attr, "%Y-%m-%dT%H:%M:%S.%fZ") return endacq @property def start_time(self): """Get the start time.""" attr = self['/attr/Beginning_Acquisition_Date'] if isinstance(attr, np.ndarray): attr = str(attr.astype(str)) # In some versions milliseconds are present, sometimes not. try: stacq = datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ") except ValueError: stacq = datetime.strptime(attr, "%Y-%m-%dT%H:%M:%S.%fZ") return stacq @property def alt(self): """Get the altitude.""" attr = self['/attr/Altitude'] if isinstance(attr, np.ndarray): attr = attr.astype(str) attr = float(attr) # This is stored in km, convert to m attr = attr * 1000. return attr @property def geoloc(self): """Get the geolocation.""" attr = self['/attr/Geolocation'] if isinstance(attr, np.ndarray): attr = attr.astype(str) cfac = float(attr[0]) coff = float(attr[1]) lfac = float(attr[2]) loff = float(attr[3]) return [cfac, lfac, coff, loff]
[docs] def get_metadata(self, data, ds_info): """Get the metadata.""" mda = {} mda.update(data.attrs) mda.update(ds_info) geoloc = self.geoloc mda.update({ 'start_time': self.start_time, 'end_time': self.end_time, 'platform_name': self.sensor_name[0], 'sensor': self.sensor_name[1], 'zone':, 'projection_altitude': self.alt, 'cfac': geoloc[0], 'lfac': geoloc[1], 'coff': geoloc[2], 'loff': geoloc[3], 'resolution': self.res, 'satellite_actual_longitude': self.satlon, 'projection_longitude': self.projlon, 'projection_type': self.projection }) return mda
[docs] def _get_dsname(self, ds_id): """Return the correct dataset name based on requested band.""" if ds_id['name'] in self.ref_bands: ds_get_name = 'Normalized_Radiance' elif ds_id['name'] in self.bt_bands: ds_get_name = 'Brightness_Temperature' else: raise NameError("Datset type "+ds_id['name']+" is not supported.") return ds_get_name
[docs] def get_dataset(self, ds_id, ds_info): """Get the dataset.""" ds_get_name = self._get_dsname(ds_id) data = self[ds_get_name] data.attrs = self.get_metadata(data, ds_info) fill = data.attrs.pop('_FillValue') offset = data.attrs.get('add_offset') scale_factor = data.attrs.get('scale_factor') data = data.where(data != fill) data = data.astype(np.float32) if scale_factor is not None and offset is not None: data = data * scale_factor data = data + offset # Now we correct range from 0-1 to 0-100 for VIS: if ds_id['name'] in self.ref_bands: data = data * 100. return data
[docs] def get_area_def(self, ds_id): """Get the area def.""" ds_get_name = self._get_dsname(ds_id) ds_shape = self[ds_get_name + '/shape'] geoloc = self.geoloc pdict = {} pdict['cfac'] = np.int32(geoloc[0]) pdict['lfac'] = np.int32(geoloc[1]) pdict['coff'] = np.float32(geoloc[2]) pdict['loff'] = -np.float32(geoloc[3]) # Unfortunately this dataset does not store a, b or h. # We assume a and b here, and calculate h from altitude # a and b are from SEVIRI data HRIT header (201912101300) pdict['a'] = 6378169 pdict['b'] = 6356583.8 pdict['h'] = self.alt - pdict['a'] pdict['ssp_lon'] = self.projlon pdict['ncols'] = int(ds_shape[0]) pdict['nlines'] = int(ds_shape[1]) # Force scandir to SEVIRI default, not known from file pdict['scandir'] = 'S2N' pdict['a_name'] = 'geosmsg' if ds_id['name'] == 'HRV': pdict['a_desc'] = 'MSG/SEVIRI HRV channel area' pdict['p_id'] = 'msg_hires' else: pdict['a_desc'] = 'MSG/SEVIRI low resolution channel area' pdict['p_id'] = 'msg_lowres' aex = get_area_extent(pdict) area = get_area_definition(pdict, aex) return area