From f8699acba7511038f5133950d8a22390f89a81d6 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 6 Nov 2020 15:52:17 +0000
Subject: [PATCH 01/47] Add reader for FIDUCEO MVIRI FCDR data

---
 doc/source/api/satpy.readers.rst            |   8 +
 satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml | 101 ++++
 satpy/readers/_geos_area.py                 |  21 +
 satpy/readers/mviri_l1b_fiduceo_nc.py       | 523 ++++++++++++++++++++
 4 files changed, 653 insertions(+)
 create mode 100644 satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
 create mode 100644 satpy/readers/mviri_l1b_fiduceo_nc.py

diff --git a/doc/source/api/satpy.readers.rst b/doc/source/api/satpy.readers.rst
index ad79cd7dd4..cdbb7988cc 100644
--- a/doc/source/api/satpy.readers.rst
+++ b/doc/source/api/satpy.readers.rst
@@ -364,6 +364,14 @@ satpy.readers.msi\_safe module
    :undoc-members:
    :show-inheritance:
 
+satpy.readers.mviri\_l1b\_fiduceo\_nc module
+------------------------------
+
+.. automodule:: satpy.readers.mviri_l1b_fiduceo_nc
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 satpy.readers.netcdf\_utils module
 ----------------------------------
 
diff --git a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
new file mode 100644
index 0000000000..736307f7ba
--- /dev/null
+++ b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
@@ -0,0 +1,101 @@
+# References:
+#    - MFG User Handbook
+#    - FIDUCEO MVIRI FCDR Product User Guide
+
+reader:
+  name: mviri_l1b_fiduceo_nc
+  short_name: FIDUCEO MVIRI FCDR
+  long_name: >
+    Fundamental Climate Data Record of re-calibrated Level 1.5 Infrared, Water Vapour, and Visible radiances
+    from the Meteosat Visible Infra-Red Imager (MVIRI) instrument onboard the Meteosat First Generation satellites
+  description: >
+    Reader for FIDUCEO MVIRI FCDR data in netCDF format. For documentation see
+    http://doi.org/10.15770/EUM_SEC_CLM_0009 .
+  sensors: [mviri]
+  default_channels: [VIS, WV, IR]
+  reader: !!python/name:satpy.readers.yaml_reader.GEOFlippableFileYAMLReader
+
+file_types:
+    nc_easy:
+        file_reader: !!python/name:satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriEasyFcdrFileHandler
+        file_patterns: [
+          'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_EASY_{processor_version}_{format_version}.nc'
+        ]
+    nc_full:
+        file_reader: !!python/name:satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriFullFcdrFileHandler
+        file_patterns: [
+          'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_FULL_{processor_version}_{format_version}.nc'
+        ]
+
+datasets:
+  VIS:
+    name: VIS
+    resolution: 2250
+    wavelength: [0.5, 0.7, 0.9]
+    calibration:
+      reflectance:
+        standard_name: toa_bidirectional_reflectance
+        units: "%"
+      radiance:
+        # standard_name: TODO
+        units: W m-2 sr-1  # TODO: Per unit wave number/length
+      counts:
+        standard_name: counts
+        units: count
+    file_type: [nc_easy, nc_full]
+
+  WV:
+    name: WV
+    resolution: 4500
+    wavelength: [5.7, 6.4, 7.1]
+    calibration:
+      brightness_temperature:
+        standard_name: toa_brightness_temperature
+        units: K
+      radiance:
+        standard_name: toa_outgoing_radiance_per_unit_wavelength
+        units: mW m-2 sr-1 (cm-1)-1
+      counts:
+        standard_name: counts
+        units: count
+    file_type: [nc_easy, nc_full]
+
+  IR:
+    name: IR
+    resolution: 4500
+    wavelength: [10.5, 11.5, 12.5]
+    calibration:
+      brightness_temperature:
+        standard_name: toa_brightness_temperature
+        units: K
+      radiance:
+        standard_name: toa_outgoing_radiance_per_unit_wavelength
+        units: mW m-2 sr-1 (cm-1)-1
+      counts:
+        standard_name: counts
+        units: count
+    file_type: [nc_easy, nc_full]
+
+  quality_pixel_bitmask:
+    name: quality_pixel_bitmask
+    resolution: 2250
+    file_type: [nc_easy, nc_full]
+
+  data_quality_bitmask:
+    name: data_quality_bitmask
+    resolution: 2250
+    file_type: [nc_easy, nc_full]
+
+  u_independent_toa_bidirectional_reflectance:
+    name: u_independent_toa_bidirectional_reflectance
+    long_name: "independent uncertainty per pixel"
+    units: "%"
+    resolution: 2250
+    file_type: [nc_easy]
+
+  u_structured_toa_bidirectional_reflectance:
+    name: u_structured_toa_bidirectional_reflectance
+    long_name: "structured uncertainty per pixel"
+    units: "%"
+    resolution: 2250
+    file_type: [nc_easy]
diff --git a/satpy/readers/_geos_area.py b/satpy/readers/_geos_area.py
index b340f5fca0..d6632ab8f6 100644
--- a/satpy/readers/_geos_area.py
+++ b/satpy/readers/_geos_area.py
@@ -145,3 +145,24 @@ def get_area_definition(pdict, a_ext):
         a_ext)
 
     return a_def
+
+
+def ang2fac(ang):
+    """Convert angular sampling to line/column scaling factor (aka LFAC/CFAC).
+
+    Reference: `MSG Ground Segment LRIT HRIT Mission Specific Implementation`_,
+    Appendix E.2.
+
+
+    .. _MSG Ground Segment LRIT HRIT Mission Specific Implementation: http://www.eumetsat.int/\
+website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_TEN_05057_SPE_MSG_LRIT_HRI\
+&RevisionSelectionMethod=LatestReleased&Rendition=Web
+
+    Args:
+        ang: float
+            Angular sampling (rad)
+
+    Returns:
+        Line/column scaling factor (deg-1)
+    """
+    return 2.0 ** 16 / np.rad2deg(ang)
diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
new file mode 100644
index 0000000000..ea75acf782
--- /dev/null
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -0,0 +1,523 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2020 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 <http://www.gnu.org/licenses/>.
+"""FIDUCEO MVIRI FCDR Reader.
+
+Introduction
+------------
+The FIDUCEO MVIRI FCDR is a Fundamental Climate Data Record (FCDR) of
+re-calibrated Level 1.5 Infrared, Water Vapour, and Visible radiances from
+the Meteosat Visible Infra-Red Imager (MVIRI) instrument onboard the
+Meteosat First Generation satellites. There are two variants of the dataset:
+The *full FCDR* and a simplified version called *easy FCDR*. Some datasets are
+only available in one of the two variants, see the corresponding YAML
+definition in ``satpy/etc/readers/``.
+
+Dataset Names
+-------------
+The FIDUCEO MVIRI readers use names ``VIS``, ``WV`` and ``IR`` for the visible,
+water vapor and infrared channels, respectively. These are different from
+the original netCDF variable names for the following reasons:
+
+- VIS channel is named differently in full FCDR (``counts_vis``) and easy FCDR
+  (``toa_bidirectional_reflectance_vis``)
+- netCDF variable names contain the calibration level (e.g. ``counts_...``),
+  which might be confusing for satpy users if a different calibration level
+  is chosen.
+
+Remaining datasets (such as quality flags and uncertainties) have the same
+name in the reader as in the netCDF file.
+
+
+Example
+-------
+This is how to read FIDUCEO MVIRI FCDR data in satpy:
+
+.. code-block:: python
+
+    from satpy import Scene
+
+    scn = Scene(filenames=['FIDUCEO_FCDR_L15_MVIRI_MET7-57.0...'],
+                reader='mviri_l1b_fiduceo_nc',
+                reader_kwargs={'mask_bad_quality': True)
+    scn.load(['VIS', 'WV', 'IR'])
+
+In the above example pixels considered bad quality are masked, see
+:class:`FiduceoMviriBase` for a keyword argument description. Global
+netCDF attributes are available in the ``raw_metadata`` attribute of
+each loaded dataset.
+
+Geolocation
+-----------
+In addition to the image data, FIDUCEO also provides so called *static FCDRs*
+containing latitude and longitude coordinates. In order to simplify their
+usage, the FIDUCEO MVIRI readers do not make use of these static files, but
+instead provide an area definition that can be used to compute longitude and
+latitude coordinates on demand.
+
+.. code-block:: python
+
+    area = scn['VIS'].attrs['area']
+    lons, lats = area.get_lonlats()
+
+Those were compared to the static FCDR and they agree very well, however there
+are small differences. The mean difference is < 1E3 degrees for all channels
+and projection longitudes.
+
+
+References:
+    - `[Handbook]`_ MFG User Handbook
+    - `[PUG]`_ FIDUCEO MVIRI FCDR Product User Guide
+
+.. _[Handbook]: http://www.eumetsat.int/\
+website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_TD06_MARF&\
+RevisionSelectionMethod=LatestReleased&Rendition=Web
+.. _[PUG]: http://doi.org/10.15770/EUM_SEC_CLM_0009
+"""
+
+import abc
+import functools
+import numpy as np
+import xarray as xr
+
+
+from satpy import CHUNK_SIZE
+from satpy.readers.file_handlers import BaseFileHandler
+from satpy.readers._geos_area import (ang2fac, get_area_definition,
+                                      get_area_extent)
+
+
+EQUATOR_RADIUS = 6378140.0
+POLE_RADIUS = 6356755.0
+ALTITUDE = 42164000.0 - EQUATOR_RADIUS
+"""[Handbook] section 5.2.1."""
+
+MVIRI_FIELD_OF_VIEW = 18.0
+"""[Handbook] section 5.3.2.1."""
+
+CHANNELS = ['VIS', 'WV', 'IR']
+OTHER_REFLECTANCES = [
+    'u_independent_toa_bidirectional_reflectance',
+    'u_structured_toa_bidirectional_reflectance'
+]
+
+
+class shape_cache:
+    """Cache function call based on image shape.
+
+    Shape can be used to maintain separate caches for low resolution
+    (WV/IR) and high resolution (VIS) channels.
+    """
+
+    def __init__(self, func):
+        """Create the cache.
+
+        Args:
+            func:
+                Function to me cached.
+        """
+        self.func = func
+        self.cache = {}
+
+    def __call__(self, *args):
+        """Call the decorated function.
+
+        Dataset is expected to be the last argument. If an instance method is
+        decorated, it is preceded by a filehandler instance.
+        """
+        ds = args[-1]
+        shape = ds.coords['y'].shape
+        if shape not in self.cache:
+            self.cache[shape] = self.func(*args)
+        return self.cache[shape]
+
+    def __get__(self, obj, objtype):
+        """To support instance methods."""
+        return functools.partial(self.__call__, obj)
+
+
+class FiduceoMviriBase(BaseFileHandler):
+    """Baseclass for FIDUCEO MVIRI file handlers."""
+
+    def __init__(self, filename, filename_info, filetype_info,
+                 mask_bad_quality=False):
+        """Initialize the file handler.
+
+        Args:
+             mask_bad_quality: Mask pixels with bad quality, that means
+                 everything else than "ok" or "use with caution". If you need
+                 more control, use the ``quality_pixel_bitmask`` and
+                 ``data_quality_bitmask`` datasets.
+        """
+        super(FiduceoMviriBase, self).__init__(
+            filename, filename_info, filetype_info)
+        self.mask_bad_quality = mask_bad_quality
+        self.nc = xr.open_dataset(filename, chunks=CHUNK_SIZE)
+
+        # Projection longitude is not provided in the file, read it from the
+        # filename.
+        self.projection_longitude = filename_info['projection_longitude']
+
+    @abc.abstractproperty
+    def nc_keys(self):
+        """Map satpy dataset names to netCDF variables."""
+        raise NotImplementedError
+
+    def get_dataset(self, dataset_id, info):
+        """Get the dataset."""
+        name = dataset_id['name']
+        ds = self._read_dataset(name)
+        if dataset_id['name'] in CHANNELS:
+            ds = self.calibrate(ds, channel=name,
+                                calibration=dataset_id['calibration'])
+            if self.mask_bad_quality:
+                ds = self._mask(ds)
+            ds['acq_time'] = ('y', self._get_acq_time(ds))
+        elif dataset_id['name'] in OTHER_REFLECTANCES:
+            ds = ds * 100  # conversion to percent
+        self._update_attrs(ds)
+        ds.attrs['orbital_parameters'] = self._get_orbital_parameters()
+        return ds
+
+    def _get_nc_key(self, name):
+        """Get netCDF key corresponding to the given dataset name."""
+        return
+
+    def _read_dataset(self, name):
+        """Read a dataset from the file."""
+        nc_key = self.nc_keys.get(name, name)
+        ds = self.nc[nc_key]
+        if 'y_ir_wv' in ds.dims:
+            ds = ds.rename({'y_ir_wv': 'y', 'x_ir_wv': 'x'})
+        return ds
+
+    def _update_attrs(self, ds):
+        """Update dataset attributes."""
+        ds.attrs.update({'platform': self.filename_info['platform'],
+                         'sensor': self.filename_info['sensor']})
+        ds.attrs['raw_metadata'] = self.nc.attrs
+
+    def get_area_def(self, dataset_id):
+        """Get area definition of the given dataset."""
+        ds = self._read_dataset(dataset_id['name'])
+        im_size = ds.coords['y'].size
+
+        # Determine line/column offsets and scaling factors. For offsets
+        # see variables "asamp" and "aline" of subroutine "refgeo" in
+        # [Handbook] and in
+        # https://github.com/FIDUCEO/FCDR_MVIRI/blob/master/lib/nrCrunch/cruncher.f
+        loff = coff = im_size / 2 + 0.5
+        lfac = cfac = ang2fac(np.deg2rad(MVIRI_FIELD_OF_VIEW) / im_size)
+
+        pdict = {
+            'ssp_lon': self.projection_longitude,
+            'a': EQUATOR_RADIUS,
+            'b': POLE_RADIUS,
+            'h': ALTITUDE,
+            'units': 'm',
+            'loff': loff - im_size,
+            'coff': coff,
+            'lfac': -lfac,
+            'cfac': -cfac,
+            'nlines': im_size,
+            'ncols': im_size,
+            'scandir': 'S2N',  # Reference: [PUG] section 2.
+            'p_id': 'geos_mviri',
+            'a_name': 'geos_mviri',
+            'a_desc': 'MVIRI Geostationary Projection'
+        }
+        extent = get_area_extent(pdict)
+        area_def = get_area_definition(pdict, extent)
+        return area_def
+
+    def calibrate(self, ds, channel, calibration):
+        """Calibrate the given dataset."""
+        if channel == 'VIS':
+            return self._calibrate_vis(ds, calibration)
+        elif channel in ['WV', 'IR']:
+            return self._calibrate_ir_wv(ds, channel, calibration)
+        else:
+            raise KeyError('Don\'t know how to calibrate channel {}'.format(
+                channel))
+
+    @abc.abstractmethod
+    def _calibrate_vis(self, ds, calibration):
+        """Calibrate VIS channel. To be implemented by subclasses."""
+        raise NotImplementedError
+
+    def _update_refl_attrs(self, refl):
+        """Update attributes of reflectance datasets."""
+        refl.attrs['sun_earth_distance_correction_applied'] = True
+        refl.attrs['sun_earth_distance_correction_factor'] = self.nc[
+            'distance_sun_earth'].values
+        return refl
+
+    def _calibrate_ir_wv(self, ds, channel, calibration):
+        """Calibrate IR and WV channel."""
+        if calibration == 'counts':
+            return ds
+        elif calibration in ('radiance', 'brightness_temperature'):
+            rad = self._ir_wv_counts_to_radiance(ds, channel)
+            if calibration == 'radiance':
+                return rad
+            bt = self._ir_wv_radiance_to_brightness_temperature(rad, channel)
+            return bt
+
+    def _ir_wv_counts_to_radiance(self, counts, channel):
+        """Convert IR/WV counts to radiance.
+
+        Reference: [PUG], equations (4.1) and (4.2).
+        """
+        if channel == 'WV':
+            a, b = self.nc['a_wv'], self.nc['b_wv']
+        else:
+            a, b = self.nc['a_ir'], self.nc['b_ir']
+        rad = a + b * counts
+        return rad.where(rad > 0)
+
+    def _ir_wv_radiance_to_brightness_temperature(self, rad, channel):
+        """Convert IR/WV radiance to brightness temperature."""
+        if channel == 'WV':
+            a, b = self.nc['bt_a_wv'], self.nc['bt_b_wv']
+        else:
+            a, b = self.nc['bt_a_ir'], self.nc['bt_b_ir']
+        bt = b / (np.log(rad) - a)
+        return bt.where(bt > 0)
+
+    def _mask(self, ds):
+        """Mask pixels with bad quality.
+
+        Pixels are considered bad quality if the "quality_pixel_bitmask" is
+        everything else than 0 (no flag set) or 2 ("use_with_caution" and no
+        other flag set).
+        """
+        mask = self._get_mask(ds)
+        return ds.where(np.logical_or(mask == 0, mask == 2))
+
+    @shape_cache
+    def _get_mask(self, ds):
+        """Get quality bitmask for the given dataset."""
+        mask = self.nc['quality_pixel_bitmask']
+
+        # Mask has high (VIS) resolution. Downsample to low (IR/WV) resolution
+        # if required.
+        if mask.coords['y'].size > ds.coords['y'].size:
+            mask = mask.isel(y=slice(None, None, 2),
+                             x=slice(None, None, 2))
+            mask = mask.assign_coords(y=ds.coords['y'].values,
+                                      x=ds.coords['x'].values)
+        return mask
+
+    @shape_cache
+    def _get_acq_time(self, ds):
+        """Get scanline acquisition time.
+
+        The files provide timestamps per pixel for the low resolution
+        channels (IR/WV) only.
+
+        1) Average values in each line to obtain one timestamp per line.
+        2) For the VIS channel duplicate values in y-direction (as
+           advised by [PUG]).
+
+        Note that the timestamps do not increase monotonically with the
+        line number in some cases.
+
+        Returns:
+            Mean scanline acquisition timestamps
+        """
+        # Variable is sometimes named "time" and sometimes "time_ir_wv".
+        try:
+            time_lores = self.nc['time_ir_wv']
+        except KeyError:
+            time_lores = self.nc['time']
+
+        # Compute mean timestamp per scanline
+        time_lores = time_lores.mean(dim='x_ir_wv').rename({'y_ir_wv': 'y'})
+
+        # If required, repeat timestamps in y-direction to obtain higher
+        # resolution
+        y_lores = time_lores.coords['y'].values
+        y_target = ds.coords['y'].values
+        if y_lores.size < y_target.size:
+            reps = y_target.size // y_lores.size
+            y_lores_rep = np.repeat(y_lores, reps)
+            time_hires = time_lores.reindex(y=y_lores_rep)
+            time_hires = time_hires.assign_coords(y=y_target)
+            return time_hires
+
+        return time_lores
+
+    def _get_orbital_parameters(self):
+        """Get the orbital parameters."""
+        ssp_lon, ssp_lat = self._get_ssp_lonlat()
+        orbital_parameters = {
+            'projection_longitude': self.projection_longitude,
+            'projection_latitude': 0.0,
+            'projection_altitude': ALTITUDE
+        }
+        if not np.isnan(ssp_lon) and not np.isnan(ssp_lat):
+            orbital_parameters.update({
+                'satellite_actual_longitude': ssp_lon,
+                'satellite_actual_latitude': ssp_lat,
+                # altitude not available
+            })
+
+    def _get_ssp_lonlat(self):
+        """Get longitude and latitude at the subsatellite point.
+
+        Easy FCDR files provide satellite position at the beginning and
+        end of the scan. This method computes the mean of those two values.
+        In the full FCDR the information seems to be missing.
+
+        Returns:
+            Subsatellite longitude and latitude
+        """
+        ssp_lon = self._get_ssp('longitude')
+        ssp_lat = self._get_ssp('latitude')
+        return ssp_lon, ssp_lat
+
+    def _get_ssp(self, coord):
+        key_start = 'sub_satellite_{}_start'.format(coord)
+        key_end = 'sub_satellite_{}_end'.format(coord)
+        try:
+            sub_lonlat = np.nanmean(
+                [self.nc[key_start].values,
+                 self.nc[key_end].values]
+            )
+        except KeyError:
+            # Variables seem to be missing in full FCDR
+            sub_lonlat = np.nan
+        return sub_lonlat
+
+    @shape_cache
+    def _get_solar_zenith_angle(self, ds):
+        """Get solar zenith angle for the given dataset.
+
+        Files provide solar zenith angle, but at a coarser resolution.
+        Interpolate to the resolution of the given dataset.
+        """
+        return self._interp_tiepoints(self.nc['solar_zenith_angle'],
+                                      ds.coords['x'],
+                                      ds.coords['y'])
+
+    def _interp_tiepoints(self, ds, target_x, target_y):
+        """Interpolate dataset between tiepoints.
+
+        Uses linear interpolation.
+
+        FUTURE: [PUG] recommends cubic spline interpolation.
+
+        Args:
+            ds:
+                Dataset to be interpolated
+            target_x:
+                Target x coordinates
+            target_y:
+                Target y coordinates
+        """
+        if 'x_tie' not in ds.dims:
+            raise ValueError('Dataset has no tiepoints to be interpolated.')
+
+        # No tiepoint coordinates specified in the files. Use dimensions
+        # to calculate tiepoint sampling and assign tiepoint coordinates
+        # accordingly.
+        sampling = target_x.size // ds.coords['x_tie'].size
+        ds = ds.assign_coords(x_tie=target_x.values[::sampling],
+                              y_tie=target_y.values[::sampling])
+
+        ds_interp = ds.interp(x_tie=target_x, y_tie=target_y)
+        return ds_interp.rename({'x_tie': 'x', 'y_tie': 'y'})
+
+
+class FiduceoMviriEasyFcdrFileHandler(FiduceoMviriBase):
+    """File handler for FIDUCEO MVIRI Easy FCDR."""
+
+    nc_keys = {
+        'VIS': 'toa_bidirectional_reflectance_vis',
+        'WV': 'count_wv',
+        'IR': 'count_ir'
+    }
+
+    def _calibrate_vis(self, ds, calibration):
+        """Calibrate VIS channel.
+
+        Easy FCDR provides reflectance only, no counts or radiance.
+        """
+        if calibration == 'reflectance':
+            refl = 100 * ds  # conversion to percent
+            refl = self._update_refl_attrs(refl)
+            return refl
+        elif calibration in ('counts', 'radiance'):
+            raise ValueError('Cannot calibrate to {}. Easy FCDR provides '
+                             'reflectance only.'.format(calibration.name))
+        else:
+            raise KeyError('Invalid calibration: {}'.format(calibration.name))
+
+
+class FiduceoMviriFullFcdrFileHandler(FiduceoMviriBase):
+    """File handler for FIDUCEO MVIRI Full FCDR."""
+
+    nc_keys = {
+        'VIS': 'count_vis',
+        'WV': 'count_wv',
+        'IR': 'count_ir'
+    }
+
+    def _calibrate_vis(self, ds, calibration):
+        """Calibrate VIS channel.
+
+        All calibration levels are available here.
+        """
+        if calibration == 'counts':
+            return ds
+        elif calibration in ('radiance', 'reflectance'):
+            rad = self._vis_counts_to_radiance(ds)
+            if calibration == 'radiance':
+                return rad
+            refl = self._vis_radiance_to_reflectance(rad)
+            refl = self._update_refl_attrs(refl)
+            return refl
+        else:
+            raise KeyError('Invalid calibration: {}'.format(calibration.name))
+
+    def _vis_counts_to_radiance(self, counts):
+        """Convert VIS counts to radiance.
+
+        Reference: [PUG], equations (7) and (8).
+        """
+        years_since_launch = self.nc['years_since_launch']
+        a_cf = (self.nc['a0_vis'] +
+                self.nc['a1_vis'] * years_since_launch +
+                self.nc['a2_vis'] * years_since_launch ** 2)
+
+        rad = (counts - self.nc['mean_count_space_vis']) * a_cf
+        return rad.where(rad > 0)
+
+    def _vis_radiance_to_reflectance(self, rad):
+        """Convert VIS radiance to reflectance factor.
+
+        Reference: [PUG], equation (6).
+        """
+        sza = self._get_solar_zenith_angle(rad)
+        cos_sza = np.cos(np.deg2rad(sza))
+        refl = (
+           (np.pi * self.nc['distance_sun_earth'] ** 2) /
+           (self.nc['solar_irradiance_vis'] * cos_sza) *
+           rad
+        )
+        refl = refl * 100  # conversion to percent
+        return refl.where(refl > 0)

From f04b3bd6cf5ef81c8dfc76381e29428adb99fee1 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 6 Nov 2020 16:23:27 +0000
Subject: [PATCH 02/47] Add example filenames to YAML

---
 satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
index 736307f7ba..5ddacf574c 100644
--- a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
+++ b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
@@ -20,11 +20,13 @@ file_types:
         file_reader: !!python/name:satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriEasyFcdrFileHandler
         file_patterns: [
           'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_EASY_{processor_version}_{format_version}.nc'
+          # Example: FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc
         ]
     nc_full:
         file_reader: !!python/name:satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriFullFcdrFileHandler
         file_patterns: [
           'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_FULL_{processor_version}_{format_version}.nc'
+          # Example: FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_FULL_v2.6_fv3.1.nc
         ]
 
 datasets:

From b784bdf1581abfb67d3ed394c96d98b2edd69203 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 6 Nov 2020 16:23:49 +0000
Subject: [PATCH 03/47] Add dataset info to dataset attributes

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index ea75acf782..b8ee463856 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -166,7 +166,13 @@ def __init__(self, filename, filename_info, filetype_info,
         super(FiduceoMviriBase, self).__init__(
             filename, filename_info, filetype_info)
         self.mask_bad_quality = mask_bad_quality
-        self.nc = xr.open_dataset(filename, chunks=CHUNK_SIZE)
+        self.nc = xr.open_dataset(
+            filename,
+            chunks={'x': CHUNK_SIZE,
+                    'y': CHUNK_SIZE,
+                    'x_ir_wv': CHUNK_SIZE,
+                    'y_ir_wv': CHUNK_SIZE}
+        )
 
         # Projection longitude is not provided in the file, read it from the
         # filename.
@@ -189,8 +195,7 @@ def get_dataset(self, dataset_id, info):
             ds['acq_time'] = ('y', self._get_acq_time(ds))
         elif dataset_id['name'] in OTHER_REFLECTANCES:
             ds = ds * 100  # conversion to percent
-        self._update_attrs(ds)
-        ds.attrs['orbital_parameters'] = self._get_orbital_parameters()
+        self._update_attrs(ds, info)
         return ds
 
     def _get_nc_key(self, name):
@@ -205,11 +210,13 @@ def _read_dataset(self, name):
             ds = ds.rename({'y_ir_wv': 'y', 'x_ir_wv': 'x'})
         return ds
 
-    def _update_attrs(self, ds):
+    def _update_attrs(self, ds, info):
         """Update dataset attributes."""
+        ds.attrs.update(info)
         ds.attrs.update({'platform': self.filename_info['platform'],
                          'sensor': self.filename_info['sensor']})
         ds.attrs['raw_metadata'] = self.nc.attrs
+        ds.attrs['orbital_parameters'] = self._get_orbital_parameters()
 
     def get_area_def(self, dataset_id):
         """Get area definition of the given dataset."""

From 061d58c8fcd53cfde69b406e9c20fa3dbc274ab4 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 6 Nov 2020 16:36:36 +0000
Subject: [PATCH 04/47] Update documentation

---
 doc/source/api/satpy.readers.rst      |  2 +-
 satpy/readers/mviri_l1b_fiduceo_nc.py | 14 +++++++++++++-
 2 files changed, 14 insertions(+), 2 deletions(-)

diff --git a/doc/source/api/satpy.readers.rst b/doc/source/api/satpy.readers.rst
index cdbb7988cc..06810f740b 100644
--- a/doc/source/api/satpy.readers.rst
+++ b/doc/source/api/satpy.readers.rst
@@ -365,7 +365,7 @@ satpy.readers.msi\_safe module
    :show-inheritance:
 
 satpy.readers.mviri\_l1b\_fiduceo\_nc module
-------------------------------
+--------------------------------------------
 
 .. automodule:: satpy.readers.mviri_l1b_fiduceo_nc
    :members:
diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index b8ee463856..1a5a1e4d78 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -61,6 +61,17 @@
 netCDF attributes are available in the ``raw_metadata`` attribute of
 each loaded dataset.
 
+
+Image Orientation
+-----------------
+The images are stored in MVIRI scanning direction, that means South is up and
+East is right. This can be changed as follows:
+
+.. code-block:: python
+
+    scn.load(['VIS'], upper_right_corner='NE')
+
+
 Geolocation
 -----------
 In addition to the image data, FIDUCEO also provides so called *static FCDRs*
@@ -79,7 +90,8 @@
 and projection longitudes.
 
 
-References:
+References
+----------
     - `[Handbook]`_ MFG User Handbook
     - `[PUG]`_ FIDUCEO MVIRI FCDR Product User Guide
 

From b8a0c2aced2ff34dd0741aecf747d992226c045a Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 6 Nov 2020 16:43:48 +0000
Subject: [PATCH 05/47] Add calibration reference

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 1a5a1e4d78..896efb5242 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -309,7 +309,10 @@ def _ir_wv_counts_to_radiance(self, counts, channel):
         return rad.where(rad > 0)
 
     def _ir_wv_radiance_to_brightness_temperature(self, rad, channel):
-        """Convert IR/WV radiance to brightness temperature."""
+        """Convert IR/WV radiance to brightness temperature.
+
+        Reference: [PUG], equations (5.1) and (5.2).
+        """
         if channel == 'WV':
             a, b = self.nc['bt_a_wv'], self.nc['bt_b_wv']
         else:

From 0b6e49d4cc70733d2b4d9938b9bb41f66bade7a6 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 6 Nov 2020 16:45:06 +0000
Subject: [PATCH 06/47] Raise KeyError in case of invalid calibration

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 896efb5242..d00f849c28 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -295,6 +295,8 @@ def _calibrate_ir_wv(self, ds, channel, calibration):
                 return rad
             bt = self._ir_wv_radiance_to_brightness_temperature(rad, channel)
             return bt
+        else:
+            raise KeyError('Invalid calibration: {}'.format(calibration.name))
 
     def _ir_wv_counts_to_radiance(self, counts, channel):
         """Convert IR/WV counts to radiance.

From 0640ea36fe0e4b35a6e0628fbdc570ee6540ec02 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 18 Nov 2020 14:45:03 +0000
Subject: [PATCH 07/47] Consider only direct illumination for reflectance

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 18 ++++++++++++++++--
 1 file changed, 16 insertions(+), 2 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index d00f849c28..bbd01ccfcb 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -90,6 +90,14 @@
 and projection longitudes.
 
 
+Huge VIS Reflectances
+---------------------
+You might encounter huge VIS reflectances (10^8 percent and greater) in
+situations where both radiance and solar zenith angle are small. The reader
+certainly needs some improvement in this regard. Maybe the corresponding
+uncertainties can be used to filter these cases before calculating reflectances.
+
+
 References
 ----------
     - `[Handbook]`_ MFG User Handbook
@@ -102,11 +110,11 @@
 """
 
 import abc
+import dask.array as da
 import functools
 import numpy as np
 import xarray as xr
 
-
 from satpy import CHUNK_SIZE
 from satpy.readers.file_handlers import BaseFileHandler
 from satpy.readers._geos_area import (ang2fac, get_area_definition,
@@ -534,9 +542,15 @@ def _vis_counts_to_radiance(self, counts):
     def _vis_radiance_to_reflectance(self, rad):
         """Convert VIS radiance to reflectance factor.
 
+        Note: Produces huge reflectances in situations where both radiance and
+        solar zenith angle are small. Maybe the corresponding uncertainties
+        can be used to filter these cases before calculating reflectances.
+
         Reference: [PUG], equation (6).
         """
+
         sza = self._get_solar_zenith_angle(rad)
+        sza = sza.where(da.fabs(sza) < 90)  # direct illumination only
         cos_sza = np.cos(np.deg2rad(sza))
         refl = (
            (np.pi * self.nc['distance_sun_earth'] ** 2) /
@@ -544,4 +558,4 @@ def _vis_radiance_to_reflectance(self, rad):
            rad
         )
         refl = refl * 100  # conversion to percent
-        return refl.where(refl > 0)
+        return refl

From f375f7f7b3ad78eeafed1af7c93c7049155e973f Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 18 Nov 2020 14:46:45 +0000
Subject: [PATCH 08/47] Make sure projection longitude is a float

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index bbd01ccfcb..c247ff57ce 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -196,12 +196,7 @@ def __init__(self, filename, filename_info, filetype_info,
 
         # Projection longitude is not provided in the file, read it from the
         # filename.
-        self.projection_longitude = filename_info['projection_longitude']
-
-    @abc.abstractproperty
-    def nc_keys(self):
-        """Map satpy dataset names to netCDF variables."""
-        raise NotImplementedError
+        self.projection_longitude = float(filename_info['projection_longitude'])
 
     def get_dataset(self, dataset_id, info):
         """Get the dataset."""

From 13144978525d496b956546d6e297418c82a9f4be Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 18 Nov 2020 14:48:02 +0000
Subject: [PATCH 09/47] Make sure sun earth dist corr factor is scalar

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index c247ff57ce..6fe033fed7 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -285,7 +285,7 @@ def _update_refl_attrs(self, refl):
         """Update attributes of reflectance datasets."""
         refl.attrs['sun_earth_distance_correction_applied'] = True
         refl.attrs['sun_earth_distance_correction_factor'] = self.nc[
-            'distance_sun_earth'].values
+            'distance_sun_earth'].item()
         return refl
 
     def _calibrate_ir_wv(self, ds, channel, calibration):

From bfe6787009035d4e5c635681fa9ab47bbcf67738 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 18 Nov 2020 14:48:40 +0000
Subject: [PATCH 10/47] Fix orbital parameters

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 6fe033fed7..2a145fc32e 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -390,18 +390,19 @@ def _get_acq_time(self, ds):
 
     def _get_orbital_parameters(self):
         """Get the orbital parameters."""
-        ssp_lon, ssp_lat = self._get_ssp_lonlat()
         orbital_parameters = {
             'projection_longitude': self.projection_longitude,
             'projection_latitude': 0.0,
             'projection_altitude': ALTITUDE
         }
+        ssp_lon, ssp_lat = self._get_ssp_lonlat()
         if not np.isnan(ssp_lon) and not np.isnan(ssp_lat):
             orbital_parameters.update({
                 'satellite_actual_longitude': ssp_lon,
                 'satellite_actual_latitude': ssp_lat,
                 # altitude not available
             })
+        return orbital_parameters
 
     def _get_ssp_lonlat(self):
         """Get longitude and latitude at the subsatellite point.

From 34776410985bab44d5d3603ec7304f00cb4ba5ec Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 18 Nov 2020 14:49:26 +0000
Subject: [PATCH 11/47] Fix tiepoint interpolation

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 2a145fc32e..1bd9a47436 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -467,7 +467,8 @@ def _interp_tiepoints(self, ds, target_x, target_y):
         ds = ds.assign_coords(x_tie=target_x.values[::sampling],
                               y_tie=target_y.values[::sampling])
 
-        ds_interp = ds.interp(x_tie=target_x, y_tie=target_y)
+        ds_interp = ds.interp(x_tie=target_x.values,
+                              y_tie=target_y.values)
         return ds_interp.rename({'x_tie': 'x', 'y_tie': 'y'})
 
 

From 585f9711bd303933d4356bc31b3689e1a201c41f Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 18 Nov 2020 14:50:08 +0000
Subject: [PATCH 12/47] Define IR/WV netCDF keys in the baseclass

This facilitates testing.
---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 38 ++++++++++++++++++++-------
 1 file changed, 28 insertions(+), 10 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 1bd9a47436..e2ac043581 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -173,6 +173,30 @@ def __get__(self, obj, objtype):
 class FiduceoMviriBase(BaseFileHandler):
     """Baseclass for FIDUCEO MVIRI file handlers."""
 
+    nc_keys = {'WV': 'count_wv', 'IR': 'count_ir'}
+    nc_keys_coefs = {
+        'WV': {
+            'radiance': {
+                'a': 'a_wv',
+                'b': 'b_wv'
+            },
+            'brightness_temperature': {
+                'a': 'bt_a_wv',
+                'b': 'bt_b_wv'
+            }
+        },
+        'IR': {
+            'radiance': {
+                'a': 'a_ir',
+                'b': 'b_ir'
+            },
+            'brightness_temperature': {
+                'a': 'bt_a_ir',
+                'b': 'bt_b_ir'
+            }
+        },
+    }
+
     def __init__(self, filename, filename_info, filetype_info,
                  mask_bad_quality=False):
         """Initialize the file handler.
@@ -475,11 +499,8 @@ def _interp_tiepoints(self, ds, target_x, target_y):
 class FiduceoMviriEasyFcdrFileHandler(FiduceoMviriBase):
     """File handler for FIDUCEO MVIRI Easy FCDR."""
 
-    nc_keys = {
-        'VIS': 'toa_bidirectional_reflectance_vis',
-        'WV': 'count_wv',
-        'IR': 'count_ir'
-    }
+    nc_keys = FiduceoMviriBase.nc_keys.copy()
+    nc_keys['VIS'] = 'toa_bidirectional_reflectance_vis'
 
     def _calibrate_vis(self, ds, calibration):
         """Calibrate VIS channel.
@@ -500,11 +521,8 @@ def _calibrate_vis(self, ds, calibration):
 class FiduceoMviriFullFcdrFileHandler(FiduceoMviriBase):
     """File handler for FIDUCEO MVIRI Full FCDR."""
 
-    nc_keys = {
-        'VIS': 'count_vis',
-        'WV': 'count_wv',
-        'IR': 'count_ir'
-    }
+    nc_keys = FiduceoMviriBase.nc_keys.copy()
+    nc_keys['VIS'] = 'count_vis'
 
     def _calibrate_vis(self, ds, calibration):
         """Calibrate VIS channel.

From 6a1a7685b1b929f43f40468f110792ec3e40467e Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 18 Nov 2020 15:40:44 +0000
Subject: [PATCH 13/47] Mask VIS channel only

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 29 ++++++++-------------------
 1 file changed, 8 insertions(+), 21 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index e2ac043581..ad10ff2e28 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -202,7 +202,7 @@ def __init__(self, filename, filename_info, filetype_info,
         """Initialize the file handler.
 
         Args:
-             mask_bad_quality: Mask pixels with bad quality, that means
+             mask_bad_quality: Mask VIS pixels with bad quality, that means
                  everything else than "ok" or "use with caution". If you need
                  more control, use the ``quality_pixel_bitmask`` and
                  ``data_quality_bitmask`` datasets.
@@ -229,8 +229,8 @@ def get_dataset(self, dataset_id, info):
         if dataset_id['name'] in CHANNELS:
             ds = self.calibrate(ds, channel=name,
                                 calibration=dataset_id['calibration'])
-            if self.mask_bad_quality:
-                ds = self._mask(ds)
+            if self.mask_bad_quality and name == 'VIS':
+                ds = self._mask_vis(ds)
             ds['acq_time'] = ('y', self._get_acq_time(ds))
         elif dataset_id['name'] in OTHER_REFLECTANCES:
             ds = ds * 100  # conversion to percent
@@ -349,29 +349,16 @@ def _ir_wv_radiance_to_brightness_temperature(self, rad, channel):
         bt = b / (np.log(rad) - a)
         return bt.where(bt > 0)
 
-    def _mask(self, ds):
-        """Mask pixels with bad quality.
+    def _mask_vis(self, ds):
+        """Mask VIS pixels with bad quality.
 
         Pixels are considered bad quality if the "quality_pixel_bitmask" is
         everything else than 0 (no flag set) or 2 ("use_with_caution" and no
-        other flag set).
+        other flag set). According to [PUG] that bitmask is only applicable to
+        the VIS channel.
         """
-        mask = self._get_mask(ds)
-        return ds.where(np.logical_or(mask == 0, mask == 2))
-
-    @shape_cache
-    def _get_mask(self, ds):
-        """Get quality bitmask for the given dataset."""
         mask = self.nc['quality_pixel_bitmask']
-
-        # Mask has high (VIS) resolution. Downsample to low (IR/WV) resolution
-        # if required.
-        if mask.coords['y'].size > ds.coords['y'].size:
-            mask = mask.isel(y=slice(None, None, 2),
-                             x=slice(None, None, 2))
-            mask = mask.assign_coords(y=ds.coords['y'].values,
-                                      x=ds.coords['x'].values)
-        return mask
+        return ds.where(np.logical_or(mask == 0, mask == 2))
 
     @shape_cache
     def _get_acq_time(self, ds):

From 6983fafd0c7e29cd3ca83de0da9da262fc011376 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 19 Nov 2020 10:15:05 +0000
Subject: [PATCH 14/47] Add more angle datasets

---
 satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml |  56 +++++++
 satpy/readers/mviri_l1b_fiduceo_nc.py       | 172 +++++++++++++-------
 2 files changed, 170 insertions(+), 58 deletions(-)

diff --git a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
index 5ddacf574c..c62ec6ab6c 100644
--- a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
+++ b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
@@ -101,3 +101,59 @@ datasets:
     units: "%"
     resolution: 2250
     file_type: [nc_easy]
+
+  solar_zenith_angle_vis:
+    name: solar_zenith_angle_vis
+    long_name: "Solar zenith angle VIS channel"
+    units: "degrees"
+    resolution: 2250
+    file_type: [nc_easy, nc_full]
+
+  solar_azimuth_angle_vis:
+    name: solar_azimuth_angle_vis
+    long_name: "Solar azimuth angle VIS channel"
+    units: "degrees"
+    resolution: 2250
+    file_type: [nc_easy, nc_full]
+
+  satellite_zenith_angle_vis:
+    name: satellite_zenith_angle_vis
+    long_name: "Satellite zenith angle VIS channel"
+    units: "degrees"
+    resolution: 2250
+    file_type: [nc_easy, nc_full]
+
+  satellite_azimuth_angle_vis:
+    name: satellite_azimuth_angle_vis
+    long_name: "Satellite azimuth angle VIS channel"
+    units: "degrees"
+    resolution: 2250
+    file_type: [nc_easy, nc_full]
+
+  solar_zenith_angle_ir_wv:
+    name: solar_zenith_angle_ir_wv
+    long_name: "Solar zenith angle IR/WV channels"
+    units: "degrees"
+    resolution: 2250
+    file_type: [nc_easy, nc_full]
+
+  solar_azimuth_angle_ir_wv:
+    name: solar_azimuth_angle_ir_wv
+    long_name: "Solar azimuth angle IR/WV channels"
+    units: "degrees"
+    resolution: 2250
+    file_type: [nc_easy, nc_full]
+
+  satellite_zenith_angle_ir_wv:
+    name: satellite_zenith_angle_ir_wv
+    long_name: "Satellite zenith angle IR/WV channels"
+    units: "degrees"
+    resolution: 2250
+    file_type: [nc_easy, nc_full]
+
+  satellite_azimuth_angle_ir_wv:
+    name: satellite_azimuth_angle_ir_wv
+    long_name: "Satellite azimuth angle IR/WV channels"
+    units: "degrees"
+    resolution: 2250
+    file_type: [nc_easy, nc_full]
\ No newline at end of file
diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index ad10ff2e28..983324d72b 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -119,6 +119,7 @@
 from satpy.readers.file_handlers import BaseFileHandler
 from satpy.readers._geos_area import (ang2fac, get_area_definition,
                                       get_area_extent)
+from satpy.tests.utils import make_dataid
 
 
 EQUATOR_RADIUS = 6378140.0
@@ -130,50 +131,79 @@
 """[Handbook] section 5.3.2.1."""
 
 CHANNELS = ['VIS', 'WV', 'IR']
+ANGLES = [
+    'solar_zenith_angle_vis',
+    'solar_azimuth_angle_vis',
+    'satellite_zenith_angle_vis',
+    'satellite_azimuth_angle_vis',
+    'solar_zenith_angle_ir_wv',
+    'solar_azimuth_angle_ir_wv',
+    'satellite_zenith_angle_ir_wv',
+    'satellite_azimuth_angle_ir_wv'
+]
 OTHER_REFLECTANCES = [
     'u_independent_toa_bidirectional_reflectance',
     'u_structured_toa_bidirectional_reflectance'
 ]
 
 
-class shape_cache:
-    """Cache function call based on image shape.
-
-    Shape can be used to maintain separate caches for low resolution
-    (WV/IR) and high resolution (VIS) channels.
-    """
+class InterpCache:
+    """Interpolation cache."""
 
-    def __init__(self, func):
+    def __init__(self, func, keys, hash_funcs):
         """Create the cache.
 
         Args:
             func:
-                Function to me cached.
+                Interpolation function to be cached.
+            keys:
+                Function arguments serving as cache key.
+            hash_funcs (dict):
+                For each key provides a method that extracts a hashable
+                value from the given keyword argument. For example, use
+                image shape to maintain separate caches for low
+                resolution (WV/IR) and high resolution (VIS) channels.
         """
         self.func = func
+        self.keys = keys
+        self.hash_funcs = hash_funcs
         self.cache = {}
 
-    def __call__(self, *args):
-        """Call the decorated function.
-
-        Dataset is expected to be the last argument. If an instance method is
-        decorated, it is preceded by a filehandler instance.
-        """
-        ds = args[-1]
-        shape = ds.coords['y'].shape
-        if shape not in self.cache:
-            self.cache[shape] = self.func(*args)
-        return self.cache[shape]
+    def __call__(self, *args, **kwargs):
+        """Call the interpolation function."""
+        key = tuple(
+            [self.hash_funcs[key](kwargs[key]) for key in self.keys]
+        )
+        if key not in self.cache:
+            self.cache[key] = self.func(*args, **kwargs)
+        return self.cache[key]
 
     def __get__(self, obj, objtype):
         """To support instance methods."""
         return functools.partial(self.__call__, obj)
 
 
+def interp_cache(keys, hash_funcs):
+    """Interpolation cache."""
+    def wrapper(func):
+        return InterpCache(func, keys, hash_funcs)
+    return wrapper
+
+
 class FiduceoMviriBase(BaseFileHandler):
     """Baseclass for FIDUCEO MVIRI file handlers."""
-
-    nc_keys = {'WV': 'count_wv', 'IR': 'count_ir'}
+    nc_keys = {
+        'WV': 'count_wv',
+        'IR': 'count_ir',
+        'solar_zenith_angle_vis': 'solar_zenith_angle',
+        'solar_azimuth_angle_vis': 'solar_azimuth_angle',
+        'satellite_zenith_angle_vis': 'satellite_zenith_angle',
+        'satellite_azimuth_angle_vis': 'satellite_azimuth_angle',
+        'solar_zenith_angle_ir_wv': 'solar_zenith_angle',
+        'solar_azimuth_angle_ir_wv': 'solar_azimuth_angle',
+        'satellite_zenith_angle_ir_wv': 'satellite_zenith_angle',
+        'satellite_azimuth_angle_ir_wv': 'satellite_azimuth_angle',
+    }
     nc_keys_coefs = {
         'WV': {
             'radiance': {
@@ -234,6 +264,8 @@ def get_dataset(self, dataset_id, info):
             ds['acq_time'] = ('y', self._get_acq_time(ds))
         elif dataset_id['name'] in OTHER_REFLECTANCES:
             ds = ds * 100  # conversion to percent
+        elif dataset_id['name'] in ANGLES:
+            ds = self._interp_angles(ds, dataset_id['name'])
         self._update_attrs(ds, info)
         return ds
 
@@ -247,6 +279,8 @@ def _read_dataset(self, name):
         ds = self.nc[nc_key]
         if 'y_ir_wv' in ds.dims:
             ds = ds.rename({'y_ir_wv': 'y', 'x_ir_wv': 'x'})
+        if 'y_tie' in ds.dims:
+            ds = ds.rename({'y_tie': 'y', 'x_tie': 'x'})
         return ds
 
     def _update_attrs(self, ds, info):
@@ -360,9 +394,21 @@ def _mask_vis(self, ds):
         mask = self.nc['quality_pixel_bitmask']
         return ds.where(np.logical_or(mask == 0, mask == 2))
 
-    @shape_cache
     def _get_acq_time(self, ds):
-        """Get scanline acquisition time.
+        """Get scanline acquisition time for the given dataset."""
+        # Variable is sometimes named "time" and sometimes "time_ir_wv".
+        try:
+            time2d = self.nc['time_ir_wv']
+        except KeyError:
+            time2d = self.nc['time']
+        return self._get_acq_time_cached(time2d, target_y=ds.coords['y'])
+
+    @interp_cache(
+        keys=('target_y',),
+        hash_funcs={'target_y': lambda y: y.size}
+    )
+    def _get_acq_time_cached(self, time2d, target_y):
+        """Get scanline acquisition time for the given image coordinates.
 
         The files provide timestamps per pixel for the low resolution
         channels (IR/WV) only.
@@ -377,27 +423,19 @@ def _get_acq_time(self, ds):
         Returns:
             Mean scanline acquisition timestamps
         """
-        # Variable is sometimes named "time" and sometimes "time_ir_wv".
-        try:
-            time_lores = self.nc['time_ir_wv']
-        except KeyError:
-            time_lores = self.nc['time']
-
         # Compute mean timestamp per scanline
-        time_lores = time_lores.mean(dim='x_ir_wv').rename({'y_ir_wv': 'y'})
+        time = time2d.mean(dim='x_ir_wv').rename({'y_ir_wv': 'y'})
 
         # If required, repeat timestamps in y-direction to obtain higher
         # resolution
-        y_lores = time_lores.coords['y'].values
-        y_target = ds.coords['y'].values
-        if y_lores.size < y_target.size:
-            reps = y_target.size // y_lores.size
-            y_lores_rep = np.repeat(y_lores, reps)
-            time_hires = time_lores.reindex(y=y_lores_rep)
-            time_hires = time_hires.assign_coords(y=y_target)
+        y = time.coords['y'].values
+        if y.size < target_y.size:
+            reps = target_y.size // y.size
+            y_rep = np.repeat(y, reps)
+            time_hires = time.reindex(y=y_rep)
+            time_hires = time_hires.assign_coords(y=target_y)
             return time_hires
-
-        return time_lores
+        return time
 
     def _get_orbital_parameters(self):
         """Get the orbital parameters."""
@@ -442,16 +480,38 @@ def _get_ssp(self, coord):
             sub_lonlat = np.nan
         return sub_lonlat
 
-    @shape_cache
-    def _get_solar_zenith_angle(self, ds):
-        """Get solar zenith angle for the given dataset.
+    def _interp_angles(self, angles, name):
+        """Get angle dataset.
 
-        Files provide solar zenith angle, but at a coarser resolution.
-        Interpolate to the resolution of the given dataset.
+        Files provide angles (solar/satellite zenith & azimuth) at a coarser
+        resolution. Interpolate them to the desired resolution.
         """
-        return self._interp_tiepoints(self.nc['solar_zenith_angle'],
-                                      ds.coords['x'],
-                                      ds.coords['y'])
+        if name.endswith('_vis'):
+            target_x = self.nc.coords['x']
+            target_y = self.nc.coords['y']
+        else:
+            target_x = self.nc.coords['x_ir_wv']
+            target_y = self.nc.coords['y_ir_wv']
+        return self._interp_angles_cached(
+            angles=angles,
+            nc_key=self.nc_keys[name],
+            target_x=target_x,
+            target_y=target_y
+        )
+
+    @interp_cache(
+        keys=('nc_key', 'target_x', 'target_y'),
+        hash_funcs={
+            'nc_key': lambda nc_key: nc_key,
+            'target_x': lambda x: x.size,
+            'target_y': lambda y: y.size
+        }
+    )
+    def _interp_angles_cached(self, angles, nc_key, target_x, target_y):
+        """Interpolate angles to the given resolution."""
+        return self._interp_tiepoints(angles,
+                                      target_x,
+                                      target_y)
 
     def _interp_tiepoints(self, ds, target_x, target_y):
         """Interpolate dataset between tiepoints.
@@ -468,19 +528,14 @@ def _interp_tiepoints(self, ds, target_x, target_y):
             target_y:
                 Target y coordinates
         """
-        if 'x_tie' not in ds.dims:
-            raise ValueError('Dataset has no tiepoints to be interpolated.')
-
         # No tiepoint coordinates specified in the files. Use dimensions
         # to calculate tiepoint sampling and assign tiepoint coordinates
         # accordingly.
-        sampling = target_x.size // ds.coords['x_tie'].size
-        ds = ds.assign_coords(x_tie=target_x.values[::sampling],
-                              y_tie=target_y.values[::sampling])
+        sampling = target_x.size // ds.coords['x'].size
+        ds = ds.assign_coords(x=target_x.values[::sampling],
+                              y=target_y.values[::sampling])
 
-        ds_interp = ds.interp(x_tie=target_x.values,
-                              y_tie=target_y.values)
-        return ds_interp.rename({'x_tie': 'x', 'y_tie': 'y'})
+        return ds.interp(x=target_x.values, y=target_y.values)
 
 
 class FiduceoMviriEasyFcdrFileHandler(FiduceoMviriBase):
@@ -550,8 +605,9 @@ def _vis_radiance_to_reflectance(self, rad):
 
         Reference: [PUG], equation (6).
         """
-
-        sza = self._get_solar_zenith_angle(rad)
+        sza = self.get_dataset(
+            make_dataid(name='solar_zenith_angle_vis'), info={}
+        )
         sza = sza.where(da.fabs(sza) < 90)  # direct illumination only
         cos_sza = np.cos(np.deg2rad(sza))
         refl = (

From 6168bb4d33a16a99090801984acead6a5d0812a7 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 19 Nov 2020 10:18:43 +0000
Subject: [PATCH 15/47] Add unit tests

---
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py | 392 ++++++++++++++++++
 1 file changed, 392 insertions(+)
 create mode 100644 satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py

diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
new file mode 100644
index 0000000000..928face502
--- /dev/null
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -0,0 +1,392 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2020 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 <http://www.gnu.org/licenses/>.
+"""Unit tests for the FIDUCEO MVIRI FCDR Reader."""
+
+import dask.array as da
+import numpy as np
+import pytest
+from unittest import mock
+import xarray as xr
+
+from satpy.tests.utils import make_dataid
+from satpy.readers.mviri_l1b_fiduceo_nc import (
+    FiduceoMviriEasyFcdrFileHandler,
+    FiduceoMviriFullFcdrFileHandler
+)
+
+
+attrs_exp = {
+    'platform': 'MET7',
+    'raw_metadata': {'foo': 'bar'},
+    'sensor': 'MVIRI',
+    'orbital_parameters': {
+        'projection_longitude': 57.0,
+        'projection_latitude': 0.0,
+        'projection_altitude': 35785860.0,
+        'satellite_actual_longitude': 57.1,
+        'satellite_actual_latitude': 0.1,
+    }
+
+}
+attrs_refl_exp = attrs_exp.copy()
+attrs_refl_exp.update(
+    {'sun_earth_distance_correction_applied': True,
+     'sun_earth_distance_correction_factor': 1.}
+)
+acq_time_hires_exp = [np.datetime64('1970-01-01 00:30'),
+                      np.datetime64('1970-01-01 00:30'),
+                      np.datetime64('1970-01-01 02:30'),
+                      np.datetime64('1970-01-01 02:30')]
+vis_counts_exp = xr.DataArray(
+    [[0., 17., 34., 51.],
+     [68., 85., 102., 119.],
+     [136., 153., np.nan, 187.],
+     [204., 221., 238., 255]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2, 3, 4],
+        'x': [1, 2, 3, 4],
+        'acq_time': ('y', acq_time_hires_exp),
+    },
+    attrs=attrs_exp
+)
+vis_rad_exp = xr.DataArray(
+    [[np.nan, 18.56, 38.28, 58.],
+     [77.72, 97.44, 117.16, 136.88],
+     [156.6, 176.32, np.nan, 215.76],
+     [235.48, 255.2, 274.92, 294.64]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2, 3, 4],
+        'x': [1, 2, 3, 4],
+        'acq_time': ('y', acq_time_hires_exp),
+    },
+    attrs=attrs_exp
+)
+vis_refl_exp = xr.DataArray(
+    [[np.nan, 23.440929, np.nan, np.nan],
+     [40.658744, 66.602233, 147.970867, np.nan],
+     [75.688217, 92.240733, np.nan, np.nan],
+     [np.nan, np.nan, np.nan, np.nan]],
+    # (0, 0) and (2, 2) are NaN because radiance is NaN
+    # (0, 2) is NaN because SZA >= 90 degrees
+    # Last row/col is NaN due to SZA interpolation
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2, 3, 4],
+        'x': [1, 2, 3, 4],
+        'acq_time': ('y', acq_time_hires_exp),
+    },
+    attrs=attrs_refl_exp
+)
+acq_time_lores_exp = [np.datetime64('1970-01-01 00:30'),
+                      np.datetime64('1970-01-01 02:30')]
+wv_counts_exp = xr.DataArray(
+    [[0, 85],
+     [170, 255]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2],
+        'x': [1, 2],
+        'acq_time': ('y', acq_time_lores_exp),
+    },
+    attrs=attrs_exp
+)
+wv_rad_exp = xr.DataArray(
+    [[np.nan, 3.75],
+     [8, 12.25]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2],
+        'x': [1, 2],
+        'acq_time': ('y', acq_time_lores_exp),
+    },
+    attrs=attrs_exp
+)
+wv_bt_exp = xr.DataArray(
+    [[np.nan, 230.461366],
+     [252.507448, 266.863289]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2],
+        'x': [1, 2],
+        'acq_time': ('y', acq_time_lores_exp),
+    },
+    attrs=attrs_exp
+)
+ir_counts_exp = xr.DataArray(
+    [[0, 85],
+     [170, 255]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2],
+        'x': [1, 2],
+        'acq_time': ('y', acq_time_lores_exp),
+    },
+    attrs=attrs_exp
+)
+ir_rad_exp = xr.DataArray(
+    [[np.nan, 80],
+     [165, 250]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2],
+        'x': [1, 2],
+        'acq_time': ('y', acq_time_lores_exp),
+    },
+    attrs=attrs_exp
+)
+ir_bt_exp = xr.DataArray(
+    [[np.nan, 178.00013189],
+     [204.32955838, 223.28709913]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2],
+        'x': [1, 2],
+        'acq_time': ('y', acq_time_lores_exp),
+    },
+    attrs=attrs_exp
+)
+quality_pixel_bitmask_exp = xr.DataArray(
+    [[2, 0, 0, 0],
+     [0, 0, 0, 0],
+     [0, 0, 1, 0],
+     [0, 0, 0, 0]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2, 3, 4],
+        'x': [1, 2, 3, 4],
+    },
+    attrs=attrs_exp
+)
+sza_vis_exp = xr.DataArray(
+    [[45., 67.5, 90., np.nan],
+     [22.5, 45., 67.5, np.nan],
+     [0., 22.5, 45., np.nan],
+     [np.nan, np.nan, np.nan, np.nan]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2, 3, 4],
+        'x': [1, 2, 3, 4],
+    },
+    attrs=attrs_exp
+)
+sza_ir_wv_exp = xr.DataArray(
+    [[45, 90],
+     [0, 45]],
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2],
+        'x': [1, 2],
+    },
+    attrs=attrs_exp
+)
+
+
+@pytest.fixture()
+def fake_dataset():
+    """Create fake dataset."""
+    count_ir = da.linspace(0, 255, 4, dtype=np.float32).reshape(2, 2)
+    count_wv = da.linspace(0, 255, 4, dtype=np.float32).reshape(2, 2)
+    count_vis = da.linspace(0, 255, 16, dtype=np.float32).reshape(4, 4)
+    sza = da.from_array(
+        [[45, 90],
+         [0, 45]]
+    ).astype(np.float32)
+    mask = da.from_array(
+        [[2, 0, 0, 0],  # 2 = "use with caution"
+         [0, 0, 0, 0],
+         [0, 0, 1, 0],  # 1 = "invalid"
+         [0, 0, 0, 0]]
+    )
+    time = np.arange(4).astype('datetime64[h]').reshape(2, 2)
+    ds = xr.Dataset(
+        data_vars={
+            'count_vis': (('y', 'x'), count_vis),
+            'count_wv': (('y_ir_wv', 'x_ir_wv'), count_wv),
+            'count_ir': (('y_ir_wv', 'x_ir_wv'), count_ir),
+            'toa_bidirectional_reflectance_vis': (
+                ('y', 'x'), vis_refl_exp / 100),
+            'quality_pixel_bitmask': (('y', 'x'), mask),
+            'solar_zenith_angle': (('y_tie', 'x_tie'), sza),
+            'time_ir_wv': (('y_ir_wv', 'x_ir_wv'), time),
+            'a_ir': -5.0,
+            'b_ir': 1.0,
+            'bt_a_ir': 10.0,
+            'bt_b_ir': -1000.0,
+            'a_wv': -0.5,
+            'b_wv': 0.05,
+            'bt_a_wv': 10.0,
+            'bt_b_wv': -2000.0,
+            'years_since_launch': 20.0,
+            'a0_vis': 1.0,
+            'a1_vis': 0.01,
+            'a2_vis': -0.0001,
+            'mean_count_space_vis': 1.0,
+            'distance_sun_earth': 1.0,
+            'solar_irradiance_vis': 650.0,
+            'sub_satellite_longitude_start': 57.1,
+            'sub_satellite_longitude_end': np.nan,
+            'sub_satellite_latitude_start': np.nan,
+            'sub_satellite_latitude_end': 0.1,
+        },
+        coords={
+            'y': [1, 2, 3, 4],
+            'x': [1, 2, 3, 4],
+            'y_ir_wv': [1, 2],
+            'x_ir_wv': [1, 2],
+            'y_tie': [1, 2],
+            'x_tie': [1, 2]
+
+        },
+        attrs={'foo': 'bar'}
+    )
+    return ds
+
+
+@pytest.fixture(params=[FiduceoMviriEasyFcdrFileHandler,
+                        FiduceoMviriFullFcdrFileHandler])
+def file_handler(fake_dataset, request):
+    """Create mocked file handler."""
+    fh_class = request.param
+    with mock.patch('satpy.readers.mviri_l1b_fiduceo_nc.xr.open_dataset') as open_dataset:
+        open_dataset.return_value = fake_dataset
+        return fh_class(
+            filename='filename',
+            filename_info={'platform': 'MET7',
+                           'sensor': 'MVIRI',
+                           'projection_longitude': '57.0'},
+            filetype_info={'foo': 'bar'},
+            mask_bad_quality=True
+        )
+
+
+
+
+
+# TODO: Test file patterns
+
+
+
+
+
+
+class TestFiduceoMviriFileHandlers:
+    """Unit tests for FIDUCEO MVIRI file handlers."""
+
+    def test_init(self, file_handler):
+        """Test file handler initialization."""
+        assert file_handler.projection_longitude == 57.0
+        assert file_handler.mask_bad_quality is True
+
+    @pytest.mark.parametrize(
+        ('name', 'calibration', 'expected'),
+        [
+            ('VIS', 'counts', vis_counts_exp),
+            ('VIS', 'radiance', vis_rad_exp),
+            ('VIS', 'reflectance', vis_refl_exp),
+            ('WV', 'counts', wv_counts_exp),
+            ('WV', 'radiance', wv_rad_exp),
+            ('WV', 'brightness_temperature', wv_bt_exp),
+            ('IR', 'counts', ir_counts_exp),
+            ('IR', 'radiance', ir_rad_exp),
+            ('IR', 'brightness_temperature', ir_bt_exp),
+            ('quality_pixel_bitmask', None, quality_pixel_bitmask_exp),
+            ('solar_zenith_angle_vis', None, sza_vis_exp),
+            ('solar_zenith_angle_ir_wv', None, sza_ir_wv_exp)
+        ]
+    )
+    def test_get_dataset(self, file_handler, name, calibration, expected):
+        """Test getting datasets."""
+        id_keys = {'name': name}
+        if calibration:
+            id_keys['calibration'] = calibration
+        data_id = make_dataid(**id_keys)
+        info = {'platform': 'MET7'}
+
+        is_easy = isinstance(file_handler, FiduceoMviriEasyFcdrFileHandler)
+        is_vis = name == 'VIS'
+        is_refl = calibration == 'reflectance'
+        if is_easy and is_vis and not is_refl:
+            # VIS counts/radiance not available in easy FCDR
+            with pytest.raises(ValueError):
+                file_handler.get_dataset(data_id, info)
+        else:
+            ds = file_handler.get_dataset(data_id, info)
+            xr.testing.assert_allclose(ds, expected)
+            assert ds.attrs == expected.attrs
+
+    def test_time_cache(self, file_handler):
+        """Test caching of acquisition times."""
+        time2d = xr.DataArray(
+            np.array([[1, 2],
+                      [3, 4]],
+                     dtype='datetime64[h]'),
+            dims=('y_ir_wv', 'x_ir_wv')
+        )
+        y1 = xr.DataArray([1, 2, 3, 4])
+        t1 = file_handler._get_acq_time_cached(time2d, target_y=y1)
+
+        # Change 2d timestamps. If the cache works correctly, the second call
+        # should not average/interpolate them again.
+        t2 = file_handler._get_acq_time_cached(
+            time2d + np.timedelta64(1, 'h'), target_y=y1)
+        xr.testing.assert_equal(t2, t1)
+
+        # With new target coordinates we shouldn't hit the cache
+        y2 = xr.DataArray([1, 2])
+        t3 = file_handler._get_acq_time_cached(target_y=y2)
+        with pytest.raises(AssertionError):
+            xr.testing.assert_equal(t3, t1)
+
+    def test_angle_cache(self, file_handler):
+        """Test caching of angle datasets."""
+        nc_key = 'mykey'
+        x1 = y1 = xr.DataArray([1, 2, 3, 4])
+        sza_coarse = xr.DataArray(
+            [[45, 90],
+             [0, 45]],
+            dims=('y', 'x')
+        )
+        sza1 = file_handler._interp_angles_cached(
+            angles=sza_coarse,
+            nc_key=nc_key,
+            target_x=x1,
+            target_y=y1
+        )
+
+        # Change coarse angles. If the cache works correctly, the second call
+        # should not interpolate them again.
+        sza2 = file_handler._interp_angles_cached(
+            angles=sza_coarse - 10,
+            nc_key=nc_key,
+            target_x=x1,
+            target_y=y1
+        )
+        xr.testing.assert_equal(sza2, sza1)
+
+        # With new target coordinates we shouldn't hit the cache
+        x2 = y2 = xr.DataArray([1, 2])
+        sza3 = file_handler._interp_angles_cached(
+            angles=sza_coarse,
+            nc_key=nc_key,
+            target_x=x2,
+            target_y=y2
+        )
+        with pytest.raises(AssertionError):
+            xr.testing.assert_equal(sza3, sza1)

From 499b52852102b246b69b23f845391ab1d88496d5 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 19 Nov 2020 13:30:12 +0000
Subject: [PATCH 16/47] Fix standard names

---
 satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml | 32 +++++++++++++--------
 1 file changed, 20 insertions(+), 12 deletions(-)

diff --git a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
index c62ec6ab6c..2ed3d7063f 100644
--- a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
+++ b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
@@ -39,8 +39,8 @@ datasets:
         standard_name: toa_bidirectional_reflectance
         units: "%"
       radiance:
-        # standard_name: TODO
-        units: W m-2 sr-1  # TODO: Per unit wave number/length
+        # Confirmed by EUM: No (1/wavenumber) here. Hence no standard name.
+        units: W m-2 sr-1
       counts:
         standard_name: counts
         units: count
@@ -55,7 +55,7 @@ datasets:
         standard_name: toa_brightness_temperature
         units: K
       radiance:
-        standard_name: toa_outgoing_radiance_per_unit_wavelength
+        standard_name: toa_outgoing_radiance_per_unit_wavenumber
         units: mW m-2 sr-1 (cm-1)-1
       counts:
         standard_name: counts
@@ -71,7 +71,7 @@ datasets:
         standard_name: toa_brightness_temperature
         units: K
       radiance:
-        standard_name: toa_outgoing_radiance_per_unit_wavelength
+        standard_name: toa_outgoing_radiance_per_unit_wavenumber
         units: mW m-2 sr-1 (cm-1)-1
       counts:
         standard_name: counts
@@ -104,56 +104,64 @@ datasets:
 
   solar_zenith_angle_vis:
     name: solar_zenith_angle_vis
+    standard_name: solar_zenith_angle
     long_name: "Solar zenith angle VIS channel"
-    units: "degrees"
+    units: degree
     resolution: 2250
     file_type: [nc_easy, nc_full]
 
   solar_azimuth_angle_vis:
     name: solar_azimuth_angle_vis
+    standard_name: solar_azimuth_angle
     long_name: "Solar azimuth angle VIS channel"
-    units: "degrees"
+    units: degree
     resolution: 2250
     file_type: [nc_easy, nc_full]
 
   satellite_zenith_angle_vis:
     name: satellite_zenith_angle_vis
+    standard_name: sensor_zenith_angle
     long_name: "Satellite zenith angle VIS channel"
-    units: "degrees"
+    units: degree
     resolution: 2250
     file_type: [nc_easy, nc_full]
 
   satellite_azimuth_angle_vis:
     name: satellite_azimuth_angle_vis
+    standard_name: sensor_azimuth_angle
     long_name: "Satellite azimuth angle VIS channel"
-    units: "degrees"
+    units: degree
     resolution: 2250
     file_type: [nc_easy, nc_full]
 
   solar_zenith_angle_ir_wv:
     name: solar_zenith_angle_ir_wv
+    standard_name: solar_zenith_angle
     long_name: "Solar zenith angle IR/WV channels"
-    units: "degrees"
+    units: degree
     resolution: 2250
     file_type: [nc_easy, nc_full]
 
   solar_azimuth_angle_ir_wv:
     name: solar_azimuth_angle_ir_wv
+    standard_name: solar_azimuth_angle
     long_name: "Solar azimuth angle IR/WV channels"
-    units: "degrees"
+    units: degree
     resolution: 2250
     file_type: [nc_easy, nc_full]
 
   satellite_zenith_angle_ir_wv:
     name: satellite_zenith_angle_ir_wv
+    standard_name: sensor_zenith_angle
     long_name: "Satellite zenith angle IR/WV channels"
-    units: "degrees"
+    units: degree
     resolution: 2250
     file_type: [nc_easy, nc_full]
 
   satellite_azimuth_angle_ir_wv:
     name: satellite_azimuth_angle_ir_wv
+    standard_name: sensor_azimuth_angle
     long_name: "Satellite azimuth angle IR/WV channels"
-    units: "degrees"
+    units: degree
     resolution: 2250
     file_type: [nc_easy, nc_full]
\ No newline at end of file

From 599d54297902462abe0847d5b650d1edf451129a Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 19 Nov 2020 15:43:12 +0000
Subject: [PATCH 17/47] Add test for area defs

---
 satpy/readers/mviri_l1b_fiduceo_nc.py         |  35 ++++--
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py | 119 +++++++++++++-----
 2 files changed, 109 insertions(+), 45 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 983324d72b..f1923441ae 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -119,7 +119,7 @@
 from satpy.readers.file_handlers import BaseFileHandler
 from satpy.readers._geos_area import (ang2fac, get_area_definition,
                                       get_area_extent)
-from satpy.tests.utils import make_dataid
+from satpy.dataset.dataid import DataQuery
 
 
 EQUATOR_RADIUS = 6378140.0
@@ -145,6 +145,7 @@
     'u_independent_toa_bidirectional_reflectance',
     'u_structured_toa_bidirectional_reflectance'
 ]
+HIGH_RESOL = 2250
 
 
 class InterpCache:
@@ -252,7 +253,7 @@ def __init__(self, filename, filename_info, filetype_info,
         # filename.
         self.projection_longitude = float(filename_info['projection_longitude'])
 
-    def get_dataset(self, dataset_id, info):
+    def get_dataset(self, dataset_id, dataset_info):
         """Get the dataset."""
         name = dataset_id['name']
         ds = self._read_dataset(name)
@@ -265,8 +266,8 @@ def get_dataset(self, dataset_id, info):
         elif dataset_id['name'] in OTHER_REFLECTANCES:
             ds = ds * 100  # conversion to percent
         elif dataset_id['name'] in ANGLES:
-            ds = self._interp_angles(ds, dataset_id['name'])
-        self._update_attrs(ds, info)
+            ds = self._interp_angles(ds, dataset_id)
+        self._update_attrs(ds, dataset_info)
         return ds
 
     def _get_nc_key(self, name):
@@ -303,6 +304,9 @@ def get_area_def(self, dataset_id):
         loff = coff = im_size / 2 + 0.5
         lfac = cfac = ang2fac(np.deg2rad(MVIRI_FIELD_OF_VIEW) / im_size)
 
+        area_name = 'geos_mviri_{}'.format(
+            'vis' if self._is_high_resol(dataset_id) else 'ir_wv'
+        )
         pdict = {
             'ssp_lon': self.projection_longitude,
             'a': EQUATOR_RADIUS,
@@ -316,8 +320,8 @@ def get_area_def(self, dataset_id):
             'nlines': im_size,
             'ncols': im_size,
             'scandir': 'S2N',  # Reference: [PUG] section 2.
-            'p_id': 'geos_mviri',
-            'a_name': 'geos_mviri',
+            'p_id': area_name,
+            'a_name': area_name,
             'a_desc': 'MVIRI Geostationary Projection'
         }
         extent = get_area_extent(pdict)
@@ -395,7 +399,11 @@ def _mask_vis(self, ds):
         return ds.where(np.logical_or(mask == 0, mask == 2))
 
     def _get_acq_time(self, ds):
-        """Get scanline acquisition time for the given dataset."""
+        """Get scanline acquisition time for the given dataset.
+
+        Note that the acquisition time does not increase monotonically
+        with the scanline number due to the scan pattern and rectification.
+        """
         # Variable is sometimes named "time" and sometimes "time_ir_wv".
         try:
             time2d = self.nc['time_ir_wv']
@@ -480,13 +488,13 @@ def _get_ssp(self, coord):
             sub_lonlat = np.nan
         return sub_lonlat
 
-    def _interp_angles(self, angles, name):
+    def _interp_angles(self, angles, dataset_id):
         """Get angle dataset.
 
         Files provide angles (solar/satellite zenith & azimuth) at a coarser
         resolution. Interpolate them to the desired resolution.
         """
-        if name.endswith('_vis'):
+        if self._is_high_resol(dataset_id):
             target_x = self.nc.coords['x']
             target_y = self.nc.coords['y']
         else:
@@ -494,7 +502,7 @@ def _interp_angles(self, angles, name):
             target_y = self.nc.coords['y_ir_wv']
         return self._interp_angles_cached(
             angles=angles,
-            nc_key=self.nc_keys[name],
+            nc_key=self.nc_keys[dataset_id['name']],
             target_x=target_x,
             target_y=target_y
         )
@@ -537,6 +545,9 @@ def _interp_tiepoints(self, ds, target_x, target_y):
 
         return ds.interp(x=target_x.values, y=target_y.values)
 
+    def _is_high_resol(self, dataset_id):
+        return dataset_id['resolution'] == HIGH_RESOL
+
 
 class FiduceoMviriEasyFcdrFileHandler(FiduceoMviriBase):
     """File handler for FIDUCEO MVIRI Easy FCDR."""
@@ -606,7 +617,9 @@ def _vis_radiance_to_reflectance(self, rad):
         Reference: [PUG], equation (6).
         """
         sza = self.get_dataset(
-            make_dataid(name='solar_zenith_angle_vis'), info={}
+            DataQuery(name='solar_zenith_angle_vis',
+                      resolution=HIGH_RESOL),
+            dataset_info={}
         )
         sza = sza.where(da.fabs(sza) < 90)  # direct illumination only
         cos_sza = np.cos(np.deg2rad(sza))
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 928face502..54aea15988 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -23,10 +23,15 @@
 from unittest import mock
 import xarray as xr
 
+from pyresample.geometry import AreaDefinition
+from pyresample.utils import proj4_radius_parameters
 from satpy.tests.utils import make_dataid
 from satpy.readers.mviri_l1b_fiduceo_nc import (
     FiduceoMviriEasyFcdrFileHandler,
-    FiduceoMviriFullFcdrFileHandler
+    FiduceoMviriFullFcdrFileHandler,
+    POLE_RADIUS,
+    EQUATOR_RADIUS,
+    ALTITUDE
 )
 
 
@@ -48,10 +53,10 @@
     {'sun_earth_distance_correction_applied': True,
      'sun_earth_distance_correction_factor': 1.}
 )
-acq_time_hires_exp = [np.datetime64('1970-01-01 00:30'),
-                      np.datetime64('1970-01-01 00:30'),
-                      np.datetime64('1970-01-01 02:30'),
-                      np.datetime64('1970-01-01 02:30')]
+acq_time_vis_exp = [np.datetime64('1970-01-01 00:30'),
+                    np.datetime64('1970-01-01 00:30'),
+                    np.datetime64('1970-01-01 02:30'),
+                    np.datetime64('1970-01-01 02:30')]
 vis_counts_exp = xr.DataArray(
     [[0., 17., 34., 51.],
      [68., 85., 102., 119.],
@@ -61,7 +66,7 @@
     coords={
         'y': [1, 2, 3, 4],
         'x': [1, 2, 3, 4],
-        'acq_time': ('y', acq_time_hires_exp),
+        'acq_time': ('y', acq_time_vis_exp),
     },
     attrs=attrs_exp
 )
@@ -74,7 +79,7 @@
     coords={
         'y': [1, 2, 3, 4],
         'x': [1, 2, 3, 4],
-        'acq_time': ('y', acq_time_hires_exp),
+        'acq_time': ('y', acq_time_vis_exp),
     },
     attrs=attrs_exp
 )
@@ -90,11 +95,11 @@
     coords={
         'y': [1, 2, 3, 4],
         'x': [1, 2, 3, 4],
-        'acq_time': ('y', acq_time_hires_exp),
+        'acq_time': ('y', acq_time_vis_exp),
     },
     attrs=attrs_refl_exp
 )
-acq_time_lores_exp = [np.datetime64('1970-01-01 00:30'),
+acq_time_ir_wv_exp = [np.datetime64('1970-01-01 00:30'),
                       np.datetime64('1970-01-01 02:30')]
 wv_counts_exp = xr.DataArray(
     [[0, 85],
@@ -103,7 +108,7 @@
     coords={
         'y': [1, 2],
         'x': [1, 2],
-        'acq_time': ('y', acq_time_lores_exp),
+        'acq_time': ('y', acq_time_ir_wv_exp),
     },
     attrs=attrs_exp
 )
@@ -114,7 +119,7 @@
     coords={
         'y': [1, 2],
         'x': [1, 2],
-        'acq_time': ('y', acq_time_lores_exp),
+        'acq_time': ('y', acq_time_ir_wv_exp),
     },
     attrs=attrs_exp
 )
@@ -125,7 +130,7 @@
     coords={
         'y': [1, 2],
         'x': [1, 2],
-        'acq_time': ('y', acq_time_lores_exp),
+        'acq_time': ('y', acq_time_ir_wv_exp),
     },
     attrs=attrs_exp
 )
@@ -136,7 +141,7 @@
     coords={
         'y': [1, 2],
         'x': [1, 2],
-        'acq_time': ('y', acq_time_lores_exp),
+        'acq_time': ('y', acq_time_ir_wv_exp),
     },
     attrs=attrs_exp
 )
@@ -147,7 +152,7 @@
     coords={
         'y': [1, 2],
         'x': [1, 2],
-        'acq_time': ('y', acq_time_lores_exp),
+        'acq_time': ('y', acq_time_ir_wv_exp),
     },
     attrs=attrs_exp
 )
@@ -158,7 +163,7 @@
     coords={
         'y': [1, 2],
         'x': [1, 2],
-        'acq_time': ('y', acq_time_lores_exp),
+        'acq_time': ('y', acq_time_ir_wv_exp),
     },
     attrs=attrs_exp
 )
@@ -196,6 +201,27 @@
     },
     attrs=attrs_exp
 )
+area_vis_exp = AreaDefinition(
+    area_id='geos_mviri_vis',
+    proj_id='geos_mviri_vis',
+    description='MVIRI Geostationary Projection',
+    projection={
+        'proj': 'geos',
+        'lon_0': 57.0,
+        'h': ALTITUDE,
+        'a': EQUATOR_RADIUS,
+        'b': POLE_RADIUS
+    },
+    width=5000,
+    height=5000,
+    area_extent=[5621229.74392, 5621229.74392, -5621229.74392, -5621229.74392]
+)
+area_ir_wv_exp = area_vis_exp.copy(
+    area_id='geos_mviri_ir_wv',
+    proj_id='geos_mviri_ir_wv',
+    width=2500,
+    height=2500
+)
 
 
 @pytest.fixture()
@@ -295,29 +321,30 @@ def test_init(self, file_handler):
         assert file_handler.mask_bad_quality is True
 
     @pytest.mark.parametrize(
-        ('name', 'calibration', 'expected'),
+        ('name', 'calibration', 'resolution', 'expected'),
         [
-            ('VIS', 'counts', vis_counts_exp),
-            ('VIS', 'radiance', vis_rad_exp),
-            ('VIS', 'reflectance', vis_refl_exp),
-            ('WV', 'counts', wv_counts_exp),
-            ('WV', 'radiance', wv_rad_exp),
-            ('WV', 'brightness_temperature', wv_bt_exp),
-            ('IR', 'counts', ir_counts_exp),
-            ('IR', 'radiance', ir_rad_exp),
-            ('IR', 'brightness_temperature', ir_bt_exp),
-            ('quality_pixel_bitmask', None, quality_pixel_bitmask_exp),
-            ('solar_zenith_angle_vis', None, sza_vis_exp),
-            ('solar_zenith_angle_ir_wv', None, sza_ir_wv_exp)
+            ('VIS', 'counts', 2250, vis_counts_exp),
+            ('VIS', 'radiance', 2250, vis_rad_exp),
+            ('VIS', 'reflectance', 2250, vis_refl_exp),
+            ('WV', 'counts', 4500, wv_counts_exp),
+            ('WV', 'radiance', 4500, wv_rad_exp),
+            ('WV', 'brightness_temperature', 4500, wv_bt_exp),
+            ('IR', 'counts', 4500, ir_counts_exp),
+            ('IR', 'radiance', 4500, ir_rad_exp),
+            ('IR', 'brightness_temperature', 4500, ir_bt_exp),
+            ('quality_pixel_bitmask', None, 2250, quality_pixel_bitmask_exp),
+            ('solar_zenith_angle_vis', None, 2250, sza_vis_exp),
+            ('solar_zenith_angle_ir_wv', None, 4500, sza_ir_wv_exp)
         ]
     )
-    def test_get_dataset(self, file_handler, name, calibration, expected):
+    def test_get_dataset(self, file_handler, name, calibration, resolution,
+                         expected):
         """Test getting datasets."""
-        id_keys = {'name': name}
+        id_keys = {'name': name, 'resolution': resolution}
         if calibration:
             id_keys['calibration'] = calibration
-        data_id = make_dataid(**id_keys)
-        info = {'platform': 'MET7'}
+        dataset_id = make_dataid(**id_keys)
+        dataset_info = {'platform': 'MET7'}
 
         is_easy = isinstance(file_handler, FiduceoMviriEasyFcdrFileHandler)
         is_vis = name == 'VIS'
@@ -325,9 +352,9 @@ def test_get_dataset(self, file_handler, name, calibration, expected):
         if is_easy and is_vis and not is_refl:
             # VIS counts/radiance not available in easy FCDR
             with pytest.raises(ValueError):
-                file_handler.get_dataset(data_id, info)
+                file_handler.get_dataset(dataset_id, dataset_info)
         else:
-            ds = file_handler.get_dataset(data_id, info)
+            ds = file_handler.get_dataset(dataset_id, dataset_info)
             xr.testing.assert_allclose(ds, expected)
             assert ds.attrs == expected.attrs
 
@@ -390,3 +417,27 @@ def test_angle_cache(self, file_handler):
         )
         with pytest.raises(AssertionError):
             xr.testing.assert_equal(sza3, sza1)
+
+    @pytest.mark.parametrize(
+        ('name', 'resolution', 'area_exp'),
+        [
+            ('VIS', 2250, area_vis_exp),
+            ('WV', 4500, area_ir_wv_exp),
+            ('IR', 4500, area_ir_wv_exp),
+            ('quality_pixel_bitmask', 4500, area_vis_exp),
+            ('solar_zenith_angle_vis', 4500, area_vis_exp),
+            ('solar_zenith_angle_ir_wv', 2250, area_ir_wv_exp)
+        ]
+    )
+    def test_get_area_definition(self, file_handler, name, resolution,
+                                 area_exp):
+        """Test getting area definitions."""
+        dataset_id = make_dataid(name=name, resolution=resolution)
+        area = file_handler.get_area_def(dataset_id)
+        a, b = proj4_radius_parameters(area.proj_dict)
+        a_exp, b_exp = proj4_radius_parameters(area_exp.proj_dict)
+        assert a == a_exp
+        assert b == b_exp
+        for key in ['h', 'lon_0', 'proj', 'units']:
+            assert area.proj_dict[key] == area_exp.proj_dict[key]
+        np.testing.assert_allclose(area.area_extent, area_exp.area_extent)

From e9f7f0ef37f2e36a61f99b1c191e18cdaee46272 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 19 Nov 2020 16:04:18 +0000
Subject: [PATCH 18/47] Make sure returned datasets are float32

---
 satpy/readers/mviri_l1b_fiduceo_nc.py         |  27 ++--
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py | 131 ++++++++++++------
 2 files changed, 105 insertions(+), 53 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index f1923441ae..b34d4d8a32 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -372,8 +372,10 @@ def _ir_wv_counts_to_radiance(self, counts, channel):
             a, b = self.nc['a_wv'], self.nc['b_wv']
         else:
             a, b = self.nc['a_ir'], self.nc['b_ir']
+        a = np.float32(a)
+        b = np.float32(b)
         rad = a + b * counts
-        return rad.where(rad > 0)
+        return rad.where(rad > 0, np.float32(np.nan))
 
     def _ir_wv_radiance_to_brightness_temperature(self, rad, channel):
         """Convert IR/WV radiance to brightness temperature.
@@ -384,8 +386,10 @@ def _ir_wv_radiance_to_brightness_temperature(self, rad, channel):
             a, b = self.nc['bt_a_wv'], self.nc['bt_b_wv']
         else:
             a, b = self.nc['bt_a_ir'], self.nc['bt_b_ir']
+        a = np.float32(a)
+        b = np.float32(b)
         bt = b / (np.log(rad) - a)
-        return bt.where(bt > 0)
+        return bt.where(bt > 0, np.float32(np.nan))
 
     def _mask_vis(self, ds):
         """Mask VIS pixels with bad quality.
@@ -396,7 +400,8 @@ def _mask_vis(self, ds):
         the VIS channel.
         """
         mask = self.nc['quality_pixel_bitmask']
-        return ds.where(np.logical_or(mask == 0, mask == 2))
+        return ds.where(np.logical_or(mask == 0, mask == 2),
+                        np.float32(np.nan))
 
     def _get_acq_time(self, ds):
         """Get scanline acquisition time for the given dataset.
@@ -603,9 +608,10 @@ def _vis_counts_to_radiance(self, counts):
         a_cf = (self.nc['a0_vis'] +
                 self.nc['a1_vis'] * years_since_launch +
                 self.nc['a2_vis'] * years_since_launch ** 2)
-
-        rad = (counts - self.nc['mean_count_space_vis']) * a_cf
-        return rad.where(rad > 0)
+        mean_count_space_vis = np.float32(self.nc['mean_count_space_vis'])
+        a_cf = np.float32(a_cf)
+        rad = (counts - mean_count_space_vis) * a_cf
+        return rad.where(rad > 0, np.float32(np.nan))
 
     def _vis_radiance_to_reflectance(self, rad):
         """Convert VIS radiance to reflectance factor.
@@ -621,11 +627,14 @@ def _vis_radiance_to_reflectance(self, rad):
                       resolution=HIGH_RESOL),
             dataset_info={}
         )
-        sza = sza.where(da.fabs(sza) < 90)  # direct illumination only
+        sza = sza.where(da.fabs(sza) < 90,
+                        np.float32(np.nan))  # direct illumination only
         cos_sza = np.cos(np.deg2rad(sza))
+        distance_sun_earth2 = np.float32(self.nc['distance_sun_earth'] ** 2)
+        solar_irradiance_vis = np.float32(self.nc['solar_irradiance_vis'])
         refl = (
-           (np.pi * self.nc['distance_sun_earth'] ** 2) /
-           (self.nc['solar_irradiance_vis'] * cos_sza) *
+           (np.pi * distance_sun_earth2) /
+           (solar_irradiance_vis * cos_sza) *
            rad
         )
         refl = refl * 100  # conversion to percent
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 54aea15988..4236b7ee5f 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -58,10 +58,13 @@
                     np.datetime64('1970-01-01 02:30'),
                     np.datetime64('1970-01-01 02:30')]
 vis_counts_exp = xr.DataArray(
-    [[0., 17., 34., 51.],
-     [68., 85., 102., 119.],
-     [136., 153., np.nan, 187.],
-     [204., 221., 238., 255]],
+    np.array(
+        [[0., 17., 34., 51.],
+         [68., 85., 102., 119.],
+         [136., 153., np.nan, 187.],
+         [204., 221., 238., 255]],
+        dtype=np.float32
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2, 3, 4],
@@ -71,10 +74,13 @@
     attrs=attrs_exp
 )
 vis_rad_exp = xr.DataArray(
-    [[np.nan, 18.56, 38.28, 58.],
-     [77.72, 97.44, 117.16, 136.88],
-     [156.6, 176.32, np.nan, 215.76],
-     [235.48, 255.2, 274.92, 294.64]],
+    np.array(
+        [[np.nan, 18.56, 38.28, 58.],
+         [77.72, 97.44, 117.16, 136.88],
+         [156.6, 176.32, np.nan, 215.76],
+         [235.48, 255.2, 274.92, 294.64]],
+        dtype=np.float32
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2, 3, 4],
@@ -84,10 +90,13 @@
     attrs=attrs_exp
 )
 vis_refl_exp = xr.DataArray(
-    [[np.nan, 23.440929, np.nan, np.nan],
-     [40.658744, 66.602233, 147.970867, np.nan],
-     [75.688217, 92.240733, np.nan, np.nan],
-     [np.nan, np.nan, np.nan, np.nan]],
+    np.array(
+        [[np.nan, 23.440929, np.nan, np.nan],
+         [40.658744, 66.602233, 147.970867, np.nan],
+         [75.688217, 92.240733, np.nan, np.nan],
+         [np.nan, np.nan, np.nan, np.nan]],
+        dtype=np.float32
+    ),
     # (0, 0) and (2, 2) are NaN because radiance is NaN
     # (0, 2) is NaN because SZA >= 90 degrees
     # Last row/col is NaN due to SZA interpolation
@@ -102,8 +111,11 @@
 acq_time_ir_wv_exp = [np.datetime64('1970-01-01 00:30'),
                       np.datetime64('1970-01-01 02:30')]
 wv_counts_exp = xr.DataArray(
-    [[0, 85],
-     [170, 255]],
+    np.array(
+        [[0, 85],
+         [170, 255]],
+        dtype=np.uint8
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2],
@@ -113,8 +125,11 @@
     attrs=attrs_exp
 )
 wv_rad_exp = xr.DataArray(
-    [[np.nan, 3.75],
-     [8, 12.25]],
+    np.array(
+        [[np.nan, 3.75],
+         [8, 12.25]],
+        dtype=np.float32
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2],
@@ -124,8 +139,11 @@
     attrs=attrs_exp
 )
 wv_bt_exp = xr.DataArray(
-    [[np.nan, 230.461366],
-     [252.507448, 266.863289]],
+    np.array(
+        [[np.nan, 230.461366],
+         [252.507448, 266.863289]],
+        dtype=np.float32
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2],
@@ -135,8 +153,11 @@
     attrs=attrs_exp
 )
 ir_counts_exp = xr.DataArray(
-    [[0, 85],
-     [170, 255]],
+    np.array(
+        [[0, 85],
+         [170, 255]],
+        dtype=np.uint8
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2],
@@ -146,8 +167,11 @@
     attrs=attrs_exp
 )
 ir_rad_exp = xr.DataArray(
-    [[np.nan, 80],
-     [165, 250]],
+    np.array(
+        [[np.nan, 80],
+         [165, 250]],
+        dtype=np.float32
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2],
@@ -157,8 +181,11 @@
     attrs=attrs_exp
 )
 ir_bt_exp = xr.DataArray(
-    [[np.nan, 178.00013189],
-     [204.32955838, 223.28709913]],
+    np.array(
+        [[np.nan, 178.00013189],
+         [204.32955838, 223.28709913]],
+        dtype=np.float32
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2],
@@ -168,10 +195,13 @@
     attrs=attrs_exp
 )
 quality_pixel_bitmask_exp = xr.DataArray(
-    [[2, 0, 0, 0],
-     [0, 0, 0, 0],
-     [0, 0, 1, 0],
-     [0, 0, 0, 0]],
+    np.array(
+        [[2, 0, 0, 0],
+         [0, 0, 0, 0],
+         [0, 0, 1, 0],
+         [0, 0, 0, 0]],
+        dtype=np.uint8
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2, 3, 4],
@@ -180,10 +210,13 @@
     attrs=attrs_exp
 )
 sza_vis_exp = xr.DataArray(
-    [[45., 67.5, 90., np.nan],
-     [22.5, 45., 67.5, np.nan],
-     [0., 22.5, 45., np.nan],
-     [np.nan, np.nan, np.nan, np.nan]],
+    np.array(
+        [[45., 67.5, 90., np.nan],
+         [22.5, 45., 67.5, np.nan],
+         [0., 22.5, 45., np.nan],
+         [np.nan, np.nan, np.nan, np.nan]],
+        dtype=np.float32
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2, 3, 4],
@@ -192,8 +225,11 @@
     attrs=attrs_exp
 )
 sza_ir_wv_exp = xr.DataArray(
-    [[45, 90],
-     [0, 45]],
+    np.array(
+        [[45, 90],
+         [0, 45]],
+        dtype=np.float32
+    ),
     dims=('y', 'x'),
     coords={
         'y': [1, 2],
@@ -227,18 +263,24 @@
 @pytest.fixture()
 def fake_dataset():
     """Create fake dataset."""
-    count_ir = da.linspace(0, 255, 4, dtype=np.float32).reshape(2, 2)
-    count_wv = da.linspace(0, 255, 4, dtype=np.float32).reshape(2, 2)
-    count_vis = da.linspace(0, 255, 16, dtype=np.float32).reshape(4, 4)
+    count_ir = da.linspace(0, 255, 4, dtype=np.uint8).reshape(2, 2)
+    count_wv = da.linspace(0, 255, 4, dtype=np.uint8).reshape(2, 2)
+    count_vis = da.linspace(0, 255, 16, dtype=np.uint8).reshape(4, 4)
     sza = da.from_array(
-        [[45, 90],
-         [0, 45]]
-    ).astype(np.float32)
+        np.array(
+            [[45, 90],
+             [0, 45]],
+            dtype=np.float32
+        )
+    )
     mask = da.from_array(
-        [[2, 0, 0, 0],  # 2 = "use with caution"
-         [0, 0, 0, 0],
-         [0, 0, 1, 0],  # 1 = "invalid"
-         [0, 0, 0, 0]]
+        np.array(
+            [[2, 0, 0, 0],  # 2 = "use with caution"
+             [0, 0, 0, 0],
+             [0, 0, 1, 0],  # 1 = "invalid"
+             [0, 0, 0, 0]],
+            dtype=np.uint8
+        )
     )
     time = np.arange(4).astype('datetime64[h]').reshape(2, 2)
     ds = xr.Dataset(
@@ -356,6 +398,7 @@ def test_get_dataset(self, file_handler, name, calibration, resolution,
         else:
             ds = file_handler.get_dataset(dataset_id, dataset_info)
             xr.testing.assert_allclose(ds, expected)
+            assert ds.dtype == expected.dtype
             assert ds.attrs == expected.attrs
 
     def test_time_cache(self, file_handler):

From 48716a296cde87fd29effba737944a2fef9351df Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 19 Nov 2020 16:05:11 +0000
Subject: [PATCH 19/47] Add test for ang2fac

---
 satpy/tests/reader_tests/test_geos_area.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/satpy/tests/reader_tests/test_geos_area.py b/satpy/tests/reader_tests/test_geos_area.py
index 53811e6500..f4bf663962 100644
--- a/satpy/tests/reader_tests/test_geos_area.py
+++ b/satpy/tests/reader_tests/test_geos_area.py
@@ -20,7 +20,8 @@
 import unittest
 from satpy.readers._geos_area import (get_xy_from_linecol,
                                       get_area_extent,
-                                      get_area_definition)
+                                      get_area_definition,
+                                      ang2fac)
 
 import numpy as np
 
@@ -145,3 +146,9 @@ def test_get_area_definition(self):
         self.assertEqual(a, 6378169)
         self.assertEqual(b, 6356583.8)
         self.assertEqual(a_def.proj_dict['h'], 35785831)
+
+    def test_ang2fac(self):
+        """Test conversion from angular sampling to line/column offset."""
+        lfac = 13642337  # SEVIRI LFAC
+        sampling = np.deg2rad(2 ** 16 / lfac)
+        np.testing.assert_allclose(ang2fac(sampling), lfac)

From 61bd85a7f1f631addae2c1c40fdc3bbd043e7939 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 19 Nov 2020 16:06:14 +0000
Subject: [PATCH 20/47] Fix angle resolution

---
 satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
index 2ed3d7063f..86b6d5923d 100644
--- a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
+++ b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
@@ -139,7 +139,7 @@ datasets:
     standard_name: solar_zenith_angle
     long_name: "Solar zenith angle IR/WV channels"
     units: degree
-    resolution: 2250
+    resolution: 4500
     file_type: [nc_easy, nc_full]
 
   solar_azimuth_angle_ir_wv:
@@ -147,7 +147,7 @@ datasets:
     standard_name: solar_azimuth_angle
     long_name: "Solar azimuth angle IR/WV channels"
     units: degree
-    resolution: 2250
+    resolution: 4500
     file_type: [nc_easy, nc_full]
 
   satellite_zenith_angle_ir_wv:
@@ -155,7 +155,7 @@ datasets:
     standard_name: sensor_zenith_angle
     long_name: "Satellite zenith angle IR/WV channels"
     units: degree
-    resolution: 2250
+    resolution: 4500
     file_type: [nc_easy, nc_full]
 
   satellite_azimuth_angle_ir_wv:
@@ -163,5 +163,5 @@ datasets:
     standard_name: sensor_azimuth_angle
     long_name: "Satellite azimuth angle IR/WV channels"
     units: degree
-    resolution: 2250
+    resolution: 4500
     file_type: [nc_easy, nc_full]
\ No newline at end of file

From 0438f2a68f4e71305beec238f0b7da760506ad83 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 19 Nov 2020 16:13:34 +0000
Subject: [PATCH 21/47] Refactor IR/WV coefficient fetching

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 20 ++++++++------------
 1 file changed, 8 insertions(+), 12 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index b34d4d8a32..327d0923b1 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -368,12 +368,10 @@ def _ir_wv_counts_to_radiance(self, counts, channel):
 
         Reference: [PUG], equations (4.1) and (4.2).
         """
-        if channel == 'WV':
-            a, b = self.nc['a_wv'], self.nc['b_wv']
-        else:
-            a, b = self.nc['a_ir'], self.nc['b_ir']
-        a = np.float32(a)
-        b = np.float32(b)
+        nc_key_a = self.nc_keys_coefs[channel]['radiance']['a']
+        nc_key_b = self.nc_keys_coefs[channel]['radiance']['b']
+        a = np.float32(self.nc[nc_key_a])
+        b = np.float32(self.nc[nc_key_b])
         rad = a + b * counts
         return rad.where(rad > 0, np.float32(np.nan))
 
@@ -382,12 +380,10 @@ def _ir_wv_radiance_to_brightness_temperature(self, rad, channel):
 
         Reference: [PUG], equations (5.1) and (5.2).
         """
-        if channel == 'WV':
-            a, b = self.nc['bt_a_wv'], self.nc['bt_b_wv']
-        else:
-            a, b = self.nc['bt_a_ir'], self.nc['bt_b_ir']
-        a = np.float32(a)
-        b = np.float32(b)
+        nc_key_a = self.nc_keys_coefs[channel]['brightness_temperature']['a']
+        nc_key_b = self.nc_keys_coefs[channel]['brightness_temperature']['b']
+        a = np.float32(self.nc[nc_key_a])
+        b = np.float32(self.nc[nc_key_b])
         bt = b / (np.log(rad) - a)
         return bt.where(bt > 0, np.float32(np.nan))
 

From 6759ad852c3c2bc509f3197d98c5a1553e4418ae Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 19 Nov 2020 16:24:47 +0000
Subject: [PATCH 22/47] Test file pattern matching

---
 satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml   |  6 +--
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py | 37 ++++++++++++++-----
 2 files changed, 30 insertions(+), 13 deletions(-)

diff --git a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
index 86b6d5923d..0fbbb3ba39 100644
--- a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
+++ b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
@@ -19,13 +19,13 @@ file_types:
     nc_easy:
         file_reader: !!python/name:satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriEasyFcdrFileHandler
         file_patterns: [
-          'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_EASY_{processor_version}_{format_version}.nc'
+          'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude:f}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_EASY_{processor_version}_{format_version}.nc'
           # Example: FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc
         ]
     nc_full:
         file_reader: !!python/name:satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriFullFcdrFileHandler
         file_patterns: [
-          'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_FULL_{processor_version}_{format_version}.nc'
+          'FIDUCEO_FCDR_{level}_{sensor}_{platform}-{projection_longitude:f}_{start_time:%Y%m%d%H%M}_{end_time:%Y%m%d%H%M}_FULL_{processor_version}_{format_version}.nc'
           # Example: FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_FULL_v2.6_fv3.1.nc
         ]
 
@@ -164,4 +164,4 @@ datasets:
     long_name: "Satellite azimuth angle IR/WV channels"
     units: degree
     resolution: 4500
-    file_type: [nc_easy, nc_full]
\ No newline at end of file
+    file_type: [nc_easy, nc_full]
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 4236b7ee5f..2e970cfdc8 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -19,6 +19,7 @@
 
 import dask.array as da
 import numpy as np
+import os
 import pytest
 from unittest import mock
 import xarray as xr
@@ -344,16 +345,6 @@ def file_handler(fake_dataset, request):
         )
 
 
-
-
-
-# TODO: Test file patterns
-
-
-
-
-
-
 class TestFiduceoMviriFileHandlers:
     """Unit tests for FIDUCEO MVIRI file handlers."""
 
@@ -484,3 +475,29 @@ def test_get_area_definition(self, file_handler, name, resolution,
         for key in ['h', 'lon_0', 'proj', 'units']:
             assert area.proj_dict[key] == area_exp.proj_dict[key]
         np.testing.assert_allclose(area.area_extent, area_exp.area_extent)
+
+
+@pytest.fixture
+def reader():
+    """Return MVIRI FIDUCEO FCDR reader."""
+    from satpy.config import config_search_paths
+    from satpy.readers import load_reader
+
+    reader_configs = config_search_paths(
+        os.path.join("readers", "mviri_l1b_fiduceo_nc.yaml"))
+    reader = load_reader(reader_configs)
+    return reader
+
+
+def test_file_pattern(reader):
+    """Test file pattern matching."""
+    filenames = [
+        "FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_FULL_v2.6_fv3.1.nc",
+        "FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc",
+        "FIDUCEO_FCDR_L15_MVIRI_MET7-00.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc",
+        "abcde",
+    ]
+
+    files = reader.select_files_from_pathnames(filenames)
+    # only 3 out of 4 above should match
+    assert len(files) == 3

From 62bfd931cb27073b0f1b59c0ab7ad61a8a62b185 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 19 Nov 2020 16:36:32 +0000
Subject: [PATCH 23/47] Add MVIRI reader to overview table

---
 doc/source/index.rst | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/doc/source/index.rst b/doc/source/index.rst
index 7f96b93106..3ec9dafdc1 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -93,6 +93,9 @@ the base Satpy installation.
     * - MSG (Meteosat 8 to 11) L2 products in GRIB2 format
       - `seviri_l2_grib`
       - In development, CLM, OCA and FIR products supported
+    * - MFG (Meteosat 2 to 7) MVIRI data in netCDF format (FIDUCEO FCDR)
+      - `mviri_l1b_fiduceo_nc`
+      - Beta
     * - Himawari 8 and 9 AHI data in HSD format
       - `ahi_hsd`
       - Nominal

From bb90e340bb4e32428654290d2af4bfefa4e45fff Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 20 Nov 2020 10:05:22 +0000
Subject: [PATCH 24/47] Factorize coefficient fetching

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 22 ++++++++++++++--------
 1 file changed, 14 insertions(+), 8 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 327d0923b1..21879846ef 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -363,15 +363,24 @@ def _calibrate_ir_wv(self, ds, channel, calibration):
         else:
             raise KeyError('Invalid calibration: {}'.format(calibration.name))
 
+    def _get_coefs_ir_wv(self, channel, calibration):
+        """Get calibration coefficients for IR/WV channels.
+
+        Returns:
+            Offset (a), Slope (b)
+        """
+        nc_key_a = self.nc_keys_coefs[channel][calibration]['a']
+        nc_key_b = self.nc_keys_coefs[channel][calibration]['b']
+        a = np.float32(self.nc[nc_key_a])
+        b = np.float32(self.nc[nc_key_b])
+        return a, b
+
     def _ir_wv_counts_to_radiance(self, counts, channel):
         """Convert IR/WV counts to radiance.
 
         Reference: [PUG], equations (4.1) and (4.2).
         """
-        nc_key_a = self.nc_keys_coefs[channel]['radiance']['a']
-        nc_key_b = self.nc_keys_coefs[channel]['radiance']['b']
-        a = np.float32(self.nc[nc_key_a])
-        b = np.float32(self.nc[nc_key_b])
+        a, b = self._get_coefs_ir_wv(channel, 'radiance')
         rad = a + b * counts
         return rad.where(rad > 0, np.float32(np.nan))
 
@@ -380,10 +389,7 @@ def _ir_wv_radiance_to_brightness_temperature(self, rad, channel):
 
         Reference: [PUG], equations (5.1) and (5.2).
         """
-        nc_key_a = self.nc_keys_coefs[channel]['brightness_temperature']['a']
-        nc_key_b = self.nc_keys_coefs[channel]['brightness_temperature']['b']
-        a = np.float32(self.nc[nc_key_a])
-        b = np.float32(self.nc[nc_key_b])
+        a, b = self._get_coefs_ir_wv(channel, 'brightness_temperature')
         bt = b / (np.log(rad) - a)
         return bt.where(bt > 0, np.float32(np.nan))
 

From abb53b2d9d43868887fd8d5006a19b0a2c12cacf Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 20 Nov 2020 12:49:28 +0000
Subject: [PATCH 25/47] Resole a couple of test misses

---
 satpy/readers/mviri_l1b_fiduceo_nc.py         | 12 ++---
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py | 50 ++++++++++++++-----
 2 files changed, 42 insertions(+), 20 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 21879846ef..85267c9bf2 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -110,17 +110,17 @@
 """
 
 import abc
-import dask.array as da
 import functools
+
+import dask.array as da
 import numpy as np
 import xarray as xr
 
 from satpy import CHUNK_SIZE
-from satpy.readers.file_handlers import BaseFileHandler
+from satpy.dataset.dataid import DataQuery
 from satpy.readers._geos_area import (ang2fac, get_area_definition,
                                       get_area_extent)
-from satpy.dataset.dataid import DataQuery
-
+from satpy.readers.file_handlers import BaseFileHandler
 
 EQUATOR_RADIUS = 6378140.0
 POLE_RADIUS = 6356755.0
@@ -270,10 +270,6 @@ def get_dataset(self, dataset_id, dataset_info):
         self._update_attrs(ds, dataset_info)
         return ds
 
-    def _get_nc_key(self, name):
-        """Get netCDF key corresponding to the given dataset name."""
-        return
-
     def _read_dataset(self, name):
         """Read a dataset from the file."""
         nc_key = self.nc_keys.get(name, name)
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 2e970cfdc8..0d3da0d54a 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -17,24 +17,20 @@
 # satpy.  If not, see <http://www.gnu.org/licenses/>.
 """Unit tests for the FIDUCEO MVIRI FCDR Reader."""
 
+import os
+from unittest import mock
+
 import dask.array as da
 import numpy as np
-import os
 import pytest
-from unittest import mock
 import xarray as xr
-
 from pyresample.geometry import AreaDefinition
 from pyresample.utils import proj4_radius_parameters
-from satpy.tests.utils import make_dataid
-from satpy.readers.mviri_l1b_fiduceo_nc import (
-    FiduceoMviriEasyFcdrFileHandler,
-    FiduceoMviriFullFcdrFileHandler,
-    POLE_RADIUS,
-    EQUATOR_RADIUS,
-    ALTITUDE
-)
 
+from satpy.readers.mviri_l1b_fiduceo_nc import (
+    ALTITUDE, EQUATOR_RADIUS, POLE_RADIUS, FiduceoMviriEasyFcdrFileHandler,
+    FiduceoMviriFullFcdrFileHandler)
+from satpy.tests.utils import make_dataid
 
 attrs_exp = {
     'platform': 'MET7',
@@ -109,6 +105,21 @@
     },
     attrs=attrs_refl_exp
 )
+u_vis_refl_exp = xr.DataArray(
+    np.array(
+        [[0.1, 0.2, 0.3, 0.4],
+         [0.5, 0.6, 0.7, 0.8],
+         [0.9, 1.0, 1.1, 1.2],
+         [1.3, 1.4, 1.5, 1.6]],
+        dtype=np.float32
+    ),
+    dims=('y', 'x'),
+    coords={
+        'y': [1, 2, 3, 4],
+        'x': [1, 2, 3, 4],
+    },
+    attrs=attrs_exp
+)
 acq_time_ir_wv_exp = [np.datetime64('1970-01-01 00:30'),
                       np.datetime64('1970-01-01 02:30')]
 wv_counts_exp = xr.DataArray(
@@ -291,6 +302,8 @@ def fake_dataset():
             'count_ir': (('y_ir_wv', 'x_ir_wv'), count_ir),
             'toa_bidirectional_reflectance_vis': (
                 ('y', 'x'), vis_refl_exp / 100),
+            'u_independent_toa_bidirectional_reflectance': (
+                ('y', 'x'), u_vis_refl_exp / 100),
             'quality_pixel_bitmask': (('y', 'x'), mask),
             'solar_zenith_angle': (('y_tie', 'x_tie'), sza),
             'time_ir_wv': (('y_ir_wv', 'x_ir_wv'), time),
@@ -367,7 +380,8 @@ def test_init(self, file_handler):
             ('IR', 'brightness_temperature', 4500, ir_bt_exp),
             ('quality_pixel_bitmask', None, 2250, quality_pixel_bitmask_exp),
             ('solar_zenith_angle_vis', None, 2250, sza_vis_exp),
-            ('solar_zenith_angle_ir_wv', None, 4500, sza_ir_wv_exp)
+            ('solar_zenith_angle_ir_wv', None, 4500, sza_ir_wv_exp),
+            ('u_independent_toa_bidirectional_reflectance', None, 4500, u_vis_refl_exp)
         ]
     )
     def test_get_dataset(self, file_handler, name, calibration, resolution,
@@ -476,6 +490,18 @@ def test_get_area_definition(self, file_handler, name, resolution,
             assert area.proj_dict[key] == area_exp.proj_dict[key]
         np.testing.assert_allclose(area.area_extent, area_exp.area_extent)
 
+    def test_calib_exceptions(self, file_handler):
+        """Test calibration exceptions."""
+        with pytest.raises(KeyError):
+            file_handler.calibrate(None, 'solar_zenith_angle', None)
+
+        calib = mock.MagicMock()
+        calib.configure_mock(name='invalid_calib')
+        with pytest.raises(KeyError):
+            file_handler.calibrate(None, 'VIS', calib)
+        with pytest.raises(KeyError):
+            file_handler.calibrate(None, 'IR', calib)
+
 
 @pytest.fixture
 def reader():

From 84cff724fa98343cf96bcc6dfc574a60d2be7417 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 20 Nov 2020 14:46:14 +0000
Subject: [PATCH 26/47] Simplify angle datasets

---
 satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml   | 64 ++++------------
 satpy/readers/mviri_l1b_fiduceo_nc.py         | 73 ++++++++++++-------
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py | 20 ++---
 3 files changed, 75 insertions(+), 82 deletions(-)

diff --git a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
index 0fbbb3ba39..bf17a427ee 100644
--- a/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
+++ b/satpy/etc/readers/mviri_l1b_fiduceo_nc.yaml
@@ -102,66 +102,34 @@ datasets:
     resolution: 2250
     file_type: [nc_easy]
 
-  solar_zenith_angle_vis:
-    name: solar_zenith_angle_vis
+  solar_zenith_angle:
+    name: solar_zenith_angle
     standard_name: solar_zenith_angle
-    long_name: "Solar zenith angle VIS channel"
+    long_name: "Solar zenith angle"
     units: degree
-    resolution: 2250
-    file_type: [nc_easy, nc_full]
-
-  solar_azimuth_angle_vis:
-    name: solar_azimuth_angle_vis
-    standard_name: solar_azimuth_angle
-    long_name: "Solar azimuth angle VIS channel"
-    units: degree
-    resolution: 2250
-    file_type: [nc_easy, nc_full]
-
-  satellite_zenith_angle_vis:
-    name: satellite_zenith_angle_vis
-    standard_name: sensor_zenith_angle
-    long_name: "Satellite zenith angle VIS channel"
-    units: degree
-    resolution: 2250
-    file_type: [nc_easy, nc_full]
-
-  satellite_azimuth_angle_vis:
-    name: satellite_azimuth_angle_vis
-    standard_name: sensor_azimuth_angle
-    long_name: "Satellite azimuth angle VIS channel"
-    units: degree
-    resolution: 2250
-    file_type: [nc_easy, nc_full]
-
-  solar_zenith_angle_ir_wv:
-    name: solar_zenith_angle_ir_wv
-    standard_name: solar_zenith_angle
-    long_name: "Solar zenith angle IR/WV channels"
-    units: degree
-    resolution: 4500
+    resolution: [2250, 4500]
     file_type: [nc_easy, nc_full]
 
-  solar_azimuth_angle_ir_wv:
-    name: solar_azimuth_angle_ir_wv
+  solar_azimuth_angle:
+    name: solar_azimuth_angle
     standard_name: solar_azimuth_angle
-    long_name: "Solar azimuth angle IR/WV channels"
+    long_name: "Solar azimuth angle"
     units: degree
-    resolution: 4500
+    resolution: [2250, 4500]
     file_type: [nc_easy, nc_full]
 
-  satellite_zenith_angle_ir_wv:
-    name: satellite_zenith_angle_ir_wv
+  satellite_zenith_angle:
+    name: satellite_zenith_angle
     standard_name: sensor_zenith_angle
-    long_name: "Satellite zenith angle IR/WV channels"
+    long_name: "Satellite zenith angle"
     units: degree
-    resolution: 4500
+    resolution: [2250, 4500]
     file_type: [nc_easy, nc_full]
 
-  satellite_azimuth_angle_ir_wv:
-    name: satellite_azimuth_angle_ir_wv
+  satellite_azimuth_angle:
+    name: satellite_azimuth_angle
     standard_name: sensor_azimuth_angle
-    long_name: "Satellite azimuth angle IR/WV channels"
+    long_name: "Satellite azimuth angle"
     units: degree
-    resolution: 4500
+    resolution: [2250, 4500]
     file_type: [nc_easy, nc_full]
diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 85267c9bf2..fc4e127e32 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -98,6 +98,36 @@
 uncertainties can be used to filter these cases before calculating reflectances.
 
 
+Angles
+------
+The FIDUCEO MVIRI FCDR provides satellite and solar angles on a coarse tiepoint
+grid. By default these datasets will be interpolated to the higher VIS
+resolution. This can be changed as follows:
+
+.. code-block:: python
+
+    scn.load(['solar_zenith_angle'], resolution=4500)
+
+If you need the angles in both resolutions, use data queries:
+
+.. code-block:: python
+
+    from satpy import DataQuery
+
+    q_vis = DataQuery(
+        name='solar_zenith_angle',
+        resolution=2250
+    )
+    q_ir = DataQuery(
+        name='solar_zenith_angle',
+        resolution=4500
+    )
+    scn.load([query_vis, query_ir])
+
+    # Use the query objects to access the datasets as follows
+    sza_vis = scn[query_vis]
+
+
 References
 ----------
     - `[Handbook]`_ MFG User Handbook
@@ -132,14 +162,10 @@
 
 CHANNELS = ['VIS', 'WV', 'IR']
 ANGLES = [
-    'solar_zenith_angle_vis',
-    'solar_azimuth_angle_vis',
-    'satellite_zenith_angle_vis',
-    'satellite_azimuth_angle_vis',
-    'solar_zenith_angle_ir_wv',
-    'solar_azimuth_angle_ir_wv',
-    'satellite_zenith_angle_ir_wv',
-    'satellite_azimuth_angle_ir_wv'
+    'solar_zenith_angle',
+    'solar_azimuth_angle',
+    'satellite_zenith_angle',
+    'satellite_azimuth_angle'
 ]
 OTHER_REFLECTANCES = [
     'u_independent_toa_bidirectional_reflectance',
@@ -195,15 +221,7 @@ class FiduceoMviriBase(BaseFileHandler):
     """Baseclass for FIDUCEO MVIRI file handlers."""
     nc_keys = {
         'WV': 'count_wv',
-        'IR': 'count_ir',
-        'solar_zenith_angle_vis': 'solar_zenith_angle',
-        'solar_azimuth_angle_vis': 'solar_azimuth_angle',
-        'satellite_zenith_angle_vis': 'satellite_zenith_angle',
-        'satellite_azimuth_angle_vis': 'satellite_azimuth_angle',
-        'solar_zenith_angle_ir_wv': 'solar_zenith_angle',
-        'solar_azimuth_angle_ir_wv': 'solar_azimuth_angle',
-        'satellite_zenith_angle_ir_wv': 'satellite_zenith_angle',
-        'satellite_azimuth_angle_ir_wv': 'satellite_azimuth_angle',
+        'IR': 'count_ir'
     }
     nc_keys_coefs = {
         'WV': {
@@ -270,9 +288,13 @@ def get_dataset(self, dataset_id, dataset_info):
         self._update_attrs(ds, dataset_info)
         return ds
 
+    def _get_nc_key(self, ds_name):
+        """Get netCDF variable name for the given dataset."""
+        return self.nc_keys.get(ds_name, ds_name)
+
     def _read_dataset(self, name):
         """Read a dataset from the file."""
-        nc_key = self.nc_keys.get(name, name)
+        nc_key = self._get_nc_key(name)
         ds = self.nc[nc_key]
         if 'y_ir_wv' in ds.dims:
             ds = ds.rename({'y_ir_wv': 'y', 'x_ir_wv': 'x'})
@@ -290,8 +312,12 @@ def _update_attrs(self, ds, info):
 
     def get_area_def(self, dataset_id):
         """Get area definition of the given dataset."""
-        ds = self._read_dataset(dataset_id['name'])
-        im_size = ds.coords['y'].size
+        if self._is_high_resol(dataset_id):
+            im_size = self.nc.coords['y'].size
+            area_name = 'geos_mviri_vis'
+        else:
+            im_size = self.nc.coords['y_ir_wv'].size
+            area_name = 'geos_mviri_ir_wv'
 
         # Determine line/column offsets and scaling factors. For offsets
         # see variables "asamp" and "aline" of subroutine "refgeo" in
@@ -300,9 +326,6 @@ def get_area_def(self, dataset_id):
         loff = coff = im_size / 2 + 0.5
         lfac = cfac = ang2fac(np.deg2rad(MVIRI_FIELD_OF_VIEW) / im_size)
 
-        area_name = 'geos_mviri_{}'.format(
-            'vis' if self._is_high_resol(dataset_id) else 'ir_wv'
-        )
         pdict = {
             'ssp_lon': self.projection_longitude,
             'a': EQUATOR_RADIUS,
@@ -505,7 +528,7 @@ def _interp_angles(self, angles, dataset_id):
             target_y = self.nc.coords['y_ir_wv']
         return self._interp_angles_cached(
             angles=angles,
-            nc_key=self.nc_keys[dataset_id['name']],
+            nc_key=self._get_nc_key(dataset_id['name']),
             target_x=target_x,
             target_y=target_y
         )
@@ -621,7 +644,7 @@ def _vis_radiance_to_reflectance(self, rad):
         Reference: [PUG], equation (6).
         """
         sza = self.get_dataset(
-            DataQuery(name='solar_zenith_angle_vis',
+            DataQuery(name='solar_zenith_angle',
                       resolution=HIGH_RESOL),
             dataset_info={}
         )
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 0d3da0d54a..c75e25261e 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -260,15 +260,15 @@
         'a': EQUATOR_RADIUS,
         'b': POLE_RADIUS
     },
-    width=5000,
-    height=5000,
+    width=4,
+    height=4,
     area_extent=[5621229.74392, 5621229.74392, -5621229.74392, -5621229.74392]
 )
 area_ir_wv_exp = area_vis_exp.copy(
     area_id='geos_mviri_ir_wv',
     proj_id='geos_mviri_ir_wv',
-    width=2500,
-    height=2500
+    width=2,
+    height=2
 )
 
 
@@ -379,8 +379,8 @@ def test_init(self, file_handler):
             ('IR', 'radiance', 4500, ir_rad_exp),
             ('IR', 'brightness_temperature', 4500, ir_bt_exp),
             ('quality_pixel_bitmask', None, 2250, quality_pixel_bitmask_exp),
-            ('solar_zenith_angle_vis', None, 2250, sza_vis_exp),
-            ('solar_zenith_angle_ir_wv', None, 4500, sza_ir_wv_exp),
+            ('solar_zenith_angle', None, 2250, sza_vis_exp),
+            ('solar_zenith_angle', None, 4500, sza_ir_wv_exp),
             ('u_independent_toa_bidirectional_reflectance', None, 4500, u_vis_refl_exp)
         ]
     )
@@ -472,9 +472,9 @@ def test_angle_cache(self, file_handler):
             ('VIS', 2250, area_vis_exp),
             ('WV', 4500, area_ir_wv_exp),
             ('IR', 4500, area_ir_wv_exp),
-            ('quality_pixel_bitmask', 4500, area_vis_exp),
-            ('solar_zenith_angle_vis', 4500, area_vis_exp),
-            ('solar_zenith_angle_ir_wv', 2250, area_ir_wv_exp)
+            ('quality_pixel_bitmask', 2250, area_vis_exp),
+            ('solar_zenith_angle', 2250, area_vis_exp),
+            ('solar_zenith_angle', 4500, area_ir_wv_exp)
         ]
     )
     def test_get_area_definition(self, file_handler, name, resolution,
@@ -486,6 +486,8 @@ def test_get_area_definition(self, file_handler, name, resolution,
         a_exp, b_exp = proj4_radius_parameters(area_exp.proj_dict)
         assert a == a_exp
         assert b == b_exp
+        assert area.width == area_exp.width
+        assert area.height == area_exp.height
         for key in ['h', 'lon_0', 'proj', 'units']:
             assert area.proj_dict[key] == area_exp.proj_dict[key]
         np.testing.assert_allclose(area.area_extent, area_exp.area_extent)

From 81d6beed6804f67bee5a1e5c5d70c19852ddff3b Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 20 Nov 2020 15:14:38 +0000
Subject: [PATCH 27/47] Fix docstring

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index fc4e127e32..609eb8a54b 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -114,11 +114,11 @@
 
     from satpy import DataQuery
 
-    q_vis = DataQuery(
+    query_vis = DataQuery(
         name='solar_zenith_angle',
         resolution=2250
     )
-    q_ir = DataQuery(
+    query_ir = DataQuery(
         name='solar_zenith_angle',
         resolution=4500
     )

From 4f957f0c13a02157c8109952fdf048ec1d95fe09 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 20 Nov 2020 15:16:29 +0000
Subject: [PATCH 28/47] Simplify angle interpolation

---
 satpy/readers/mviri_l1b_fiduceo_nc.py                 | 8 ++++----
 satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py | 8 ++++----
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 609eb8a54b..80bef44af9 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -528,20 +528,20 @@ def _interp_angles(self, angles, dataset_id):
             target_y = self.nc.coords['y_ir_wv']
         return self._interp_angles_cached(
             angles=angles,
-            nc_key=self._get_nc_key(dataset_id['name']),
+            name=dataset_id['name'],
             target_x=target_x,
             target_y=target_y
         )
 
     @interp_cache(
-        keys=('nc_key', 'target_x', 'target_y'),
+        keys=('name', 'target_x', 'target_y'),
         hash_funcs={
-            'nc_key': lambda nc_key: nc_key,
+            'name': lambda name: name,
             'target_x': lambda x: x.size,
             'target_y': lambda y: y.size
         }
     )
-    def _interp_angles_cached(self, angles, nc_key, target_x, target_y):
+    def _interp_angles_cached(self, angles, name, target_x, target_y):
         """Interpolate angles to the given resolution."""
         return self._interp_tiepoints(angles,
                                       target_x,
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index c75e25261e..3bb1d48544 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -431,7 +431,7 @@ def test_time_cache(self, file_handler):
 
     def test_angle_cache(self, file_handler):
         """Test caching of angle datasets."""
-        nc_key = 'mykey'
+        name = 'my_angles'
         x1 = y1 = xr.DataArray([1, 2, 3, 4])
         sza_coarse = xr.DataArray(
             [[45, 90],
@@ -440,7 +440,7 @@ def test_angle_cache(self, file_handler):
         )
         sza1 = file_handler._interp_angles_cached(
             angles=sza_coarse,
-            nc_key=nc_key,
+            name=name,
             target_x=x1,
             target_y=y1
         )
@@ -449,7 +449,7 @@ def test_angle_cache(self, file_handler):
         # should not interpolate them again.
         sza2 = file_handler._interp_angles_cached(
             angles=sza_coarse - 10,
-            nc_key=nc_key,
+            name=name,
             target_x=x1,
             target_y=y1
         )
@@ -459,7 +459,7 @@ def test_angle_cache(self, file_handler):
         x2 = y2 = xr.DataArray([1, 2])
         sza3 = file_handler._interp_angles_cached(
             angles=sza_coarse,
-            nc_key=nc_key,
+            name=name,
             target_x=x2,
             target_y=y2
         )

From cdeee73d903bc06521d9b12c8581d8df189ddb13 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Mon, 23 Nov 2020 13:57:55 +0000
Subject: [PATCH 29/47] Fix VIS masking

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 80bef44af9..c5bb15a9e0 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -421,7 +421,13 @@ def _mask_vis(self, ds):
         the VIS channel.
         """
         mask = self.nc['quality_pixel_bitmask']
-        return ds.where(np.logical_or(mask == 0, mask == 2),
+        if 'y' not in mask.coords:
+            # For some reason xarray doesn't assign coordinates to the
+            # masks. Maybe because of the 'coordinates' attribute which
+            # points to non-existent latitude and longitude.
+            mask = mask.assign_coords({'y': ds.coords['y'],
+                                       'x': ds.coords['x']})
+        return ds.where(da.logical_or(mask == 0, mask == 2),
                         np.float32(np.nan))
 
     def _get_acq_time(self, ds):

From ffabcc318e5af0ca0c076a343a4d90cd1da31934 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Mon, 23 Nov 2020 14:03:44 +0000
Subject: [PATCH 30/47] Remove ancillary_variables attribute

Avoids misleading warnings from satpy.
---
 satpy/readers/mviri_l1b_fiduceo_nc.py                 | 1 +
 satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py | 9 ++++++---
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index c5bb15a9e0..9825d065e2 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -349,6 +349,7 @@ def get_area_def(self, dataset_id):
 
     def calibrate(self, ds, channel, calibration):
         """Calibrate the given dataset."""
+        ds.attrs.pop('ancillary_variables', None)  # to avoid satpy warnings
         if channel == 'VIS':
             return self._calibrate_vis(ds, calibration)
         elif channel in ['WV', 'IR']:
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 3bb1d48544..c4c9bfd679 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -338,6 +338,8 @@ def fake_dataset():
         },
         attrs={'foo': 'bar'}
     )
+    ds['count_ir'].attrs['ancillary_variables'] = 'a_ir b_ir'
+    ds['count_wv'].attrs['ancillary_variables'] = 'a_wv b_wv'
     return ds
 
 
@@ -494,15 +496,16 @@ def test_get_area_definition(self, file_handler, name, resolution,
 
     def test_calib_exceptions(self, file_handler):
         """Test calibration exceptions."""
+        ds = xr.Dataset()
         with pytest.raises(KeyError):
-            file_handler.calibrate(None, 'solar_zenith_angle', None)
+            file_handler.calibrate(ds, 'solar_zenith_angle', None)
 
         calib = mock.MagicMock()
         calib.configure_mock(name='invalid_calib')
         with pytest.raises(KeyError):
-            file_handler.calibrate(None, 'VIS', calib)
+            file_handler.calibrate(ds, 'VIS', calib)
         with pytest.raises(KeyError):
-            file_handler.calibrate(None, 'IR', calib)
+            file_handler.calibrate(ds, 'IR', calib)
 
 
 @pytest.fixture

From 15d465e0904f18acc041bfb320fc642f0fe096a4 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Mon, 23 Nov 2020 14:50:33 +0000
Subject: [PATCH 31/47] Fix coordinates of high resolution datasets

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 9825d065e2..6c2b3ba37a 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -298,8 +298,13 @@ def _read_dataset(self, name):
         ds = self.nc[nc_key]
         if 'y_ir_wv' in ds.dims:
             ds = ds.rename({'y_ir_wv': 'y', 'x_ir_wv': 'x'})
-        if 'y_tie' in ds.dims:
+        elif 'y_tie' in ds.dims:
             ds = ds.rename({'y_tie': 'y', 'x_tie': 'x'})
+        elif 'y' in ds.dims:
+            # For some reason xarray doesn't assign coordinates to all
+            # high resolution data variables.
+            ds = ds.assign_coords({'y': self.nc.coords['y'],
+                                   'x': self.nc.coords['x']})
         return ds
 
     def _update_attrs(self, ds, info):
@@ -421,13 +426,7 @@ def _mask_vis(self, ds):
         other flag set). According to [PUG] that bitmask is only applicable to
         the VIS channel.
         """
-        mask = self.nc['quality_pixel_bitmask']
-        if 'y' not in mask.coords:
-            # For some reason xarray doesn't assign coordinates to the
-            # masks. Maybe because of the 'coordinates' attribute which
-            # points to non-existent latitude and longitude.
-            mask = mask.assign_coords({'y': ds.coords['y'],
-                                       'x': ds.coords['x']})
+        mask = self._read_dataset('quality_pixel_bitmask')
         return ds.where(da.logical_or(mask == 0, mask == 2),
                         np.float32(np.nan))
 

From 20d6c187ad876ac765806cbb2faca6c71062bceb Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Mon, 23 Nov 2020 14:54:32 +0000
Subject: [PATCH 32/47] Make coordinate fix more robust

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 6c2b3ba37a..c88aa2458a 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -300,7 +300,7 @@ def _read_dataset(self, name):
             ds = ds.rename({'y_ir_wv': 'y', 'x_ir_wv': 'x'})
         elif 'y_tie' in ds.dims:
             ds = ds.rename({'y_tie': 'y', 'x_tie': 'x'})
-        elif 'y' in ds.dims:
+        elif 'y' in ds.dims and 'y' not in ds.coords:
             # For some reason xarray doesn't assign coordinates to all
             # high resolution data variables.
             ds = ds.assign_coords({'y': self.nc.coords['y'],

From da3df9be661c8a33cdfa7363d301634b7f7702d7 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Mon, 23 Nov 2020 16:34:35 +0000
Subject: [PATCH 33/47] Simplify masking

Furthermore, issue a warning if all pixels are flagged as
"use with caution".
---
 satpy/readers/mviri_l1b_fiduceo_nc.py         | 55 ++++++++++++++-----
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py |  4 +-
 2 files changed, 43 insertions(+), 16 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index c88aa2458a..e9c30dc81a 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -52,13 +52,10 @@
     from satpy import Scene
 
     scn = Scene(filenames=['FIDUCEO_FCDR_L15_MVIRI_MET7-57.0...'],
-                reader='mviri_l1b_fiduceo_nc',
-                reader_kwargs={'mask_bad_quality': True)
+                reader='mviri_l1b_fiduceo_nc')
     scn.load(['VIS', 'WV', 'IR'])
 
-In the above example pixels considered bad quality are masked, see
-:class:`FiduceoMviriBase` for a keyword argument description. Global
-netCDF attributes are available in the ``raw_metadata`` attribute of
+Global netCDF attributes are available in the ``raw_metadata`` attribute of
 each loaded dataset.
 
 
@@ -98,6 +95,23 @@
 uncertainties can be used to filter these cases before calculating reflectances.
 
 
+VIS Channel Quality Flags
+-------------------------
+Quality flags are available for the VIS channel only. A simple approach for
+masking bad quality pixels is to set the ``mask_bad_quality`` keyword argument
+to ``True``:
+
+.. code-block:: python
+    scn = Scene(filenames=['FIDUCEO_FCDR_L15_MVIRI_MET7-57.0...'],
+                reader='mviri_l1b_fiduceo_nc',
+                reader_kwargs={'mask_bad_quality': True})
+
+See :class:`FiduceoMviriBase` for an argument description. In some situations
+however the entire image can be flagged (look out for warnings). In that case
+check out the ``quality_pixel_bitmask`` and ``data_quality_bitmask`` datasets
+to find out why.
+
+
 Angles
 ------
 The FIDUCEO MVIRI FCDR provides satellite and solar angles on a coarse tiepoint
@@ -141,6 +155,7 @@
 
 import abc
 import functools
+import warnings
 
 import dask.array as da
 import numpy as np
@@ -252,9 +267,9 @@ def __init__(self, filename, filename_info, filetype_info,
 
         Args:
              mask_bad_quality: Mask VIS pixels with bad quality, that means
-                 everything else than "ok" or "use with caution". If you need
-                 more control, use the ``quality_pixel_bitmask`` and
-                 ``data_quality_bitmask`` datasets.
+                 any quality flag except "ok". If you need more control, use
+                 the ``quality_pixel_bitmask`` and ``data_quality_bitmask``
+                 datasets.
         """
         super(FiduceoMviriBase, self).__init__(
             filename, filename_info, filetype_info)
@@ -278,8 +293,11 @@ def get_dataset(self, dataset_id, dataset_info):
         if dataset_id['name'] in CHANNELS:
             ds = self.calibrate(ds, channel=name,
                                 calibration=dataset_id['calibration'])
-            if self.mask_bad_quality and name == 'VIS':
-                ds = self._mask_vis(ds)
+            if name == 'VIS':
+                if self.mask_bad_quality:
+                    ds = self._mask_vis(ds)
+                else:
+                    self._check_vis_quality(ds)
             ds['acq_time'] = ('y', self._get_acq_time(ds))
         elif dataset_id['name'] in OTHER_REFLECTANCES:
             ds = ds * 100  # conversion to percent
@@ -418,16 +436,25 @@ def _ir_wv_radiance_to_brightness_temperature(self, rad, channel):
         bt = b / (np.log(rad) - a)
         return bt.where(bt > 0, np.float32(np.nan))
 
+    def _check_vis_quality(self, ds):
+        """Check VIS channel quality and issue a warning if it's bad."""
+        mask = self._read_dataset('quality_pixel_bitmask')
+        use_with_caution = da.bitwise_and(mask, 2)
+        if use_with_caution.all():
+            warnings.warn(
+                'All pixels of the VIS channel are flagged as "use with '
+                'caution". Use datasets "quality_pixel_bitmask" and '
+                '"data_quality_bitmask" to find out why.'
+            )
+
     def _mask_vis(self, ds):
         """Mask VIS pixels with bad quality.
 
         Pixels are considered bad quality if the "quality_pixel_bitmask" is
-        everything else than 0 (no flag set) or 2 ("use_with_caution" and no
-        other flag set). According to [PUG] that bitmask is only applicable to
-        the VIS channel.
+        everything else than 0 (no flag set).
         """
         mask = self._read_dataset('quality_pixel_bitmask')
-        return ds.where(da.logical_or(mask == 0, mask == 2),
+        return ds.where(mask == 0,
                         np.float32(np.nan))
 
     def _get_acq_time(self, ds):
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index c4c9bfd679..5b2ed25ae2 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -208,7 +208,7 @@
 )
 quality_pixel_bitmask_exp = xr.DataArray(
     np.array(
-        [[2, 0, 0, 0],
+        [[0, 0, 0, 0],
          [0, 0, 0, 0],
          [0, 0, 1, 0],
          [0, 0, 0, 0]],
@@ -287,7 +287,7 @@ def fake_dataset():
     )
     mask = da.from_array(
         np.array(
-            [[2, 0, 0, 0],  # 2 = "use with caution"
+            [[0, 0, 0, 0],
              [0, 0, 0, 0],
              [0, 0, 1, 0],  # 1 = "invalid"
              [0, 0, 0, 0]],

From e40dd3bed99c0097dedf9661fb9658494cfa3dcc Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Mon, 23 Nov 2020 16:54:12 +0000
Subject: [PATCH 34/47] Fix documentation

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index e9c30dc81a..a5c2705a48 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -102,6 +102,7 @@
 to ``True``:
 
 .. code-block:: python
+
     scn = Scene(filenames=['FIDUCEO_FCDR_L15_MVIRI_MET7-57.0...'],
                 reader='mviri_l1b_fiduceo_nc',
                 reader_kwargs={'mask_bad_quality': True})

From 8633801a7a33fa3a19a32167e4df9a5649436c0f Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Mon, 23 Nov 2020 17:06:15 +0000
Subject: [PATCH 35/47] Add test for warning

---
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py     | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 5b2ed25ae2..2d59f9629d 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -347,6 +347,10 @@ def fake_dataset():
                         FiduceoMviriFullFcdrFileHandler])
 def file_handler(fake_dataset, request):
     """Create mocked file handler."""
+    marker = request.node.get_closest_marker("file_handler_data")
+    mask_bad_quality = True
+    if marker:
+        mask_bad_quality = marker.kwargs['mask_bad_quality']
     fh_class = request.param
     with mock.patch('satpy.readers.mviri_l1b_fiduceo_nc.xr.open_dataset') as open_dataset:
         open_dataset.return_value = fake_dataset
@@ -356,7 +360,7 @@ def file_handler(fake_dataset, request):
                            'sensor': 'MVIRI',
                            'projection_longitude': '57.0'},
             filetype_info={'foo': 'bar'},
-            mask_bad_quality=True
+            mask_bad_quality=mask_bad_quality
         )
 
 
@@ -507,6 +511,15 @@ def test_calib_exceptions(self, file_handler):
         with pytest.raises(KeyError):
             file_handler.calibrate(ds, 'IR', calib)
 
+    @pytest.mark.file_handler_data(mask_bad_quality=False)
+    def test_bad_quality_warning(self, file_handler):
+        """Test warning about bad VIS quality."""
+        file_handler.nc['quality_pixel_bitmask'] = 2
+        vis = make_dataid(name='VIS', resolution=2250,
+                          calibration='reflectance')
+        with pytest.warns(UserWarning):
+            file_handler.get_dataset(vis, {})
+
 
 @pytest.fixture
 def reader():

From 22aa270811832d687610abd629fc788e1ee9628f Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Tue, 24 Nov 2020 15:42:06 +0000
Subject: [PATCH 36/47] Use lru_cache instead of custom cache

---
 satpy/readers/mviri_l1b_fiduceo_nc.py         | 132 +++++----------
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py | 155 +++++++++---------
 2 files changed, 113 insertions(+), 174 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index a5c2705a48..aa6046ce24 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -155,7 +155,7 @@
 """
 
 import abc
-import functools
+from functools import lru_cache
 import warnings
 
 import dask.array as da
@@ -163,7 +163,6 @@
 import xarray as xr
 
 from satpy import CHUNK_SIZE
-from satpy.dataset.dataid import DataQuery
 from satpy.readers._geos_area import (ang2fac, get_area_definition,
                                       get_area_extent)
 from satpy.readers.file_handlers import BaseFileHandler
@@ -190,49 +189,6 @@
 HIGH_RESOL = 2250
 
 
-class InterpCache:
-    """Interpolation cache."""
-
-    def __init__(self, func, keys, hash_funcs):
-        """Create the cache.
-
-        Args:
-            func:
-                Interpolation function to be cached.
-            keys:
-                Function arguments serving as cache key.
-            hash_funcs (dict):
-                For each key provides a method that extracts a hashable
-                value from the given keyword argument. For example, use
-                image shape to maintain separate caches for low
-                resolution (WV/IR) and high resolution (VIS) channels.
-        """
-        self.func = func
-        self.keys = keys
-        self.hash_funcs = hash_funcs
-        self.cache = {}
-
-    def __call__(self, *args, **kwargs):
-        """Call the interpolation function."""
-        key = tuple(
-            [self.hash_funcs[key](kwargs[key]) for key in self.keys]
-        )
-        if key not in self.cache:
-            self.cache[key] = self.func(*args, **kwargs)
-        return self.cache[key]
-
-    def __get__(self, obj, objtype):
-        """To support instance methods."""
-        return functools.partial(self.__call__, obj)
-
-
-def interp_cache(keys, hash_funcs):
-    """Interpolation cache."""
-    def wrapper(func):
-        return InterpCache(func, keys, hash_funcs)
-    return wrapper
-
-
 class FiduceoMviriBase(BaseFileHandler):
     """Baseclass for FIDUCEO MVIRI file handlers."""
     nc_keys = {
@@ -290,20 +246,22 @@ def __init__(self, filename, filename_info, filetype_info,
     def get_dataset(self, dataset_id, dataset_info):
         """Get the dataset."""
         name = dataset_id['name']
-        ds = self._read_dataset(name)
-        if dataset_id['name'] in CHANNELS:
-            ds = self.calibrate(ds, channel=name,
-                                calibration=dataset_id['calibration'])
-            if name == 'VIS':
-                if self.mask_bad_quality:
-                    ds = self._mask_vis(ds)
-                else:
-                    self._check_vis_quality(ds)
-            ds['acq_time'] = ('y', self._get_acq_time(ds))
-        elif dataset_id['name'] in OTHER_REFLECTANCES:
-            ds = ds * 100  # conversion to percent
-        elif dataset_id['name'] in ANGLES:
-            ds = self._interp_angles(ds, dataset_id)
+        resolution = dataset_id['resolution']
+        if name in ANGLES:
+            ds = self._get_angles(name, resolution)
+        else:
+            ds = self._read_dataset(name)
+            if name in CHANNELS:
+                ds = self.calibrate(ds, channel=name,
+                                    calibration=dataset_id['calibration'])
+                if name == 'VIS':
+                    if self.mask_bad_quality:
+                        ds = self._mask_vis(ds)
+                    else:
+                        self._check_vis_quality(ds)
+                ds['acq_time'] = ('y', self._get_acq_time(resolution))
+            elif name in OTHER_REFLECTANCES:
+                ds = ds * 100  # conversion to percent
         self._update_attrs(ds, dataset_info)
         return ds
 
@@ -336,7 +294,7 @@ def _update_attrs(self, ds, info):
 
     def get_area_def(self, dataset_id):
         """Get area definition of the given dataset."""
-        if self._is_high_resol(dataset_id):
+        if self._is_high_resol(dataset_id['resolution']):
             im_size = self.nc.coords['y'].size
             area_name = 'geos_mviri_vis'
         else:
@@ -458,8 +416,9 @@ def _mask_vis(self, ds):
         return ds.where(mask == 0,
                         np.float32(np.nan))
 
-    def _get_acq_time(self, ds):
-        """Get scanline acquisition time for the given dataset.
+    @lru_cache(maxsize=3)  # Three channels
+    def _get_acq_time(self, resolution):
+        """Get scanline acquisition time for the given resolution.
 
         Note that the acquisition time does not increase monotonically
         with the scanline number due to the scan pattern and rectification.
@@ -469,14 +428,14 @@ def _get_acq_time(self, ds):
             time2d = self.nc['time_ir_wv']
         except KeyError:
             time2d = self.nc['time']
-        return self._get_acq_time_cached(time2d, target_y=ds.coords['y'])
+        if self._is_high_resol(resolution):
+            target_y = self.nc.coords['x']
+        else:
+            target_y = self.nc.coords['x_ir_wv']
+        return self._interp_acq_time(time2d, target_y=target_y.values)
 
-    @interp_cache(
-        keys=('target_y',),
-        hash_funcs={'target_y': lambda y: y.size}
-    )
-    def _get_acq_time_cached(self, time2d, target_y):
-        """Get scanline acquisition time for the given image coordinates.
+    def _interp_acq_time(self, time2d, target_y):
+        """Interpolate scanline acquisition time to the given coordinates.
 
         The files provide timestamps per pixel for the low resolution
         channels (IR/WV) only.
@@ -548,39 +507,26 @@ def _get_ssp(self, coord):
             sub_lonlat = np.nan
         return sub_lonlat
 
-    def _interp_angles(self, angles, dataset_id):
+    @lru_cache(maxsize=8)  # 4 angle datasets with two resolutions each
+    def _get_angles(self, name, resolution):
         """Get angle dataset.
 
         Files provide angles (solar/satellite zenith & azimuth) at a coarser
         resolution. Interpolate them to the desired resolution.
         """
-        if self._is_high_resol(dataset_id):
+        angles = self._read_dataset(name)
+        if self._is_high_resol(resolution):
             target_x = self.nc.coords['x']
             target_y = self.nc.coords['y']
         else:
             target_x = self.nc.coords['x_ir_wv']
             target_y = self.nc.coords['y_ir_wv']
-        return self._interp_angles_cached(
-            angles=angles,
-            name=dataset_id['name'],
+        return self._interp_tiepoints(
+            angles,
             target_x=target_x,
             target_y=target_y
         )
 
-    @interp_cache(
-        keys=('name', 'target_x', 'target_y'),
-        hash_funcs={
-            'name': lambda name: name,
-            'target_x': lambda x: x.size,
-            'target_y': lambda y: y.size
-        }
-    )
-    def _interp_angles_cached(self, angles, name, target_x, target_y):
-        """Interpolate angles to the given resolution."""
-        return self._interp_tiepoints(angles,
-                                      target_x,
-                                      target_y)
-
     def _interp_tiepoints(self, ds, target_x, target_y):
         """Interpolate dataset between tiepoints.
 
@@ -605,8 +551,8 @@ def _interp_tiepoints(self, ds, target_x, target_y):
 
         return ds.interp(x=target_x.values, y=target_y.values)
 
-    def _is_high_resol(self, dataset_id):
-        return dataset_id['resolution'] == HIGH_RESOL
+    def _is_high_resol(self, resolution):
+        return resolution == HIGH_RESOL
 
 
 class FiduceoMviriEasyFcdrFileHandler(FiduceoMviriBase):
@@ -677,11 +623,7 @@ def _vis_radiance_to_reflectance(self, rad):
 
         Reference: [PUG], equation (6).
         """
-        sza = self.get_dataset(
-            DataQuery(name='solar_zenith_angle',
-                      resolution=HIGH_RESOL),
-            dataset_info={}
-        )
+        sza = self._get_angles('solar_zenith_angle', HIGH_RESOL)
         sza = sza.where(da.fabs(sza) < 90,
                         np.float32(np.nan))  # direct illumination only
         cos_sza = np.cos(np.deg2rad(sza))
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 2d59f9629d..229b37a2dd 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -364,6 +364,18 @@ def file_handler(fake_dataset, request):
         )
 
 
+@pytest.fixture
+def reader():
+    """Return MVIRI FIDUCEO FCDR reader."""
+    from satpy.config import config_search_paths
+    from satpy.readers import load_reader
+
+    reader_configs = config_search_paths(
+        os.path.join("readers", "mviri_l1b_fiduceo_nc.yaml"))
+    reader = load_reader(reader_configs)
+    return reader
+
+
 class TestFiduceoMviriFileHandlers:
     """Unit tests for FIDUCEO MVIRI file handlers."""
 
@@ -412,65 +424,63 @@ def test_get_dataset(self, file_handler, name, calibration, resolution,
             assert ds.dtype == expected.dtype
             assert ds.attrs == expected.attrs
 
-    def test_time_cache(self, file_handler):
+    @mock.patch(
+        'satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriBase._interp_acq_time'
+    )
+    def test_time_cache(self, interp_acq_time, file_handler):
         """Test caching of acquisition times."""
-        time2d = xr.DataArray(
-            np.array([[1, 2],
-                      [3, 4]],
-                     dtype='datetime64[h]'),
-            dims=('y_ir_wv', 'x_ir_wv')
-        )
-        y1 = xr.DataArray([1, 2, 3, 4])
-        t1 = file_handler._get_acq_time_cached(time2d, target_y=y1)
-
-        # Change 2d timestamps. If the cache works correctly, the second call
-        # should not average/interpolate them again.
-        t2 = file_handler._get_acq_time_cached(
-            time2d + np.timedelta64(1, 'h'), target_y=y1)
-        xr.testing.assert_equal(t2, t1)
-
-        # With new target coordinates we shouldn't hit the cache
-        y2 = xr.DataArray([1, 2])
-        t3 = file_handler._get_acq_time_cached(target_y=y2)
-        with pytest.raises(AssertionError):
-            xr.testing.assert_equal(t3, t1)
-
-    def test_angle_cache(self, file_handler):
-        """Test caching of angle datasets."""
-        name = 'my_angles'
-        x1 = y1 = xr.DataArray([1, 2, 3, 4])
-        sza_coarse = xr.DataArray(
-            [[45, 90],
-             [0, 45]],
-            dims=('y', 'x')
-        )
-        sza1 = file_handler._interp_angles_cached(
-            angles=sza_coarse,
-            name=name,
-            target_x=x1,
-            target_y=y1
+        dataset_id = make_dataid(
+            name='VIS',
+            resolution=2250,
+            calibration='reflectance'
         )
+        info = {}
+        interp_acq_time.return_value = xr.DataArray([1, 2, 3, 4], dims='y')
 
-        # Change coarse angles. If the cache works correctly, the second call
-        # should not interpolate them again.
-        sza2 = file_handler._interp_angles_cached(
-            angles=sza_coarse - 10,
-            name=name,
-            target_x=x1,
-            target_y=y1
-        )
-        xr.testing.assert_equal(sza2, sza1)
-
-        # With new target coordinates we shouldn't hit the cache
-        x2 = y2 = xr.DataArray([1, 2])
-        sza3 = file_handler._interp_angles_cached(
-            angles=sza_coarse,
-            name=name,
-            target_x=x2,
-            target_y=y2
+        # Cache init
+        file_handler.get_dataset(dataset_id, info)
+        interp_acq_time.assert_called()
+
+        # Cache hit
+        interp_acq_time.reset_mock()
+        file_handler.get_dataset(dataset_id, info)
+        interp_acq_time.assert_not_called()
+
+        # Cache miss
+        interp_acq_time.return_value = xr.DataArray([1, 2], dims='y')
+        another_id = make_dataid(
+            name='IR',
+            resolution=4500,
+            calibration='brightness_temperature'
         )
-        with pytest.raises(AssertionError):
-            xr.testing.assert_equal(sza3, sza1)
+        interp_acq_time.reset_mock()
+        file_handler.get_dataset(another_id, info)
+        interp_acq_time.assert_called()
+
+    @mock.patch(
+        'satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriBase._interp_tiepoints'
+    )
+    def test_angle_cache(self, interp_tiepoints, file_handler):
+        """Test caching of angle datasets."""
+        dataset_id = make_dataid(name='solar_zenith_angle',
+                                 resolution=2250)
+        info = {}
+
+        # Cache init
+        file_handler.get_dataset(dataset_id, info)
+        interp_tiepoints.assert_called()
+
+        # Cache hit
+        interp_tiepoints.reset_mock()
+        file_handler.get_dataset(dataset_id, info)
+        interp_tiepoints.assert_not_called()
+
+        # Cache miss
+        another_id = make_dataid(name='solar_zenith_angle',
+                                 resolution=4500)
+        interp_tiepoints.reset_mock()
+        file_handler.get_dataset(another_id, info)
+        interp_tiepoints.assert_called()
 
     @pytest.mark.parametrize(
         ('name', 'resolution', 'area_exp'),
@@ -520,28 +530,15 @@ def test_bad_quality_warning(self, file_handler):
         with pytest.warns(UserWarning):
             file_handler.get_dataset(vis, {})
 
+    def test_file_pattern(self, reader):
+        """Test file pattern matching."""
+        filenames = [
+            "FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_FULL_v2.6_fv3.1.nc",
+            "FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc",
+            "FIDUCEO_FCDR_L15_MVIRI_MET7-00.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc",
+            "abcde",
+        ]
 
-@pytest.fixture
-def reader():
-    """Return MVIRI FIDUCEO FCDR reader."""
-    from satpy.config import config_search_paths
-    from satpy.readers import load_reader
-
-    reader_configs = config_search_paths(
-        os.path.join("readers", "mviri_l1b_fiduceo_nc.yaml"))
-    reader = load_reader(reader_configs)
-    return reader
-
-
-def test_file_pattern(reader):
-    """Test file pattern matching."""
-    filenames = [
-        "FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_FULL_v2.6_fv3.1.nc",
-        "FIDUCEO_FCDR_L15_MVIRI_MET7-57.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc",
-        "FIDUCEO_FCDR_L15_MVIRI_MET7-00.0_201701201000_201701201030_EASY_v2.6_fv3.1.nc",
-        "abcde",
-    ]
-
-    files = reader.select_files_from_pathnames(filenames)
-    # only 3 out of 4 above should match
-    assert len(files) == 3
+        files = reader.select_files_from_pathnames(filenames)
+        # only 3 out of 4 above should match
+        assert len(files) == 3

From 2c74ec1f44027170a47dbd3e8f5b62c928267965 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Tue, 24 Nov 2020 15:48:01 +0000
Subject: [PATCH 37/47] Clarify method name

---
 satpy/readers/_geos_area.py           |  6 +++---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 11 ++++++++---
 2 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/satpy/readers/_geos_area.py b/satpy/readers/_geos_area.py
index d6632ab8f6..51b39dd3c8 100644
--- a/satpy/readers/_geos_area.py
+++ b/satpy/readers/_geos_area.py
@@ -147,7 +147,7 @@ def get_area_definition(pdict, a_ext):
     return a_def
 
 
-def ang2fac(ang):
+def sampling_to_lfac_cfac(sampling):
     """Convert angular sampling to line/column scaling factor (aka LFAC/CFAC).
 
     Reference: `MSG Ground Segment LRIT HRIT Mission Specific Implementation`_,
@@ -159,10 +159,10 @@ def ang2fac(ang):
 &RevisionSelectionMethod=LatestReleased&Rendition=Web
 
     Args:
-        ang: float
+        sampling: float
             Angular sampling (rad)
 
     Returns:
         Line/column scaling factor (deg-1)
     """
-    return 2.0 ** 16 / np.rad2deg(ang)
+    return 2.0 ** 16 / np.rad2deg(sampling)
diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index aa6046ce24..1b00d93864 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -163,8 +163,11 @@
 import xarray as xr
 
 from satpy import CHUNK_SIZE
-from satpy.readers._geos_area import (ang2fac, get_area_definition,
-                                      get_area_extent)
+from satpy.readers._geos_area import (
+    sampling_to_lfac_cfac,
+    get_area_definition,
+    get_area_extent
+)
 from satpy.readers.file_handlers import BaseFileHandler
 
 EQUATOR_RADIUS = 6378140.0
@@ -306,7 +309,9 @@ def get_area_def(self, dataset_id):
         # [Handbook] and in
         # https://github.com/FIDUCEO/FCDR_MVIRI/blob/master/lib/nrCrunch/cruncher.f
         loff = coff = im_size / 2 + 0.5
-        lfac = cfac = ang2fac(np.deg2rad(MVIRI_FIELD_OF_VIEW) / im_size)
+        lfac = cfac = sampling_to_lfac_cfac(
+            np.deg2rad(MVIRI_FIELD_OF_VIEW) / im_size
+        )
 
         pdict = {
             'ssp_lon': self.projection_longitude,

From 35877dc003d69ad0285184ccf1f0467c62c15b2e Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Tue, 24 Nov 2020 15:52:41 +0000
Subject: [PATCH 38/47] Rename fixtures to avoid arguments being masked

---
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py   | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 229b37a2dd..1eba781de4 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -272,8 +272,8 @@
 )
 
 
-@pytest.fixture()
-def fake_dataset():
+@pytest.fixture(name='fake_dataset')
+def fixture_fake_dataset():
     """Create fake dataset."""
     count_ir = da.linspace(0, 255, 4, dtype=np.uint8).reshape(2, 2)
     count_wv = da.linspace(0, 255, 4, dtype=np.uint8).reshape(2, 2)
@@ -343,9 +343,12 @@ def fake_dataset():
     return ds
 
 
-@pytest.fixture(params=[FiduceoMviriEasyFcdrFileHandler,
-                        FiduceoMviriFullFcdrFileHandler])
-def file_handler(fake_dataset, request):
+@pytest.fixture(
+    name='file_handler',
+    params=[FiduceoMviriEasyFcdrFileHandler,
+            FiduceoMviriFullFcdrFileHandler]
+)
+def fixture_file_handler(fake_dataset, request):
     """Create mocked file handler."""
     marker = request.node.get_closest_marker("file_handler_data")
     mask_bad_quality = True
@@ -364,8 +367,8 @@ def file_handler(fake_dataset, request):
         )
 
 
-@pytest.fixture
-def reader():
+@pytest.fixture(name='reader')
+def fixture_reader():
     """Return MVIRI FIDUCEO FCDR reader."""
     from satpy.config import config_search_paths
     from satpy.readers import load_reader

From eb2ebf51e03671dd52f2e19085db70682b710e15 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Tue, 24 Nov 2020 16:12:54 +0000
Subject: [PATCH 39/47] Simplify get_dataset

---
 satpy/readers/mviri_l1b_fiduceo_nc.py         | 134 ++++++++++--------
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py |  19 ++-
 2 files changed, 86 insertions(+), 67 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 1b00d93864..1b4d2e1151 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -252,49 +252,13 @@ def get_dataset(self, dataset_id, dataset_info):
         resolution = dataset_id['resolution']
         if name in ANGLES:
             ds = self._get_angles(name, resolution)
+        elif name in CHANNELS:
+            ds = self._get_channel(name, resolution, dataset_id['calibration'])
         else:
-            ds = self._read_dataset(name)
-            if name in CHANNELS:
-                ds = self.calibrate(ds, channel=name,
-                                    calibration=dataset_id['calibration'])
-                if name == 'VIS':
-                    if self.mask_bad_quality:
-                        ds = self._mask_vis(ds)
-                    else:
-                        self._check_vis_quality(ds)
-                ds['acq_time'] = ('y', self._get_acq_time(resolution))
-            elif name in OTHER_REFLECTANCES:
-                ds = ds * 100  # conversion to percent
+            ds = self._get_other_dataset(name)
         self._update_attrs(ds, dataset_info)
         return ds
 
-    def _get_nc_key(self, ds_name):
-        """Get netCDF variable name for the given dataset."""
-        return self.nc_keys.get(ds_name, ds_name)
-
-    def _read_dataset(self, name):
-        """Read a dataset from the file."""
-        nc_key = self._get_nc_key(name)
-        ds = self.nc[nc_key]
-        if 'y_ir_wv' in ds.dims:
-            ds = ds.rename({'y_ir_wv': 'y', 'x_ir_wv': 'x'})
-        elif 'y_tie' in ds.dims:
-            ds = ds.rename({'y_tie': 'y', 'x_tie': 'x'})
-        elif 'y' in ds.dims and 'y' not in ds.coords:
-            # For some reason xarray doesn't assign coordinates to all
-            # high resolution data variables.
-            ds = ds.assign_coords({'y': self.nc.coords['y'],
-                                   'x': self.nc.coords['x']})
-        return ds
-
-    def _update_attrs(self, ds, info):
-        """Update dataset attributes."""
-        ds.attrs.update(info)
-        ds.attrs.update({'platform': self.filename_info['platform'],
-                         'sensor': self.filename_info['sensor']})
-        ds.attrs['raw_metadata'] = self.nc.attrs
-        ds.attrs['orbital_parameters'] = self._get_orbital_parameters()
-
     def get_area_def(self, dataset_id):
         """Get area definition of the given dataset."""
         if self._is_high_resol(dataset_id['resolution']):
@@ -334,7 +298,77 @@ def get_area_def(self, dataset_id):
         area_def = get_area_definition(pdict, extent)
         return area_def
 
-    def calibrate(self, ds, channel, calibration):
+    def _get_channel(self, name, resolution, calibration):
+        """Get and calibrate channel data."""
+        ds = self._read_dataset(name)
+        ds = self._calibrate(
+            ds,
+            channel=name,
+            calibration=calibration
+        )
+        if name == 'VIS':
+            if self.mask_bad_quality:
+                ds = self._mask_vis(ds)
+            else:
+                self._check_vis_quality(ds)
+        ds['acq_time'] = ('y', self._get_acq_time(resolution))
+        return ds
+
+    @lru_cache(maxsize=8)  # 4 angle datasets with two resolutions each
+    def _get_angles(self, name, resolution):
+        """Get angle dataset.
+
+        Files provide angles (solar/satellite zenith & azimuth) at a coarser
+        resolution. Interpolate them to the desired resolution.
+        """
+        angles = self._read_dataset(name)
+        if self._is_high_resol(resolution):
+            target_x = self.nc.coords['x']
+            target_y = self.nc.coords['y']
+        else:
+            target_x = self.nc.coords['x_ir_wv']
+            target_y = self.nc.coords['y_ir_wv']
+        return self._interp_tiepoints(
+            angles,
+            target_x=target_x,
+            target_y=target_y
+        )
+
+    def _get_other_dataset(self, name):
+        """Get other datasets such as uncertainties."""
+        ds = self._read_dataset(name)
+        if name in OTHER_REFLECTANCES:
+            ds = ds * 100  # conversion to percent
+        return ds
+
+    def _get_nc_key(self, ds_name):
+        """Get netCDF variable name for the given dataset."""
+        return self.nc_keys.get(ds_name, ds_name)
+
+    def _read_dataset(self, name):
+        """Read a dataset from the file."""
+        nc_key = self._get_nc_key(name)
+        ds = self.nc[nc_key]
+        if 'y_ir_wv' in ds.dims:
+            ds = ds.rename({'y_ir_wv': 'y', 'x_ir_wv': 'x'})
+        elif 'y_tie' in ds.dims:
+            ds = ds.rename({'y_tie': 'y', 'x_tie': 'x'})
+        elif 'y' in ds.dims and 'y' not in ds.coords:
+            # For some reason xarray doesn't assign coordinates to all
+            # high resolution data variables.
+            ds = ds.assign_coords({'y': self.nc.coords['y'],
+                                   'x': self.nc.coords['x']})
+        return ds
+
+    def _update_attrs(self, ds, info):
+        """Update dataset attributes."""
+        ds.attrs.update(info)
+        ds.attrs.update({'platform': self.filename_info['platform'],
+                         'sensor': self.filename_info['sensor']})
+        ds.attrs['raw_metadata'] = self.nc.attrs
+        ds.attrs['orbital_parameters'] = self._get_orbital_parameters()
+
+    def _calibrate(self, ds, channel, calibration):
         """Calibrate the given dataset."""
         ds.attrs.pop('ancillary_variables', None)  # to avoid satpy warnings
         if channel == 'VIS':
@@ -512,26 +546,6 @@ def _get_ssp(self, coord):
             sub_lonlat = np.nan
         return sub_lonlat
 
-    @lru_cache(maxsize=8)  # 4 angle datasets with two resolutions each
-    def _get_angles(self, name, resolution):
-        """Get angle dataset.
-
-        Files provide angles (solar/satellite zenith & azimuth) at a coarser
-        resolution. Interpolate them to the desired resolution.
-        """
-        angles = self._read_dataset(name)
-        if self._is_high_resol(resolution):
-            target_x = self.nc.coords['x']
-            target_y = self.nc.coords['y']
-        else:
-            target_x = self.nc.coords['x_ir_wv']
-            target_y = self.nc.coords['y_ir_wv']
-        return self._interp_tiepoints(
-            angles,
-            target_x=target_x,
-            target_y=target_y
-        )
-
     def _interp_tiepoints(self, ds, target_x, target_y):
         """Interpolate dataset between tiepoints.
 
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 1eba781de4..09571e0ccc 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -513,16 +513,21 @@ def test_get_area_definition(self, file_handler, name, resolution,
 
     def test_calib_exceptions(self, file_handler):
         """Test calibration exceptions."""
-        ds = xr.Dataset()
         with pytest.raises(KeyError):
-            file_handler.calibrate(ds, 'solar_zenith_angle', None)
-
-        calib = mock.MagicMock()
-        calib.configure_mock(name='invalid_calib')
+            file_handler.get_dataset(
+                make_dataid(name='solar_zenith_angle', calibration='counts'),
+                {}
+            )
         with pytest.raises(KeyError):
-            file_handler.calibrate(ds, 'VIS', calib)
+            file_handler.get_dataset(
+                {'name': 'VIS', 'calibration': 'invalid'},
+                {}
+            )
         with pytest.raises(KeyError):
-            file_handler.calibrate(ds, 'IR', calib)
+            file_handler.get_dataset(
+                {'name': 'IR', 'calibration': 'invalid'},
+                {}
+            )
 
     @pytest.mark.file_handler_data(mask_bad_quality=False)
     def test_bad_quality_warning(self, file_handler):

From 6f51ca938c3361702b81f6edd945ef3a72f1bdab Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Tue, 24 Nov 2020 16:23:54 +0000
Subject: [PATCH 40/47] Fix geos_area test

---
 satpy/tests/reader_tests/test_geos_area.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/satpy/tests/reader_tests/test_geos_area.py b/satpy/tests/reader_tests/test_geos_area.py
index f4bf663962..87e0c2cf46 100644
--- a/satpy/tests/reader_tests/test_geos_area.py
+++ b/satpy/tests/reader_tests/test_geos_area.py
@@ -21,7 +21,7 @@
 from satpy.readers._geos_area import (get_xy_from_linecol,
                                       get_area_extent,
                                       get_area_definition,
-                                      ang2fac)
+                                      sampling_to_lfac_cfac)
 
 import numpy as np
 
@@ -147,8 +147,8 @@ def test_get_area_definition(self):
         self.assertEqual(b, 6356583.8)
         self.assertEqual(a_def.proj_dict['h'], 35785831)
 
-    def test_ang2fac(self):
+    def test_sampling_to_lfac_cfac(self):
         """Test conversion from angular sampling to line/column offset."""
         lfac = 13642337  # SEVIRI LFAC
         sampling = np.deg2rad(2 ** 16 / lfac)
-        np.testing.assert_allclose(ang2fac(sampling), lfac)
+        np.testing.assert_allclose(sampling_to_lfac_cfac(sampling), lfac)

From 180929bc0ebf1ce9d0bd026de00f811fcfe88c42 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 25 Nov 2020 09:50:39 +0000
Subject: [PATCH 41/47] Reduce complexity of the file handler base class

---
 satpy/readers/mviri_l1b_fiduceo_nc.py         | 412 +++++++++++-------
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py |   8 +-
 2 files changed, 248 insertions(+), 172 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 1b4d2e1151..aedbbb07d4 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -192,8 +192,193 @@
 HIGH_RESOL = 2250
 
 
+class IRWVCalibrator:
+    """Calibrate IR & WV channels."""
+
+    def __init__(self, coefs):
+        """Initialize the calibrator.
+
+        Args:
+            coefs: Calibration coefficients.
+        """
+        self.coefs = coefs
+
+    def calibrate(self, counts, calibration):
+        """Calibrate IR/WV counts to the given calibration."""
+        if calibration == 'counts':
+            return counts
+        elif calibration in ('radiance', 'brightness_temperature'):
+            rad = self._counts_to_radiance(counts)
+            if calibration == 'radiance':
+                return rad
+            bt = self._radiance_to_brightness_temperature(rad)
+            return bt
+        else:
+            raise KeyError('Invalid calibration: {}'.format(calibration.name))
+
+    def _counts_to_radiance(self, counts):
+        """Convert IR/WV counts to radiance.
+
+        Reference: [PUG], equations (4.1) and (4.2).
+        """
+        a, b = self.coefs['radiance']
+        rad = a + b * counts
+        return rad.where(rad > 0, np.float32(np.nan))
+
+    def _radiance_to_brightness_temperature(self, rad):
+        """Convert IR/WV radiance to brightness temperature.
+
+        Reference: [PUG], equations (5.1) and (5.2).
+        """
+        a, b = self.coefs['brightness_temperature']
+        bt = b / (np.log(rad) - a)
+        return bt.where(bt > 0, np.float32(np.nan))
+
+
+class VISCalibrator:
+    """Calibrate VIS channel."""
+
+    def __init__(self, coefs, solar_zenith_angle=None):
+        """Initialize the calibrator.
+
+        Args:
+            coefs:
+                Calibration coefficients.
+            solar_zenith_angle (optional):
+                Solar zenith angle. Only required for calibration to
+                reflectance.
+        """
+        self.coefs = coefs
+        self.solar_zenith_angle = solar_zenith_angle
+
+    def calibrate(self, counts, calibration):
+        """Calibrate VIS counts."""
+        if calibration == 'counts':
+            return counts
+        elif calibration in ('radiance', 'reflectance'):
+            rad = self._counts_to_radiance(counts)
+            if calibration == 'radiance':
+                return rad
+            refl = self._radiance_to_reflectance(rad)
+            refl = self.update_refl_attrs(refl)
+            return refl
+        else:
+            raise KeyError('Invalid calibration: {}'.format(calibration.name))
+
+    def _counts_to_radiance(self, counts):
+        """Convert VIS counts to radiance.
+
+        Reference: [PUG], equations (7) and (8).
+        """
+        coefs = self.coefs['radiance']
+        years_since_launch = coefs['years_since_launch']
+        a_cf = (coefs['a0_vis'] +
+                coefs['a1_vis'] * years_since_launch +
+                coefs['a2_vis'] * years_since_launch ** 2)
+        mean_count_space_vis = coefs['mean_count_space_vis']
+        rad = (counts - mean_count_space_vis) * a_cf
+        return rad.where(rad > 0, np.float32(np.nan))
+
+    def _radiance_to_reflectance(self, rad):
+        """Convert VIS radiance to reflectance factor.
+
+        Note: Produces huge reflectances in situations where both radiance and
+        solar zenith angle are small. Maybe the corresponding uncertainties
+        can be used to filter these cases before calculating reflectances.
+
+        Reference: [PUG], equation (6).
+        """
+        coefs = self.coefs['reflectance']
+        sza = self.solar_zenith_angle
+        sza = sza.where(da.fabs(sza) < 90,
+                        np.float32(np.nan))  # direct illumination only
+        cos_sza = np.cos(np.deg2rad(sza))
+        distance_sun_earth2 = coefs['distance_sun_earth'] ** 2
+        solar_irradiance_vis = coefs['solar_irradiance_vis']
+        refl = (
+           (np.pi * distance_sun_earth2) /
+           (solar_irradiance_vis * cos_sza) *
+           rad
+        )
+        return self.refl_factor_to_percent(refl)
+
+    def update_refl_attrs(self, refl):
+        """Update attributes of reflectance datasets."""
+        refl.attrs['sun_earth_distance_correction_applied'] = True
+        refl.attrs['sun_earth_distance_correction_factor'] = self.coefs[
+            'reflectance']['distance_sun_earth'].item()
+        return refl
+
+    @staticmethod
+    def refl_factor_to_percent(refl):
+        """Convert reflectance factor to percent."""
+        return refl * 100
+
+
+def interp_tiepoints(ds, target_x, target_y):
+    """Interpolate dataset between tiepoints.
+
+    Uses linear interpolation.
+
+    FUTURE: [PUG] recommends cubic spline interpolation.
+
+    Args:
+        ds:
+            Dataset to be interpolated
+        target_x:
+            Target x coordinates
+        target_y:
+            Target y coordinates
+    """
+    # No tiepoint coordinates specified in the files. Use dimensions
+    # to calculate tiepoint sampling and assign tiepoint coordinates
+    # accordingly.
+    sampling = target_x.size // ds.coords['x'].size
+    ds = ds.assign_coords(x=target_x.values[::sampling],
+                          y=target_y.values[::sampling])
+
+    return ds.interp(x=target_x.values, y=target_y.values)
+
+
+def interp_acq_time(time2d, target_y):
+    """Interpolate scanline acquisition time to the given coordinates.
+
+    The files provide timestamps per pixel for the low resolution
+    channels (IR/WV) only.
+
+    1) Average values in each line to obtain one timestamp per line.
+    2) For the VIS channel duplicate values in y-direction (as
+       advised by [PUG]).
+
+    Note that the timestamps do not increase monotonically with the
+    line number in some cases.
+
+    Returns:
+        Mean scanline acquisition timestamps
+    """
+    # Compute mean timestamp per scanline
+    time = time2d.mean(dim='x_ir_wv').rename({'y_ir_wv': 'y'})
+
+    # If required, repeat timestamps in y-direction to obtain higher
+    # resolution
+    y = time.coords['y'].values
+    if y.size < target_y.size:
+        reps = target_y.size // y.size
+        y_rep = np.repeat(y, reps)
+        time_hires = time.reindex(y=y_rep)
+        time_hires = time_hires.assign_coords(y=target_y)
+        return time_hires
+    return time
+
+
+def is_high_resol(resolution):
+    """Identify high resolution channel."""
+    return resolution == HIGH_RESOL
+
+
 class FiduceoMviriBase(BaseFileHandler):
     """Baseclass for FIDUCEO MVIRI file handlers."""
+
     nc_keys = {
         'WV': 'count_wv',
         'IR': 'count_ir'
@@ -261,7 +446,7 @@ def get_dataset(self, dataset_id, dataset_info):
 
     def get_area_def(self, dataset_id):
         """Get area definition of the given dataset."""
-        if self._is_high_resol(dataset_id['resolution']):
+        if is_high_resol(dataset_id['resolution']):
             im_size = self.nc.coords['y'].size
             area_name = 'geos_mviri_vis'
         else:
@@ -310,7 +495,7 @@ def _get_channel(self, name, resolution, calibration):
             if self.mask_bad_quality:
                 ds = self._mask_vis(ds)
             else:
-                self._check_vis_quality(ds)
+                self._check_vis_quality()
         ds['acq_time'] = ('y', self._get_acq_time(resolution))
         return ds
 
@@ -322,13 +507,13 @@ def _get_angles(self, name, resolution):
         resolution. Interpolate them to the desired resolution.
         """
         angles = self._read_dataset(name)
-        if self._is_high_resol(resolution):
+        if is_high_resol(resolution):
             target_x = self.nc.coords['x']
             target_y = self.nc.coords['y']
         else:
             target_x = self.nc.coords['x_ir_wv']
             target_y = self.nc.coords['y_ir_wv']
-        return self._interp_tiepoints(
+        return interp_tiepoints(
             angles,
             target_x=target_x,
             target_y=target_y
@@ -338,7 +523,7 @@ def _get_other_dataset(self, name):
         """Get other datasets such as uncertainties."""
         ds = self._read_dataset(name)
         if name in OTHER_REFLECTANCES:
-            ds = ds * 100  # conversion to percent
+            ds = VISCalibrator.refl_factor_to_percent(ds)
         return ds
 
     def _get_nc_key(self, ds_name):
@@ -384,57 +569,37 @@ def _calibrate_vis(self, ds, calibration):
         """Calibrate VIS channel. To be implemented by subclasses."""
         raise NotImplementedError
 
-    def _update_refl_attrs(self, refl):
-        """Update attributes of reflectance datasets."""
-        refl.attrs['sun_earth_distance_correction_applied'] = True
-        refl.attrs['sun_earth_distance_correction_factor'] = self.nc[
-            'distance_sun_earth'].item()
-        return refl
-
     def _calibrate_ir_wv(self, ds, channel, calibration):
         """Calibrate IR and WV channel."""
-        if calibration == 'counts':
-            return ds
-        elif calibration in ('radiance', 'brightness_temperature'):
-            rad = self._ir_wv_counts_to_radiance(ds, channel)
-            if calibration == 'radiance':
-                return rad
-            bt = self._ir_wv_radiance_to_brightness_temperature(rad, channel)
-            return bt
-        else:
-            raise KeyError('Invalid calibration: {}'.format(calibration.name))
-
-    def _get_coefs_ir_wv(self, channel, calibration):
-        """Get calibration coefficients for IR/WV channels.
-
-        Returns:
-            Offset (a), Slope (b)
-        """
-        nc_key_a = self.nc_keys_coefs[channel][calibration]['a']
-        nc_key_b = self.nc_keys_coefs[channel][calibration]['b']
-        a = np.float32(self.nc[nc_key_a])
-        b = np.float32(self.nc[nc_key_b])
-        return a, b
-
-    def _ir_wv_counts_to_radiance(self, counts, channel):
-        """Convert IR/WV counts to radiance.
-
-        Reference: [PUG], equations (4.1) and (4.2).
-        """
-        a, b = self._get_coefs_ir_wv(channel, 'radiance')
-        rad = a + b * counts
-        return rad.where(rad > 0, np.float32(np.nan))
-
-    def _ir_wv_radiance_to_brightness_temperature(self, rad, channel):
-        """Convert IR/WV radiance to brightness temperature.
-
-        Reference: [PUG], equations (5.1) and (5.2).
-        """
-        a, b = self._get_coefs_ir_wv(channel, 'brightness_temperature')
-        bt = b / (np.log(rad) - a)
-        return bt.where(bt > 0, np.float32(np.nan))
+        coefs = self._get_coefs_ir_wv(channel)
+        cal = IRWVCalibrator(coefs)
+        return cal.calibrate(ds, calibration)
+
+    def _get_coefs_ir_wv(self, channel):
+        """Get calibration coefficients for IR/WV channels."""
+        coefs = {}
+        for calibration in ('radiance', 'brightness_temperature'):
+            nc_key_a = self.nc_keys_coefs[channel][calibration]['a']
+            nc_key_b = self.nc_keys_coefs[channel][calibration]['b']
+            a = np.float32(self.nc[nc_key_a])
+            b = np.float32(self.nc[nc_key_b])
+            coefs[calibration] = (a, b)
+        return coefs
+
+    def _get_coefs_vis(self):
+        """Get VIS calibration coefs present in both file types."""
+        return {
+            'reflectance': {
+                'distance_sun_earth': np.float32(
+                    self.nc['distance_sun_earth']
+                ),
+                'solar_irradiance_vis': np.float32(
+                    self.nc['solar_irradiance_vis']
+                )
+            }
+        }
 
-    def _check_vis_quality(self, ds):
+    def _check_vis_quality(self):
         """Check VIS channel quality and issue a warning if it's bad."""
         mask = self._read_dataset('quality_pixel_bitmask')
         use_with_caution = da.bitwise_and(mask, 2)
@@ -467,41 +632,11 @@ def _get_acq_time(self, resolution):
             time2d = self.nc['time_ir_wv']
         except KeyError:
             time2d = self.nc['time']
-        if self._is_high_resol(resolution):
+        if is_high_resol(resolution):
             target_y = self.nc.coords['x']
         else:
             target_y = self.nc.coords['x_ir_wv']
-        return self._interp_acq_time(time2d, target_y=target_y.values)
-
-    def _interp_acq_time(self, time2d, target_y):
-        """Interpolate scanline acquisition time to the given coordinates.
-
-        The files provide timestamps per pixel for the low resolution
-        channels (IR/WV) only.
-
-        1) Average values in each line to obtain one timestamp per line.
-        2) For the VIS channel duplicate values in y-direction (as
-           advised by [PUG]).
-
-        Note that the timestamps do not increase monotonically with the
-        line number in some cases.
-
-        Returns:
-            Mean scanline acquisition timestamps
-        """
-        # Compute mean timestamp per scanline
-        time = time2d.mean(dim='x_ir_wv').rename({'y_ir_wv': 'y'})
-
-        # If required, repeat timestamps in y-direction to obtain higher
-        # resolution
-        y = time.coords['y'].values
-        if y.size < target_y.size:
-            reps = target_y.size // y.size
-            y_rep = np.repeat(y, reps)
-            time_hires = time.reindex(y=y_rep)
-            time_hires = time_hires.assign_coords(y=target_y)
-            return time_hires
-        return time
+        return interp_acq_time(time2d, target_y=target_y.values)
 
     def _get_orbital_parameters(self):
         """Get the orbital parameters."""
@@ -546,33 +681,6 @@ def _get_ssp(self, coord):
             sub_lonlat = np.nan
         return sub_lonlat
 
-    def _interp_tiepoints(self, ds, target_x, target_y):
-        """Interpolate dataset between tiepoints.
-
-        Uses linear interpolation.
-
-        FUTURE: [PUG] recommends cubic spline interpolation.
-
-        Args:
-            ds:
-                Dataset to be interpolated
-            target_x:
-                Target x coordinates
-            target_y:
-                Target y coordinates
-        """
-        # No tiepoint coordinates specified in the files. Use dimensions
-        # to calculate tiepoint sampling and assign tiepoint coordinates
-        # accordingly.
-        sampling = target_x.size // ds.coords['x'].size
-        ds = ds.assign_coords(x=target_x.values[::sampling],
-                              y=target_y.values[::sampling])
-
-        return ds.interp(x=target_x.values, y=target_y.values)
-
-    def _is_high_resol(self, resolution):
-        return resolution == HIGH_RESOL
-
 
 class FiduceoMviriEasyFcdrFileHandler(FiduceoMviriBase):
     """File handler for FIDUCEO MVIRI Easy FCDR."""
@@ -586,8 +694,10 @@ def _calibrate_vis(self, ds, calibration):
         Easy FCDR provides reflectance only, no counts or radiance.
         """
         if calibration == 'reflectance':
-            refl = 100 * ds  # conversion to percent
-            refl = self._update_refl_attrs(refl)
+            coefs = self._get_coefs_vis()
+            cal = VISCalibrator(coefs)
+            refl = cal.refl_factor_to_percent(ds)
+            refl = cal.update_refl_attrs(refl)
             return refl
         elif calibration in ('counts', 'radiance'):
             raise ValueError('Cannot calibrate to {}. Easy FCDR provides '
@@ -602,56 +712,26 @@ class FiduceoMviriFullFcdrFileHandler(FiduceoMviriBase):
     nc_keys = FiduceoMviriBase.nc_keys.copy()
     nc_keys['VIS'] = 'count_vis'
 
-    def _calibrate_vis(self, ds, calibration):
-        """Calibrate VIS channel.
-
-        All calibration levels are available here.
-        """
-        if calibration == 'counts':
-            return ds
-        elif calibration in ('radiance', 'reflectance'):
-            rad = self._vis_counts_to_radiance(ds)
-            if calibration == 'radiance':
-                return rad
-            refl = self._vis_radiance_to_reflectance(rad)
-            refl = self._update_refl_attrs(refl)
-            return refl
-        else:
-            raise KeyError('Invalid calibration: {}'.format(calibration.name))
-
-    def _vis_counts_to_radiance(self, counts):
-        """Convert VIS counts to radiance.
-
-        Reference: [PUG], equations (7) and (8).
-        """
-        years_since_launch = self.nc['years_since_launch']
-        a_cf = (self.nc['a0_vis'] +
-                self.nc['a1_vis'] * years_since_launch +
-                self.nc['a2_vis'] * years_since_launch ** 2)
-        mean_count_space_vis = np.float32(self.nc['mean_count_space_vis'])
-        a_cf = np.float32(a_cf)
-        rad = (counts - mean_count_space_vis) * a_cf
-        return rad.where(rad > 0, np.float32(np.nan))
-
-    def _vis_radiance_to_reflectance(self, rad):
-        """Convert VIS radiance to reflectance factor.
-
-        Note: Produces huge reflectances in situations where both radiance and
-        solar zenith angle are small. Maybe the corresponding uncertainties
-        can be used to filter these cases before calculating reflectances.
+    def _get_coefs_vis(self):
+        """Get full set of calibration coefficients for VIS channel."""
+        coefs = super()._get_coefs_vis()
+        coefs.update({
+            'radiance': {
+                'years_since_launch': np.float32(self.nc['years_since_launch']),
+                'a0_vis': np.float32(self.nc['a0_vis']),
+                'a1_vis': np.float32(self.nc['a1_vis']),
+                'a2_vis': np.float32(self.nc['a2_vis']),
+                'mean_count_space_vis': np.float32(
+                    self.nc['mean_count_space_vis']
+                )
+            },
+        })
+        return coefs
 
-        Reference: [PUG], equation (6).
-        """
-        sza = self._get_angles('solar_zenith_angle', HIGH_RESOL)
-        sza = sza.where(da.fabs(sza) < 90,
-                        np.float32(np.nan))  # direct illumination only
-        cos_sza = np.cos(np.deg2rad(sza))
-        distance_sun_earth2 = np.float32(self.nc['distance_sun_earth'] ** 2)
-        solar_irradiance_vis = np.float32(self.nc['solar_irradiance_vis'])
-        refl = (
-           (np.pi * distance_sun_earth2) /
-           (solar_irradiance_vis * cos_sza) *
-           rad
-        )
-        refl = refl * 100  # conversion to percent
-        return refl
+    def _calibrate_vis(self, ds, calibration):
+        coefs = self._get_coefs_vis()
+        sza = None
+        if calibration == 'reflectance':
+            sza = self._get_angles('solar_zenith_angle', HIGH_RESOL)
+        cal = VISCalibrator(coefs, sza)
+        return cal.calibrate(ds, calibration)
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 09571e0ccc..390bdfee44 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -427,9 +427,7 @@ def test_get_dataset(self, file_handler, name, calibration, resolution,
             assert ds.dtype == expected.dtype
             assert ds.attrs == expected.attrs
 
-    @mock.patch(
-        'satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriBase._interp_acq_time'
-    )
+    @mock.patch('satpy.readers.mviri_l1b_fiduceo_nc.interp_acq_time')
     def test_time_cache(self, interp_acq_time, file_handler):
         """Test caching of acquisition times."""
         dataset_id = make_dataid(
@@ -460,9 +458,7 @@ def test_time_cache(self, interp_acq_time, file_handler):
         file_handler.get_dataset(another_id, info)
         interp_acq_time.assert_called()
 
-    @mock.patch(
-        'satpy.readers.mviri_l1b_fiduceo_nc.FiduceoMviriBase._interp_tiepoints'
-    )
+    @mock.patch('satpy.readers.mviri_l1b_fiduceo_nc.interp_tiepoints')
     def test_angle_cache(self, interp_tiepoints, file_handler):
         """Test caching of angle datasets."""
         dataset_id = make_dataid(name='solar_zenith_angle',

From 72042bf922e75c222afdcab730450795022b2ecd Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 25 Nov 2020 12:11:22 +0000
Subject: [PATCH 42/47] Factorize again to reduce complexity even more

---
 satpy/readers/mviri_l1b_fiduceo_nc.py         | 414 +++++++++---------
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py |  16 +-
 2 files changed, 218 insertions(+), 212 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index aedbbb07d4..61c7e0525a 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -208,21 +208,24 @@ def calibrate(self, counts, calibration):
         if calibration == 'counts':
             return counts
         elif calibration in ('radiance', 'brightness_temperature'):
-            rad = self._counts_to_radiance(counts)
-            if calibration == 'radiance':
-                return rad
-            bt = self._radiance_to_brightness_temperature(rad)
-            return bt
+            return self._calibrate_rad_bt(counts, calibration)
         else:
             raise KeyError('Invalid calibration: {}'.format(calibration.name))
 
+    def _calibrate_rad_bt(self, counts, calibration):
+        """Calibrate counts to radiance or brightness temperature."""
+        rad = self._counts_to_radiance(counts)
+        if calibration == 'radiance':
+            return rad
+        bt = self._radiance_to_brightness_temperature(rad)
+        return bt
+
     def _counts_to_radiance(self, counts):
         """Convert IR/WV counts to radiance.
 
         Reference: [PUG], equations (4.1) and (4.2).
         """
-        a, b = self.coefs['radiance']
-        rad = a + b * counts
+        rad = self.coefs['a'] + self.coefs['b'] * counts
         return rad.where(rad > 0, np.float32(np.nan))
 
     def _radiance_to_brightness_temperature(self, rad):
@@ -230,8 +233,7 @@ def _radiance_to_brightness_temperature(self, rad):
 
         Reference: [PUG], equations (5.1) and (5.2).
         """
-        a, b = self.coefs['brightness_temperature']
-        bt = b / (np.log(rad) - a)
+        bt = self.coefs['bt_b'] / (np.log(rad) - self.coefs['bt_a'])
         return bt.where(bt > 0, np.float32(np.nan))
 
 
@@ -256,26 +258,29 @@ def calibrate(self, counts, calibration):
         if calibration == 'counts':
             return counts
         elif calibration in ('radiance', 'reflectance'):
-            rad = self._counts_to_radiance(counts)
-            if calibration == 'radiance':
-                return rad
-            refl = self._radiance_to_reflectance(rad)
-            refl = self.update_refl_attrs(refl)
-            return refl
+            return self._calibrate_rad_refl(counts, calibration)
         else:
             raise KeyError('Invalid calibration: {}'.format(calibration.name))
 
+    def _calibrate_rad_refl(self, counts, calibration):
+        """Calibrate counts to radiance or reflectance."""
+        rad = self._counts_to_radiance(counts)
+        if calibration == 'radiance':
+            return rad
+        refl = self._radiance_to_reflectance(rad)
+        refl = self.update_refl_attrs(refl)
+        return refl
+
     def _counts_to_radiance(self, counts):
         """Convert VIS counts to radiance.
 
         Reference: [PUG], equations (7) and (8).
         """
-        coefs = self.coefs['radiance']
-        years_since_launch = coefs['years_since_launch']
-        a_cf = (coefs['a0_vis'] +
-                coefs['a1_vis'] * years_since_launch +
-                coefs['a2_vis'] * years_since_launch ** 2)
-        mean_count_space_vis = coefs['mean_count_space_vis']
+        years_since_launch = self.coefs['years_since_launch']
+        a_cf = (self.coefs['a0'] +
+                self.coefs['a1'] * years_since_launch +
+                self.coefs['a2'] * years_since_launch ** 2)
+        mean_count_space_vis = self.coefs['mean_count_space']
         rad = (counts - mean_count_space_vis) * a_cf
         return rad.where(rad > 0, np.float32(np.nan))
 
@@ -288,16 +293,14 @@ def _radiance_to_reflectance(self, rad):
 
         Reference: [PUG], equation (6).
         """
-        coefs = self.coefs['reflectance']
-        sza = self.solar_zenith_angle
-        sza = sza.where(da.fabs(sza) < 90,
-                        np.float32(np.nan))  # direct illumination only
+        sza = self.solar_zenith_angle.where(
+            da.fabs(self.solar_zenith_angle) < 90,
+            np.float32(np.nan)
+        )  # direct illumination only
         cos_sza = np.cos(np.deg2rad(sza))
-        distance_sun_earth2 = coefs['distance_sun_earth'] ** 2
-        solar_irradiance_vis = coefs['solar_irradiance_vis']
         refl = (
-           (np.pi * distance_sun_earth2) /
-           (solar_irradiance_vis * cos_sza) *
+           (np.pi * self.coefs['distance_sun_earth'] ** 2) /
+           (self.coefs['solar_irradiance'] * cos_sza) *
            rad
         )
         return self.refl_factor_to_percent(refl)
@@ -306,7 +309,7 @@ def update_refl_attrs(self, refl):
         """Update attributes of reflectance datasets."""
         refl.attrs['sun_earth_distance_correction_applied'] = True
         refl.attrs['sun_earth_distance_correction_factor'] = self.coefs[
-            'reflectance']['distance_sun_earth'].item()
+            'distance_sun_earth'].item()
         return refl
 
     @staticmethod
@@ -315,60 +318,133 @@ def refl_factor_to_percent(refl):
         return refl * 100
 
 
-def interp_tiepoints(ds, target_x, target_y):
-    """Interpolate dataset between tiepoints.
+class Navigator:
+    """Navigate MVIRI images."""
+
+    def get_area_def(self, im_size, projection_longitude):
+        """Create MVIRI area definition."""
+        proj_params = self._get_proj_params(im_size, projection_longitude)
+        extent = get_area_extent(proj_params)
+        return get_area_definition(proj_params, extent)
+
+    def _get_proj_params(self, im_size, projection_longitude):
+        """Get projection parameters for the given settings."""
+        area_name = 'geos_mviri_{0}x{0}'.format(im_size)
+        lfac, cfac, loff, coff = self._get_factors_offsets(im_size)
+        return {
+            'ssp_lon': projection_longitude,
+            'a': EQUATOR_RADIUS,
+            'b': POLE_RADIUS,
+            'h': ALTITUDE,
+            'units': 'm',
+            'loff': loff - im_size,
+            'coff': coff,
+            'lfac': -lfac,
+            'cfac': -cfac,
+            'nlines': im_size,
+            'ncols': im_size,
+            'scandir': 'S2N',  # Reference: [PUG] section 2.
+            'p_id': area_name,
+            'a_name': area_name,
+            'a_desc': 'MVIRI Geostationary Projection'
+        }
+
+    def _get_factors_offsets(self, im_size):
+        """Determine line/column offsets and scaling factors."""
+        # For offsets see variables "asamp" and "aline" of subroutine
+        # "refgeo" in [Handbook] and in
+        # https://github.com/FIDUCEO/FCDR_MVIRI/blob/master/lib/nrCrunch/cruncher.f
+        loff = coff = im_size / 2 + 0.5
+        lfac = cfac = sampling_to_lfac_cfac(
+            np.deg2rad(MVIRI_FIELD_OF_VIEW) / im_size
+        )
+        return lfac, cfac, loff, coff
+
+
+class Interpolator:
+    """Interpolate datasets to another resolution."""
+
+    @staticmethod
+    def interp_tiepoints(ds, target_x, target_y):
+        """Interpolate dataset between tiepoints.
+
+        Uses linear interpolation.
 
-    Uses linear interpolation.
+        FUTURE: [PUG] recommends cubic spline interpolation.
+
+        Args:
+            ds:
+                Dataset to be interpolated
+            target_x:
+                Target x coordinates
+            target_y:
+                Target y coordinates
+        """
+        # No tiepoint coordinates specified in the files. Use dimensions
+        # to calculate tiepoint sampling and assign tiepoint coordinates
+        # accordingly.
+        sampling = target_x.size // ds.coords['x'].size
+        ds = ds.assign_coords(x=target_x.values[::sampling],
+                              y=target_y.values[::sampling])
+
+        return ds.interp(x=target_x.values, y=target_y.values)
+
+    @staticmethod
+    def interp_acq_time(time2d, target_y):
+        """Interpolate scanline acquisition time to the given coordinates.
 
-    FUTURE: [PUG] recommends cubic spline interpolation.
+        The files provide timestamps per pixel for the low resolution
+        channels (IR/WV) only.
 
-    Args:
-        ds:
-            Dataset to be interpolated
-        target_x:
-            Target x coordinates
-        target_y:
-            Target y coordinates
-    """
-    # No tiepoint coordinates specified in the files. Use dimensions
-    # to calculate tiepoint sampling and assign tiepoint coordinates
-    # accordingly.
-    sampling = target_x.size // ds.coords['x'].size
-    ds = ds.assign_coords(x=target_x.values[::sampling],
-                          y=target_y.values[::sampling])
+        1) Average values in each line to obtain one timestamp per line.
+        2) For the VIS channel duplicate values in y-direction (as
+           advised by [PUG]).
 
-    return ds.interp(x=target_x.values, y=target_y.values)
+        Note that the timestamps do not increase monotonically with the
+        line number in some cases.
 
+        Returns:
+            Mean scanline acquisition timestamps
+        """
+        # Compute mean timestamp per scanline
+        time = time2d.mean(dim='x_ir_wv').rename({'y_ir_wv': 'y'})
 
-def interp_acq_time(time2d, target_y):
-    """Interpolate scanline acquisition time to the given coordinates.
+        # If required, repeat timestamps in y-direction to obtain higher
+        # resolution
+        y = time.coords['y'].values
+        if y.size < target_y.size:
+            reps = target_y.size // y.size
+            y_rep = np.repeat(y, reps)
+            time_hires = time.reindex(y=y_rep)
+            time_hires = time_hires.assign_coords(y=target_y)
+            return time_hires
+        return time
 
-    The files provide timestamps per pixel for the low resolution
-    channels (IR/WV) only.
 
-    1) Average values in each line to obtain one timestamp per line.
-    2) For the VIS channel duplicate values in y-direction (as
-       advised by [PUG]).
+class VisQualityControl:
+    """Simple quality control for VIS channel."""
 
-    Note that the timestamps do not increase monotonically with the
-    line number in some cases.
+    def __init__(self, mask):
+        """Initialize the quality control."""
+        self._mask = mask
 
-    Returns:
-        Mean scanline acquisition timestamps
-    """
-    # Compute mean timestamp per scanline
-    time = time2d.mean(dim='x_ir_wv').rename({'y_ir_wv': 'y'})
+    def check(self):
+        """Check VIS channel quality and issue a warning if it's bad."""
+        use_with_caution = da.bitwise_and(self._mask, 2)
+        if use_with_caution.all():
+            warnings.warn(
+                'All pixels of the VIS channel are flagged as "use with '
+                'caution". Use datasets "quality_pixel_bitmask" and '
+                '"data_quality_bitmask" to find out why.'
+            )
 
-    # If required, repeat timestamps in y-direction to obtain higher
-    # resolution
-    y = time.coords['y'].values
-    if y.size < target_y.size:
-        reps = target_y.size // y.size
-        y_rep = np.repeat(y, reps)
-        time_hires = time.reindex(y=y_rep)
-        time_hires = time_hires.assign_coords(y=target_y)
-        return time_hires
-    return time
+    def mask(self, ds):
+        """Mask VIS pixels with bad quality.
+
+        Pixels are considered bad quality if the "quality_pixel_bitmask" is
+        everything else than 0 (no flag set).
+        """
+        return ds.where(self._mask == 0, np.float32(np.nan))
 
 
 def is_high_resol(resolution):
@@ -383,28 +459,6 @@ class FiduceoMviriBase(BaseFileHandler):
         'WV': 'count_wv',
         'IR': 'count_ir'
     }
-    nc_keys_coefs = {
-        'WV': {
-            'radiance': {
-                'a': 'a_wv',
-                'b': 'b_wv'
-            },
-            'brightness_temperature': {
-                'a': 'bt_a_wv',
-                'b': 'bt_b_wv'
-            }
-        },
-        'IR': {
-            'radiance': {
-                'a': 'a_ir',
-                'b': 'b_ir'
-            },
-            'brightness_temperature': {
-                'a': 'bt_a_ir',
-                'b': 'bt_b_ir'
-            }
-        },
-    }
 
     def __init__(self, filename, filename_info, filetype_info,
                  mask_bad_quality=False):
@@ -430,6 +484,7 @@ def __init__(self, filename, filename_info, filetype_info,
         # Projection longitude is not provided in the file, read it from the
         # filename.
         self.projection_longitude = float(filename_info['projection_longitude'])
+        self.calib_coefs = self._get_calib_coefs()
 
     def get_dataset(self, dataset_id, dataset_info):
         """Get the dataset."""
@@ -448,41 +503,14 @@ def get_area_def(self, dataset_id):
         """Get area definition of the given dataset."""
         if is_high_resol(dataset_id['resolution']):
             im_size = self.nc.coords['y'].size
-            area_name = 'geos_mviri_vis'
         else:
             im_size = self.nc.coords['y_ir_wv'].size
-            area_name = 'geos_mviri_ir_wv'
-
-        # Determine line/column offsets and scaling factors. For offsets
-        # see variables "asamp" and "aline" of subroutine "refgeo" in
-        # [Handbook] and in
-        # https://github.com/FIDUCEO/FCDR_MVIRI/blob/master/lib/nrCrunch/cruncher.f
-        loff = coff = im_size / 2 + 0.5
-        lfac = cfac = sampling_to_lfac_cfac(
-            np.deg2rad(MVIRI_FIELD_OF_VIEW) / im_size
+        nav = Navigator()
+        return nav.get_area_def(
+            im_size=im_size,
+            projection_longitude=self.projection_longitude
         )
 
-        pdict = {
-            'ssp_lon': self.projection_longitude,
-            'a': EQUATOR_RADIUS,
-            'b': POLE_RADIUS,
-            'h': ALTITUDE,
-            'units': 'm',
-            'loff': loff - im_size,
-            'coff': coff,
-            'lfac': -lfac,
-            'cfac': -cfac,
-            'nlines': im_size,
-            'ncols': im_size,
-            'scandir': 'S2N',  # Reference: [PUG] section 2.
-            'p_id': area_name,
-            'a_name': area_name,
-            'a_desc': 'MVIRI Geostationary Projection'
-        }
-        extent = get_area_extent(pdict)
-        area_def = get_area_definition(pdict, extent)
-        return area_def
-
     def _get_channel(self, name, resolution, calibration):
         """Get and calibrate channel data."""
         ds = self._read_dataset(name)
@@ -492,10 +520,11 @@ def _get_channel(self, name, resolution, calibration):
             calibration=calibration
         )
         if name == 'VIS':
+            qc = VisQualityControl(self._read_dataset('quality_pixel_bitmask'))
             if self.mask_bad_quality:
-                ds = self._mask_vis(ds)
+                ds = qc.mask(ds)
             else:
-                self._check_vis_quality()
+                qc.check()
         ds['acq_time'] = ('y', self._get_acq_time(resolution))
         return ds
 
@@ -513,7 +542,7 @@ def _get_angles(self, name, resolution):
         else:
             target_x = self.nc.coords['x_ir_wv']
             target_y = self.nc.coords['y_ir_wv']
-        return interp_tiepoints(
+        return Interpolator.interp_tiepoints(
             angles,
             target_x=target_x,
             target_y=target_y
@@ -526,13 +555,9 @@ def _get_other_dataset(self, name):
             ds = VISCalibrator.refl_factor_to_percent(ds)
         return ds
 
-    def _get_nc_key(self, ds_name):
-        """Get netCDF variable name for the given dataset."""
-        return self.nc_keys.get(ds_name, ds_name)
-
     def _read_dataset(self, name):
         """Read a dataset from the file."""
-        nc_key = self._get_nc_key(name)
+        nc_key = self.nc_keys.get(name, name)
         ds = self.nc[nc_key]
         if 'y_ir_wv' in ds.dims:
             ds = ds.rename({'y_ir_wv': 'y', 'x_ir_wv': 'x'})
@@ -543,6 +568,7 @@ def _read_dataset(self, name):
             # high resolution data variables.
             ds = ds.assign_coords({'y': self.nc.coords['y'],
                                    'x': self.nc.coords['x']})
+        ds.attrs.pop('ancillary_variables', None)  # to avoid satpy warnings
         return ds
 
     def _update_attrs(self, ds, info):
@@ -555,70 +581,48 @@ def _update_attrs(self, ds, info):
 
     def _calibrate(self, ds, channel, calibration):
         """Calibrate the given dataset."""
-        ds.attrs.pop('ancillary_variables', None)  # to avoid satpy warnings
         if channel == 'VIS':
-            return self._calibrate_vis(ds, calibration)
-        elif channel in ['WV', 'IR']:
-            return self._calibrate_ir_wv(ds, channel, calibration)
+            return self._calibrate_vis(ds, channel, calibration)
         else:
-            raise KeyError('Don\'t know how to calibrate channel {}'.format(
-                channel))
+            calib = IRWVCalibrator(self.calib_coefs[channel])
+            return calib.calibrate(ds, calibration)
 
     @abc.abstractmethod
-    def _calibrate_vis(self, ds, calibration):
+    def _calibrate_vis(self, ds, channel, calibration):
         """Calibrate VIS channel. To be implemented by subclasses."""
         raise NotImplementedError
 
-    def _calibrate_ir_wv(self, ds, channel, calibration):
-        """Calibrate IR and WV channel."""
-        coefs = self._get_coefs_ir_wv(channel)
-        cal = IRWVCalibrator(coefs)
-        return cal.calibrate(ds, calibration)
+    def _get_calib_coefs(self):
+        """Get calibration coefficients for all channels.
 
-    def _get_coefs_ir_wv(self, channel):
-        """Get calibration coefficients for IR/WV channels."""
-        coefs = {}
-        for calibration in ('radiance', 'brightness_temperature'):
-            nc_key_a = self.nc_keys_coefs[channel][calibration]['a']
-            nc_key_b = self.nc_keys_coefs[channel][calibration]['b']
-            a = np.float32(self.nc[nc_key_a])
-            b = np.float32(self.nc[nc_key_b])
-            coefs[calibration] = (a, b)
-        return coefs
-
-    def _get_coefs_vis(self):
-        """Get VIS calibration coefs present in both file types."""
-        return {
-            'reflectance': {
-                'distance_sun_earth': np.float32(
-                    self.nc['distance_sun_earth']
-                ),
-                'solar_irradiance_vis': np.float32(
-                    self.nc['solar_irradiance_vis']
-                )
-            }
+        Note: Only coefficients present in both file types.
+        """
+        coefs = {
+            'VIS': {
+                'distance_sun_earth': self.nc['distance_sun_earth'],
+                'solar_irradiance': self.nc['solar_irradiance_vis']
+            },
+            'IR': {
+                'a': self.nc['a_ir'],
+                'b': self.nc['b_ir'],
+                'bt_a': self.nc['bt_a_ir'],
+                'bt_b': self.nc['bt_b_ir']
+            },
+            'WV': {
+                'a': self.nc['a_wv'],
+                'b': self.nc['b_wv'],
+                'bt_a': self.nc['bt_a_wv'],
+                'bt_b': self.nc['bt_b_wv']
+            },
         }
 
-    def _check_vis_quality(self):
-        """Check VIS channel quality and issue a warning if it's bad."""
-        mask = self._read_dataset('quality_pixel_bitmask')
-        use_with_caution = da.bitwise_and(mask, 2)
-        if use_with_caution.all():
-            warnings.warn(
-                'All pixels of the VIS channel are flagged as "use with '
-                'caution". Use datasets "quality_pixel_bitmask" and '
-                '"data_quality_bitmask" to find out why.'
-            )
-
-    def _mask_vis(self, ds):
-        """Mask VIS pixels with bad quality.
+        # Convert coefficients to 32bit float to reduce memory footprint
+        # of calibrated data.
+        for ch in coefs:
+            for name in coefs[ch]:
+                coefs[ch][name] = np.float32(coefs[ch][name])
 
-        Pixels are considered bad quality if the "quality_pixel_bitmask" is
-        everything else than 0 (no flag set).
-        """
-        mask = self._read_dataset('quality_pixel_bitmask')
-        return ds.where(mask == 0,
-                        np.float32(np.nan))
+        return coefs
 
     @lru_cache(maxsize=3)  # Three channels
     def _get_acq_time(self, resolution):
@@ -636,7 +640,7 @@ def _get_acq_time(self, resolution):
             target_y = self.nc.coords['x']
         else:
             target_y = self.nc.coords['x_ir_wv']
-        return interp_acq_time(time2d, target_y=target_y.values)
+        return Interpolator.interp_acq_time(time2d, target_y=target_y.values)
 
     def _get_orbital_parameters(self):
         """Get the orbital parameters."""
@@ -688,13 +692,13 @@ class FiduceoMviriEasyFcdrFileHandler(FiduceoMviriBase):
     nc_keys = FiduceoMviriBase.nc_keys.copy()
     nc_keys['VIS'] = 'toa_bidirectional_reflectance_vis'
 
-    def _calibrate_vis(self, ds, calibration):
+    def _calibrate_vis(self, ds, channel, calibration):
         """Calibrate VIS channel.
 
         Easy FCDR provides reflectance only, no counts or radiance.
         """
         if calibration == 'reflectance':
-            coefs = self._get_coefs_vis()
+            coefs = self.calib_coefs[channel]
             cal = VISCalibrator(coefs)
             refl = cal.refl_factor_to_percent(ds)
             refl = cal.update_refl_attrs(refl)
@@ -712,26 +716,24 @@ class FiduceoMviriFullFcdrFileHandler(FiduceoMviriBase):
     nc_keys = FiduceoMviriBase.nc_keys.copy()
     nc_keys['VIS'] = 'count_vis'
 
-    def _get_coefs_vis(self):
-        """Get full set of calibration coefficients for VIS channel."""
-        coefs = super()._get_coefs_vis()
-        coefs.update({
-            'radiance': {
-                'years_since_launch': np.float32(self.nc['years_since_launch']),
-                'a0_vis': np.float32(self.nc['a0_vis']),
-                'a1_vis': np.float32(self.nc['a1_vis']),
-                'a2_vis': np.float32(self.nc['a2_vis']),
-                'mean_count_space_vis': np.float32(
-                    self.nc['mean_count_space_vis']
-                )
-            },
+    def _get_calib_coefs(self):
+        """Add additional VIS coefficients only present in full FCDR."""
+        coefs = super()._get_calib_coefs()
+        coefs['VIS'].update({
+            'years_since_launch': np.float32(self.nc['years_since_launch']),
+            'a0': np.float32(self.nc['a0_vis']),
+            'a1': np.float32(self.nc['a1_vis']),
+            'a2': np.float32(self.nc['a2_vis']),
+            'mean_count_space': np.float32(
+                self.nc['mean_count_space_vis']
+            )
         })
         return coefs
 
-    def _calibrate_vis(self, ds, calibration):
-        coefs = self._get_coefs_vis()
+    def _calibrate_vis(self, ds, channel, calibration):
+        """Calibrate VIS channel."""
         sza = None
         if calibration == 'reflectance':
             sza = self._get_angles('solar_zenith_angle', HIGH_RESOL)
-        cal = VISCalibrator(coefs, sza)
+        cal = VISCalibrator(self.calib_coefs[channel], sza)
         return cal.calibrate(ds, calibration)
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 390bdfee44..82a84b5944 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -250,8 +250,8 @@
     attrs=attrs_exp
 )
 area_vis_exp = AreaDefinition(
-    area_id='geos_mviri_vis',
-    proj_id='geos_mviri_vis',
+    area_id='geos_mviri_4x4',
+    proj_id='geos_mviri_4x4',
     description='MVIRI Geostationary Projection',
     projection={
         'proj': 'geos',
@@ -265,8 +265,8 @@
     area_extent=[5621229.74392, 5621229.74392, -5621229.74392, -5621229.74392]
 )
 area_ir_wv_exp = area_vis_exp.copy(
-    area_id='geos_mviri_ir_wv',
-    proj_id='geos_mviri_ir_wv',
+    area_id='geos_mviri_2x2',
+    proj_id='geos_mviri_2x2',
     width=2,
     height=2
 )
@@ -427,7 +427,9 @@ def test_get_dataset(self, file_handler, name, calibration, resolution,
             assert ds.dtype == expected.dtype
             assert ds.attrs == expected.attrs
 
-    @mock.patch('satpy.readers.mviri_l1b_fiduceo_nc.interp_acq_time')
+    @mock.patch(
+        'satpy.readers.mviri_l1b_fiduceo_nc.Interpolator.interp_acq_time'
+    )
     def test_time_cache(self, interp_acq_time, file_handler):
         """Test caching of acquisition times."""
         dataset_id = make_dataid(
@@ -458,7 +460,9 @@ def test_time_cache(self, interp_acq_time, file_handler):
         file_handler.get_dataset(another_id, info)
         interp_acq_time.assert_called()
 
-    @mock.patch('satpy.readers.mviri_l1b_fiduceo_nc.interp_tiepoints')
+    @mock.patch(
+        'satpy.readers.mviri_l1b_fiduceo_nc.Interpolator.interp_tiepoints'
+    )
     def test_angle_cache(self, interp_tiepoints, file_handler):
         """Test caching of angle datasets."""
         dataset_id = make_dataid(name='solar_zenith_angle',

From bc7b749d72e054c14ca4acce748daebd53f4b9f5 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 25 Nov 2020 12:13:11 +0000
Subject: [PATCH 43/47] Remove unnecessary else

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 61c7e0525a..d2b15e1c88 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -583,9 +583,8 @@ def _calibrate(self, ds, channel, calibration):
         """Calibrate the given dataset."""
         if channel == 'VIS':
             return self._calibrate_vis(ds, channel, calibration)
-        else:
-            calib = IRWVCalibrator(self.calib_coefs[channel])
-            return calib.calibrate(ds, calibration)
+        calib = IRWVCalibrator(self.calib_coefs[channel])
+        return calib.calibrate(ds, calibration)
 
     @abc.abstractmethod
     def _calibrate_vis(self, ds, channel, calibration):

From 1b504699f394d0b4a9ac4df98fd1f954e72c7059 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 25 Nov 2020 12:42:12 +0000
Subject: [PATCH 44/47] Resolve more test misses

---
 satpy/readers/mviri_l1b_fiduceo_nc.py         |  8 +++--
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py | 36 +++++++++++++++++--
 2 files changed, 40 insertions(+), 4 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index d2b15e1c88..0d4d65338f 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -210,7 +210,9 @@ def calibrate(self, counts, calibration):
         elif calibration in ('radiance', 'brightness_temperature'):
             return self._calibrate_rad_bt(counts, calibration)
         else:
-            raise KeyError('Invalid calibration: {}'.format(calibration.name))
+            raise KeyError(
+                'Invalid IR/WV calibration: {}'.format(calibration.name)
+            )
 
     def _calibrate_rad_bt(self, counts, calibration):
         """Calibrate counts to radiance or brightness temperature."""
@@ -260,7 +262,9 @@ def calibrate(self, counts, calibration):
         elif calibration in ('radiance', 'reflectance'):
             return self._calibrate_rad_refl(counts, calibration)
         else:
-            raise KeyError('Invalid calibration: {}'.format(calibration.name))
+            raise KeyError(
+                'Invalid VIS calibration: {}'.format(calibration.name)
+            )
 
     def _calibrate_rad_refl(self, counts, calibration):
         """Calibrate counts to radiance or reflectance."""
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 82a84b5944..42783b6fb3 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -427,6 +427,26 @@ def test_get_dataset(self, file_handler, name, calibration, resolution,
             assert ds.dtype == expected.dtype
             assert ds.attrs == expected.attrs
 
+    def test_get_dataset_corrupt(self, file_handler):
+        """Test getting datasets with known corruptions."""
+        # Time may have different names and satellite position might be missing
+        file_handler.nc = file_handler.nc.rename(
+            {'time_ir_wv': 'time'}
+        )
+        file_handler.nc = file_handler.nc.drop_vars(
+            ['sub_satellite_longitude_start']
+        )
+
+        dataset_id = make_dataid(
+            name='VIS',
+            calibration='reflectance',
+            resolution=2250
+        )
+        ds = file_handler.get_dataset(dataset_id, {'platform': 'MET7'})
+        assert 'actual_satellite_longitude' not in ds.attrs['orbital_parameters']
+        assert 'actual_satellite_latitude' not in ds.attrs['orbital_parameters']
+        xr.testing.assert_allclose(ds, vis_refl_exp)
+
     @mock.patch(
         'satpy.readers.mviri_l1b_fiduceo_nc.Interpolator.interp_acq_time'
     )
@@ -520,14 +540,26 @@ def test_calib_exceptions(self, file_handler):
             )
         with pytest.raises(KeyError):
             file_handler.get_dataset(
-                {'name': 'VIS', 'calibration': 'invalid'},
+                make_dataid(
+                    name='VIS',
+                    resolution=2250,
+                    calibration='brightness_temperature'),
                 {}
             )
         with pytest.raises(KeyError):
             file_handler.get_dataset(
-                {'name': 'IR', 'calibration': 'invalid'},
+                make_dataid(
+                    name='IR',
+                    resolution=4500,
+                    calibration='reflectance'),
                 {}
             )
+        if isinstance(file_handler, FiduceoMviriEasyFcdrFileHandler):
+            with pytest.raises(KeyError):
+                file_handler.get_dataset(
+                    {'name': 'VIS', 'calibration': 'counts'},
+                    {}
+                )  # not available in easy FCDR
 
     @pytest.mark.file_handler_data(mask_bad_quality=False)
     def test_bad_quality_warning(self, file_handler):

From ba1ee858e67c3d74724cf3eb1295f66d45817c28 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Wed, 25 Nov 2020 20:18:28 +0000
Subject: [PATCH 45/47] Add wrapper class for accessing the dataset

---
 satpy/readers/mviri_l1b_fiduceo_nc.py         | 132 ++++++++++++------
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py |   6 +-
 2 files changed, 94 insertions(+), 44 deletions(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 0d4d65338f..6ab18cf812 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -411,7 +411,7 @@ def interp_acq_time(time2d, target_y):
             Mean scanline acquisition timestamps
         """
         # Compute mean timestamp per scanline
-        time = time2d.mean(dim='x_ir_wv').rename({'y_ir_wv': 'y'})
+        time = time2d.mean(dim='x')
 
         # If required, repeat timestamps in y-direction to obtain higher
         # resolution
@@ -456,6 +456,86 @@ def is_high_resol(resolution):
     return resolution == HIGH_RESOL
 
 
+class DatasetWrapper:
+    """Helper class for accessing the dataset."""
+
+    def __init__(self, nc):
+        """Wrap the given dataset."""
+        self.nc = nc
+
+    @property
+    def attrs(self):
+        """Exposes dataset attributes."""
+        return self.nc.attrs
+
+    def __getitem__(self, item):
+        """Get a variable from the dataset."""
+        ds = self.nc[item]
+        if self._should_dims_be_renamed(ds):
+            ds = self._rename_dims(ds)
+        elif self._coordinates_not_assigned(ds):
+            ds = self._reassign_coords(ds)
+        self._cleanup_attrs(ds)
+        return ds
+
+    def _should_dims_be_renamed(self, ds):
+        """Determine whether dataset dimensions need to be renamed."""
+        return 'y_ir_wv' in ds.dims or 'y_tie' in ds.dims
+
+    def _rename_dims(self, ds):
+        """Rename dataset dimensions to match satpy's expectations."""
+        new_names = {
+            'y_ir_wv': 'y',
+            'x_ir_wv': 'x',
+            'y_tie': 'y',
+            'x_tie': 'x'
+        }
+        for old_name, new_name in new_names.items():
+            if old_name in ds.dims:
+                ds = ds.rename({old_name: new_name})
+        return ds
+
+    def _coordinates_not_assigned(self, ds):
+        return 'y' in ds.dims and 'y' not in ds.coords
+
+    def _reassign_coords(self, ds):
+        """Re-assign coordinates.
+
+        For some reason xarray doesn't assign coordinates to all high
+        resolution data variables.
+        """
+        return ds.assign_coords({'y': self.nc.coords['y'],
+                                 'x': self.nc.coords['x']})
+
+    def _cleanup_attrs(self, ds):
+        """Cleanup dataset attributes."""
+        # Remove ancillary_variables attribute to avoid downstream
+        # satpy warnings.
+        ds.attrs.pop('ancillary_variables', None)
+
+    def get_time(self):
+        """Get time coordinate.
+
+        Variable is sometimes named "time" and sometimes "time_ir_wv".
+        """
+        try:
+            return self['time_ir_wv']
+        except KeyError:
+            return self['time']
+
+    def get_xy_coords(self, resolution):
+        """Get x and y coordinates for the given resolution."""
+        if is_high_resol(resolution):
+            return self.nc.coords['x'], self.nc.coords['y']
+        return self.nc.coords['x_ir_wv'], self.nc.coords['x_ir_wv']
+
+    def get_image_size(self, resolution):
+        """Get image size for the given resolution."""
+        if is_high_resol(resolution):
+            return self.nc.coords['y'].size
+        return self.nc.coords['y_ir_wv'].size
+
+
 class FiduceoMviriBase(BaseFileHandler):
     """Baseclass for FIDUCEO MVIRI file handlers."""
 
@@ -477,13 +557,14 @@ def __init__(self, filename, filename_info, filetype_info,
         super(FiduceoMviriBase, self).__init__(
             filename, filename_info, filetype_info)
         self.mask_bad_quality = mask_bad_quality
-        self.nc = xr.open_dataset(
+        nc_raw = xr.open_dataset(
             filename,
             chunks={'x': CHUNK_SIZE,
                     'y': CHUNK_SIZE,
                     'x_ir_wv': CHUNK_SIZE,
                     'y_ir_wv': CHUNK_SIZE}
         )
+        self.nc = DatasetWrapper(nc_raw)
 
         # Projection longitude is not provided in the file, read it from the
         # filename.
@@ -505,10 +586,7 @@ def get_dataset(self, dataset_id, dataset_info):
 
     def get_area_def(self, dataset_id):
         """Get area definition of the given dataset."""
-        if is_high_resol(dataset_id['resolution']):
-            im_size = self.nc.coords['y'].size
-        else:
-            im_size = self.nc.coords['y_ir_wv'].size
+        im_size = self.nc.get_image_size(dataset_id['resolution'])
         nav = Navigator()
         return nav.get_area_def(
             im_size=im_size,
@@ -517,14 +595,14 @@ def get_area_def(self, dataset_id):
 
     def _get_channel(self, name, resolution, calibration):
         """Get and calibrate channel data."""
-        ds = self._read_dataset(name)
+        ds = self.nc[self.nc_keys[name]]
         ds = self._calibrate(
             ds,
             channel=name,
             calibration=calibration
         )
         if name == 'VIS':
-            qc = VisQualityControl(self._read_dataset('quality_pixel_bitmask'))
+            qc = VisQualityControl(self.nc['quality_pixel_bitmask'])
             if self.mask_bad_quality:
                 ds = qc.mask(ds)
             else:
@@ -539,13 +617,8 @@ def _get_angles(self, name, resolution):
         Files provide angles (solar/satellite zenith & azimuth) at a coarser
         resolution. Interpolate them to the desired resolution.
         """
-        angles = self._read_dataset(name)
-        if is_high_resol(resolution):
-            target_x = self.nc.coords['x']
-            target_y = self.nc.coords['y']
-        else:
-            target_x = self.nc.coords['x_ir_wv']
-            target_y = self.nc.coords['y_ir_wv']
+        angles = self.nc[name]
+        target_x, target_y = self.nc.get_xy_coords(resolution)
         return Interpolator.interp_tiepoints(
             angles,
             target_x=target_x,
@@ -554,27 +627,11 @@ def _get_angles(self, name, resolution):
 
     def _get_other_dataset(self, name):
         """Get other datasets such as uncertainties."""
-        ds = self._read_dataset(name)
+        ds = self.nc[name]
         if name in OTHER_REFLECTANCES:
             ds = VISCalibrator.refl_factor_to_percent(ds)
         return ds
 
-    def _read_dataset(self, name):
-        """Read a dataset from the file."""
-        nc_key = self.nc_keys.get(name, name)
-        ds = self.nc[nc_key]
-        if 'y_ir_wv' in ds.dims:
-            ds = ds.rename({'y_ir_wv': 'y', 'x_ir_wv': 'x'})
-        elif 'y_tie' in ds.dims:
-            ds = ds.rename({'y_tie': 'y', 'x_tie': 'x'})
-        elif 'y' in ds.dims and 'y' not in ds.coords:
-            # For some reason xarray doesn't assign coordinates to all
-            # high resolution data variables.
-            ds = ds.assign_coords({'y': self.nc.coords['y'],
-                                   'x': self.nc.coords['x']})
-        ds.attrs.pop('ancillary_variables', None)  # to avoid satpy warnings
-        return ds
-
     def _update_attrs(self, ds, info):
         """Update dataset attributes."""
         ds.attrs.update(info)
@@ -634,15 +691,8 @@ def _get_acq_time(self, resolution):
         Note that the acquisition time does not increase monotonically
         with the scanline number due to the scan pattern and rectification.
         """
-        # Variable is sometimes named "time" and sometimes "time_ir_wv".
-        try:
-            time2d = self.nc['time_ir_wv']
-        except KeyError:
-            time2d = self.nc['time']
-        if is_high_resol(resolution):
-            target_y = self.nc.coords['x']
-        else:
-            target_y = self.nc.coords['x_ir_wv']
+        time2d = self.nc.get_time()
+        _, target_y = self.nc.get_xy_coords(resolution)
         return Interpolator.interp_acq_time(time2d, target_y=target_y.values)
 
     def _get_orbital_parameters(self):
diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 42783b6fb3..5e8d1c29b2 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -430,10 +430,10 @@ def test_get_dataset(self, file_handler, name, calibration, resolution,
     def test_get_dataset_corrupt(self, file_handler):
         """Test getting datasets with known corruptions."""
         # Time may have different names and satellite position might be missing
-        file_handler.nc = file_handler.nc.rename(
+        file_handler.nc.nc = file_handler.nc.nc.rename(
             {'time_ir_wv': 'time'}
         )
-        file_handler.nc = file_handler.nc.drop_vars(
+        file_handler.nc.nc = file_handler.nc.nc.drop_vars(
             ['sub_satellite_longitude_start']
         )
 
@@ -564,7 +564,7 @@ def test_calib_exceptions(self, file_handler):
     @pytest.mark.file_handler_data(mask_bad_quality=False)
     def test_bad_quality_warning(self, file_handler):
         """Test warning about bad VIS quality."""
-        file_handler.nc['quality_pixel_bitmask'] = 2
+        file_handler.nc.nc['quality_pixel_bitmask'] = 2
         vis = make_dataid(name='VIS', resolution=2250,
                           calibration='reflectance')
         with pytest.warns(UserWarning):

From 4bebf0b54c5ba18d301a3ec80ea969d91b83eb92 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Thu, 26 Nov 2020 20:23:28 +0000
Subject: [PATCH 46/47] Test re-assignment of dataset coordinates

---
 .../reader_tests/test_mviri_l1b_fiduceo_nc.py | 41 ++++++++++++++++++-
 1 file changed, 40 insertions(+), 1 deletion(-)

diff --git a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
index 5e8d1c29b2..11bbc84c88 100644
--- a/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
+++ b/satpy/tests/reader_tests/test_mviri_l1b_fiduceo_nc.py
@@ -29,7 +29,8 @@
 
 from satpy.readers.mviri_l1b_fiduceo_nc import (
     ALTITUDE, EQUATOR_RADIUS, POLE_RADIUS, FiduceoMviriEasyFcdrFileHandler,
-    FiduceoMviriFullFcdrFileHandler)
+    FiduceoMviriFullFcdrFileHandler, DatasetWrapper
+)
 from satpy.tests.utils import make_dataid
 
 attrs_exp = {
@@ -582,3 +583,41 @@ def test_file_pattern(self, reader):
         files = reader.select_files_from_pathnames(filenames)
         # only 3 out of 4 above should match
         assert len(files) == 3
+
+
+class TestDatasetWrapper:
+    """Unit tests for DatasetWrapper class."""
+
+    def test_reassign_coords(self):
+        """Test reassigning of coordinates.
+
+        For some reason xarray does not always assign (y, x) coordinates to
+        the high resolution datasets, although they have dimensions (y, x) and
+        coordinates y and x exist. A dataset with these properties seems
+        impossible to create (neither dropping, resetting or deleting
+        coordinates seems to work). Instead use mock as a workaround.
+        """
+        nc = mock.MagicMock(
+            coords={
+                'y': [.1, .2],
+                'x': [.3, .4]
+            },
+            dims=('y', 'x')
+        )
+        nc.__getitem__.return_value = xr.DataArray(
+            [[1, 2],
+             [3, 4]],
+            dims=('y', 'x')
+        )
+        foo_exp = xr.DataArray(
+            [[1, 2],
+             [3, 4]],
+            dims=('y', 'x'),
+            coords={
+                'y': [.1, .2],
+                'x': [.3, .4]
+            }
+        )
+        ds = DatasetWrapper(nc)
+        foo = ds['foo']
+        xr.testing.assert_equal(foo, foo_exp)

From c5b5331af64652a936a86b99a7d87459892310e0 Mon Sep 17 00:00:00 2001
From: Stephan Finkensieper <stephan.finkensieper@dwd.de>
Date: Fri, 27 Nov 2020 07:44:41 +0000
Subject: [PATCH 47/47] Skip abstract method for coverage report

---
 satpy/readers/mviri_l1b_fiduceo_nc.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/satpy/readers/mviri_l1b_fiduceo_nc.py b/satpy/readers/mviri_l1b_fiduceo_nc.py
index 6ab18cf812..c4fbcaf717 100644
--- a/satpy/readers/mviri_l1b_fiduceo_nc.py
+++ b/satpy/readers/mviri_l1b_fiduceo_nc.py
@@ -648,7 +648,7 @@ def _calibrate(self, ds, channel, calibration):
         return calib.calibrate(ds, calibration)
 
     @abc.abstractmethod
-    def _calibrate_vis(self, ds, channel, calibration):
+    def _calibrate_vis(self, ds, channel, calibration):  # pragma: no cover
         """Calibrate VIS channel. To be implemented by subclasses."""
         raise NotImplementedError