Implementing a reader for satpy

In order to add a reader to satpy, you will need to implement two files:

• a YAML file for describing the files to read and the datasets that are available
• a python file implementing the actual reading of the datasets and metadata

For this tutorial, we will use the already implemented NWCSAF reader as an example

The YAML file

reader:
  description: NetCDF4 reader for the NWCSAF MSG 2016 format
  name: nc_nwcsaf_msg
  sensors: [seviri]
  default_channels: []
  reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader

file_types:
    nc_nwcsaf_cma:
        file_reader: !!python/name:satpy.readers.nc_nwcsaf.NcNWCSAF
        file_patterns: ['S_NWC_CMA_{platform_id}_{region_id}_{start_time:%Y%m%dT%H%M%S}Z.nc']

    nc_nwcsaf_ct:
        file_reader: !!python/name:satpy.readers.nc_nwcsaf.NcNWCSAF
        file_patterns: ['S_NWC_CT_{platform_id}_{region_id}_{start_time:%Y%m%dT%H%M%S}Z.nc']

    nc_nwcsaf_ctth:
        file_reader: !!python/name:satpy.readers.nc_nwcsaf.NcNWCSAF
        file_patterns: ['S_NWC_CTTH_{platform_id}_{region_id}_{start_time:%Y%m%dT%H%M%S}Z.nc']

    nc_nwcsaf_cmic:
        file_reader: !!python/name:satpy.readers.nc_nwcsaf.NcNWCSAF
        file_patterns: ['S_NWC_CMIC_{platform_id}_{region_id}_{start_time:%Y%m%dT%H%M%S}Z.nc']

    nc_nwcsaf_pc:
        file_reader: !!python/name:satpy.readers.nc_nwcsaf.NcNWCSAF
        file_patterns: ['S_NWC_PC_{platform_id}_{region_id}_{start_time:%Y%m%dT%H%M%S}Z.nc']

    nc_nwcsaf_crr:
        file_reader: !!python/name:satpy.readers.nc_nwcsaf.NcNWCSAF
        file_patterns: ['S_NWC_CRR_{platform_id}_{region_id}_{start_time:%Y%m%dT%H%M%S}Z.nc']

    nc_nwcsaf_ishai:
        file_reader: !!python/name:satpy.readers.nc_nwcsaf.NcNWCSAF
        file_patterns: ['S_NWC_iSHAI_{platform_id}_{region_id}_{start_time:%Y%m%dT%H%M%S}Z.nc']

    nc_nwcsaf_ci:
        file_reader: !!python/name:satpy.readers.nc_nwcsaf.NcNWCSAF
        file_patterns: ['S_NWC_CI_{platform_id}_{region_id}_{start_time:%Y%m%dT%H%M%S}Z.nc']

datasets:

# ---- CMA products ------------

  cma:
    name: cma
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cma

  cma_pal:
    name: cma_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cma

  cma_cloudsnow:
    name: cma_cloudsnow
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cma

  cma_cloudsnow_pal:
    name: cma_cloudsnow_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cma

  cma_dust:
    name: cma_dust
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cma

  cma_dust_pal:
    name: cma_dust_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cma

  cma_volcanic:
    name: cma_volcanic
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cma

  cma_volcanic_pal:
    name: cma_volcanic_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cma

# ---- CT products ------------

  ct:
    name: ct
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ct

  ct_pal:
    name: ct_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ct

# ---- CTTH products ------------

  ctth_alti:
    name: ctth_alti
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ctth

  ctth_alti_pal:
    name: ctth_alti_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ctth

  ctth_pres:
    name: ctth_pres
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ctth

  ctth_pres_pal:
    name: ctth_pres_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ctth

  ctth_tempe:
    name: ctth_tempe
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ctth

  ctth_tempe_pal:
    name: ctth_tempe_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ctth

  ctth_effectiv:
    name: ctth_effectiv
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ctth

  ctth_effectiv_pal:
    name: ctth_effectiv_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ctth

# ---- CMIC products ------------

  cmic_phase:
    name: cmic_phase
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cmic

  cmic_phase_pal:
    name: cmic_phase_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cmic

  cmic_reff:
    name: cmic_reff
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cmic

  cmic_reff_pal:
    name: cmic_reff_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cmic

  cmic_cot:
    name: cmic_cot
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cmic

  cmic_cot_pal:
    name: cmic_cot_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cmic

  cmic_lwp:
    name: cmic_lwp
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cmic

  cmic_lwp_pal:
    name: cmic_lwp_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cmic

  cmic_iwp:
    name: cmic_iwp
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cmic

  cmic_iwp_pal:
    name: cmic_iwp_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_cmic

# ---- PC products ------------

  pc:
    name: pc
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_pc

  pc_pal:
    name: pc_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_pc

# ---- CRR products ------------

  crr:
    name: crr
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_crr

  crr_pal:
    name: crr_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_crr

  crr_intensity:
    name: crr_intensity
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_crr

  crr_intensity_pal:
    name: crr_intensity_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_crr

# ----iSHAI  products ------------

  ishai_tpw:
    name: ishai_tpw
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_tpw_pal:
    name: ishai_tpw_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_shw:
    name: ishai_shw
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_shw_pal:
    name: ishai_shw_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_li:
    name: ishai_li
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_li_pal:
    name: ishai_li_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_ki:
    name: ishai_ki
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_ki_pal:
    name: ishai_ki_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_shw:
    name: ishai_shw
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_shw_pal:
    name: ishai_shw_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_bl:
    name: ishai_bl
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_bl_pal:
    name: ishai_bl_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_ml:
    name: ishai_ml
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_ml_pal:
    name: ishai_ml_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_hl:
    name: ishai_hl
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_hl_pal:
    name: ishai_hl_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_toz:
    name: ishai_toz
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_toz_pal:
    name: ishai_toz_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_skt:
    name: ishai_skt
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai

  ishai_skt_pal:
    name: ishai_skt_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ishai


# ----CI products ------------

  ci_prob30:
    name: ci_prob30
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ci

  ci_prob60:
    name: ci_prob60
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ci

  ci_prob90:
    name: ci_prob90
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ci

  ci_pal:
    name: ci_pal
    sensor: seviri
    resolution: 3000
    file_type: nc_nwcsaf_ci

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# Copyright (c) 2017, 2018 Pytroll

# Author(s):

#   Martin Raspaud <[email protected]>
#   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/>.

"""Nowcasting SAF common PPS&MSG NetCDF/CF format reader
"""

import logging
from datetime import datetime
import os

import numpy as np
import xarray as xr

from pyresample.utils import get_area_def
from satpy.readers.file_handlers import BaseFileHandler
from satpy.readers.utils import unzip_file
from satpy import CHUNK_SIZE

logger = logging.getLogger(__name__)

SENSOR = {'NOAA-19': 'avhrr/3',
          'NOAA-18': 'avhrr/3',
          'NOAA-15': 'avhrr/3',
          'Metop-A': 'avhrr/3',
          'Metop-B': 'avhrr/3',
          'Metop-C': 'avhrr/3',
          'EOS-Aqua': 'modis',
          'EOS-Terra': 'modis',
          'Suomi-NPP': 'viirs',
          'NOAA-20': 'viirs',
          'JPSS-1': 'viirs', }

PLATFORM_NAMES = {'MSG1': 'Meteosat-8',
                  'MSG2': 'Meteosat-9',
                  'MSG3': 'Meteosat-10',
                  'MSG4': 'Meteosat-11', }


class NcNWCSAF(BaseFileHandler):

    """NWCSAF PPS&MSG NetCDF reader."""

    def __init__(self, filename, filename_info, filetype_info):
        """Init method."""
        super(NcNWCSAF, self).__init__(filename, filename_info,
                                       filetype_info)

        self._unzipped = unzip_file(filename)
        if self._unzipped:
            filename = self._unzipped

        self.cache = {}
        self.nc = xr.open_dataset(filename,
                                  decode_cf=True,
                                  mask_and_scale=False,
                                  engine='h5netcdf',
                                  chunks=CHUNK_SIZE)

        self.nc = self.nc.rename({'nx': 'x', 'ny': 'y'})
        self.pps = False

        try:
            # MSG:
            sat_id = self.nc.attrs['satellite_identifier']
            try:
                self.platform_name = PLATFORM_NAMES[sat_id]
            except KeyError:
                self.platform_name = PLATFORM_NAMES[sat_id.astype(str)]
        except KeyError:
            # PPS:
            self.platform_name = self.nc.attrs['platform']
            self.pps = True

        self.sensor = SENSOR.get(self.platform_name, 'seviri')

    def remove_timedim(self, var):
        """Remove time dimension from dataset"""
        if self.pps and var.dims[0] == 'time':
            data = var[0, :, :]
            data.attrs = var.attrs
            var = data

        return var

    def get_dataset(self, dsid, info):
        """Load a dataset."""

        dsid_name = dsid.name
        if dsid_name in self.cache:
            logger.debug('Get the data set from cache: %s.', dsid_name)
            return self.cache[dsid_name]

        if dsid_name in ['lon', 'lat'] and dsid_name not in self.nc.keys():
            dsid_name = dsid_name + '_reduced'

        logger.debug('Reading %s.', dsid_name)
        variable = self.nc[dsid_name]
        variable = self.remove_timedim(variable)
        variable = self.scale_dataset(dsid, variable, info)

        if dsid_name.endswith('_reduced'):
            # Get full resolution lon,lat from the reduced (tie points) grid
            self.upsample_geolocation(dsid, info)

            return self.cache[dsid.name]

        return variable

    def scale_dataset(self, dsid, variable, info):
        """Scale the data set, applying the attributes from the netCDF file"""

        variable = remove_empties(variable)
        scale = variable.attrs.get('scale_factor', np.array(1))
        offset = variable.attrs.get('add_offset', np.array(0))
        if np.issubdtype((scale + offset).dtype, np.floating):
            if '_FillValue' in variable.attrs:
                variable = variable.where(
                    variable != variable.attrs['_FillValue'])
                variable.attrs['_FillValue'] = np.nan
            if 'valid_range' in variable.attrs:
                variable = variable.where(
                    variable <= variable.attrs['valid_range'][1])
                variable = variable.where(
                    variable >= variable.attrs['valid_range'][0])
            if 'valid_max' in variable.attrs:
                variable = variable.where(
                    variable <= variable.attrs['valid_max'])
            if 'valid_min' in variable.attrs:
                variable = variable.where(
                    variable >= variable.attrs['valid_min'])
        attrs = variable.attrs
        variable = variable * scale + offset
        variable.attrs = attrs

        variable.attrs.update({'platform_name': self.platform_name,
                               'sensor': self.sensor})

        variable.attrs.setdefault('units', '1')

        ancillary_names = variable.attrs.get('ancillary_variables', '')
        try:
            variable.attrs['ancillary_variables'] = ancillary_names.split()
        except AttributeError:
            pass

        if 'standard_name' in info:
            variable.attrs.setdefault('standard_name', info['standard_name'])

        if self.pps and dsid.name == 'ctth_alti':
            # pps valid range and palette don't match
            variable.attrs['valid_range'] = (0., 9000.)
        if self.pps and dsid.name == 'ctth_alti_pal':
            # pps palette has the nodata color (black) first
            variable = variable[1:, :]

        return variable

    def upsample_geolocation(self, dsid, info):
        """Upsample the geolocation (lon,lat) from the tiepoint grid"""
        from geotiepoints import SatelliteInterpolator
        # Read the fields needed:
        col_indices = self.nc['nx_reduced'].values
        row_indices = self.nc['ny_reduced'].values
        lat_reduced = self.scale_dataset(dsid, self.nc['lat_reduced'], info)
        lon_reduced = self.scale_dataset(dsid, self.nc['lon_reduced'], info)

        shape = (self.nc['y'].shape[0], self.nc['x'].shape[0])
        cols_full = np.arange(shape[1])
        rows_full = np.arange(shape[0])

        satint = SatelliteInterpolator((lon_reduced.values, lat_reduced.values),
                                       (row_indices,
                                        col_indices),
                                       (rows_full, cols_full))

        lons, lats = satint.interpolate()
        self.cache['lon'] = xr.DataArray(lons, attrs=lon_reduced.attrs, dims=['y', 'x'])
        self.cache['lat'] = xr.DataArray(lats, attrs=lat_reduced.attrs, dims=['y', 'x'])

        return

    def get_area_def(self, dsid):
        """Get the area definition of the datasets in the file.

        Only applicable for MSG products!
        """
        if self.pps:
            # PPS:
            raise NotImplementedError

        if dsid.name.endswith('_pal'):
            raise NotImplementedError

        try:
            proj_str = self.nc.attrs['gdal_projection'] + ' +units=km'
        except TypeError:
            proj_str = self.nc.attrs['gdal_projection'].decode() + ' +units=km'

        nlines, ncols = self.nc[dsid.name].shape

        area_extent = (float(self.nc.attrs['gdal_xgeo_up_left']) / 1000,
                       float(self.nc.attrs['gdal_ygeo_low_right']) / 1000,
                       float(self.nc.attrs['gdal_xgeo_low_right']) / 1000,
                       float(self.nc.attrs['gdal_ygeo_up_left']) / 1000)

        area = get_area_def('some_area_name',
                            "On-the-fly area",
                            'geosmsg',
                            proj_str,
                            ncols,
                            nlines,
                            area_extent)

        return area

    def __del__(self):
        if self._unzipped:
            try:
                os.remove(self._unzipped)
            except (IOError, OSError):
                pass

    @property
    def start_time(self):
        """Return the start time of the object."""
        try:
            # MSG:
            try:
                return datetime.strptime(self.nc.attrs['time_coverage_start'],
                                         '%Y-%m-%dT%H:%M:%SZ')
            except TypeError:
                return datetime.strptime(self.nc.attrs['time_coverage_start'].astype(str),
                                         '%Y-%m-%dT%H:%M:%SZ')
        except ValueError:
            # PPS:
            return datetime.strptime(self.nc.attrs['time_coverage_start'],
                                     '%Y%m%dT%H%M%S%fZ')

    @property
    def end_time(self):
        """Return the end time of the object."""
        try:
            # MSG:
            try:
                return datetime.strptime(self.nc.attrs['time_coverage_end'],
                                         '%Y-%m-%dT%H:%M:%SZ')
            except TypeError:
                return datetime.strptime(self.nc.attrs['time_coverage_end'].astype(str),
                                         '%Y-%m-%dT%H:%M:%SZ')
        except ValueError:
            # PPS:
            return datetime.strptime(self.nc.attrs['time_coverage_end'],
                                     '%Y%m%dT%H%M%S%fZ')


def remove_empties(variable):
    """Remove empty objects from the *variable*'s attrs."""
    import h5py
    for key, val in variable.attrs.items():
        if isinstance(val, h5py._hl.base.Empty):
            variable.attrs.pop(key)

    return variable