Update VIIRS/MODIS ReflecanceCorrector modifier to download DEM data at runtime

satpy/composites/viirs.py

@@ -18,133 +18,18 @@
 """Composite classes for the VIIRS instrument."""
 
 import logging
-import os
 
 import numpy as np
 import dask
 import dask.array as da
 import xarray as xr
 
-import satpy
 from satpy.composites import CompositeBase, GenericCompositor
 from satpy.dataset import combine_metadata
-from satpy.utils import get_satpos
 
 LOG = logging.getLogger(__name__)
 
 
-class ReflectanceCorrector(CompositeBase):
-    """Corrected Reflectance (crefl) modifier.
-
-    Uses a python rewrite of the C CREFL code written for VIIRS and MODIS.
-    """
-
-    def __init__(self, *args, **kwargs):
-        """Initialize the compositor with values from the user or from the configuration file.
-
-        If `dem_filename` can't be found or opened then correction is done
-        assuming TOA or sealevel options.
-
-        :param dem_filename: path to the ancillary 'averaged heights' file
-                             default: CMGDEM.hdf
-                             environment override: os.path.join(<SATPY_ANCPATH>, <CREFL_ANCFILENAME>)
-        :param dem_sds: variable name to load from the ancillary file
-        """
-        dem_filename = kwargs.pop("dem_filename",
-                                  os.environ.get("CREFL_ANCFILENAME",
-                                                 "CMGDEM.hdf"))
-        if os.path.exists(dem_filename):
-            self.dem_file = dem_filename
-        else:
-            self.dem_file = os.path.join(satpy.config.get('data_dir'),
-                                         dem_filename)
-        self.dem_sds = kwargs.pop("dem_sds", "averaged elevation")
-        super(ReflectanceCorrector, self).__init__(*args, **kwargs)
-
-    def __call__(self, datasets, optional_datasets, **info):
-        """Create modified DataArray object by applying the crefl algorithm."""
-        if not optional_datasets or len(optional_datasets) != 4:
-            vis = self.match_data_arrays([datasets[0]])[0]
-            sensor_aa, sensor_za, solar_aa, solar_za = self.get_angles(vis)
-        else:
-            vis, sensor_aa, sensor_za, solar_aa, solar_za = self.match_data_arrays(
-                datasets + optional_datasets)
-            # get the dask array underneath
-            sensor_aa = sensor_aa.data
-            sensor_za = sensor_za.data
-            solar_aa = solar_aa.data
-            solar_za = solar_za.data
-        # angles must be xarrays
-        sensor_aa = xr.DataArray(sensor_aa, dims=['y', 'x'])
-        sensor_za = xr.DataArray(sensor_za, dims=['y', 'x'])
-        solar_aa = xr.DataArray(solar_aa, dims=['y', 'x'])
-        solar_za = xr.DataArray(solar_za, dims=['y', 'x'])
-        refl_data = datasets[0]
-        if refl_data.attrs.get("rayleigh_corrected"):
-            return refl_data
-        if os.path.isfile(self.dem_file):
-            LOG.debug("Loading CREFL averaged elevation information from: %s",
-                      self.dem_file)
-            from netCDF4 import Dataset as NCDataset
-            # HDF4 file, NetCDF library needs to be compiled with HDF4 support
-            nc = NCDataset(self.dem_file, "r")
-            # average elevation is stored as a 16-bit signed integer but with
-            # scale factor 1 and offset 0, convert it to float here
-            avg_elevation = nc.variables[self.dem_sds][:].astype(np.float64)
-            if isinstance(avg_elevation, np.ma.MaskedArray):
-                avg_elevation = avg_elevation.filled(np.nan)
-        else:
-            avg_elevation = None
-
-        from satpy.composites.crefl_utils import run_crefl, get_coefficients
-
-        percent = refl_data.attrs["units"] == "%"
-
-        coefficients = get_coefficients(refl_data.attrs["sensor"],
-                                        refl_data.attrs["wavelength"],
-                                        refl_data.attrs["resolution"])
-        use_abi = vis.attrs['sensor'] == 'abi'
-        lons, lats = vis.attrs['area'].get_lonlats(chunks=vis.chunks)
-        results = run_crefl(refl_data,
-                            coefficients,
-                            lons,
-                            lats,
-                            sensor_aa,
-                            sensor_za,
-                            solar_aa,
-                            solar_za,
-                            avg_elevation=avg_elevation,
-                            percent=percent,
-                            use_abi=use_abi)
-        info.update(refl_data.attrs)
-        info["rayleigh_corrected"] = True
-        factor = 100. if percent else 1.
-        results = results * factor
-        results.attrs = info
-        self.apply_modifier_info(refl_data, results)
-        return results
-
-    def get_angles(self, vis):
-        """Get sun and satellite angles to use in crefl calculations."""
-        from pyorbital.astronomy import get_alt_az, sun_zenith_angle
-        from pyorbital.orbital import get_observer_look
-        lons, lats = vis.attrs['area'].get_lonlats(chunks=vis.data.chunks)
-        lons = da.where(lons >= 1e30, np.nan, lons)
-        lats = da.where(lats >= 1e30, np.nan, lats)
-        suna = get_alt_az(vis.attrs['start_time'], lons, lats)[1]
-        suna = np.rad2deg(suna)
-        sunz = sun_zenith_angle(vis.attrs['start_time'], lons, lats)
-        sat_lon, sat_lat, sat_alt = get_satpos(vis)
-        sata, satel = get_observer_look(
-            sat_lon,
-            sat_lat,
-            sat_alt / 1000.0,  # km
-            vis.attrs['start_time'],
-            lons, lats, 0)
-        satz = 90 - satel
-        return sata, satz, suna, sunz
-
-
 class HistogramDNB(CompositeBase):
     """Histogram equalized DNB composite.
 

satpy/etc/composites/abi.yaml

@@ -2,7 +2,7 @@ sensor_name: visir/abi
 
 modifiers:
   rayleigh_corrected_crefl:
-    modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector
+    modifier: !!python/name:satpy.modifiers.atmosphere.ReflectanceCorrector
     dem_filename: CMGDEM.hdf
     optional_prerequisites:
       - name: satellite_azimuth_angle

satpy/etc/composites/modis.yaml

@@ -1,9 +1,10 @@
 sensor_name: visir/modis
 modifiers:
   rayleigh_corrected_crefl:
-    modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector
-    dem_filename: CMGDEM.hdf
-    optional_prerequisites:
+    modifier: !!python/name:satpy.modifiers.atmosphere.ReflectanceCorrector
+    url: "https://www.ssec.wisc.edu/~davidh/polar2grid/modis_crefl/tbase.hdf"
+    known_hash: "sha256:ed5183cddce905361c1cac8ae6e3a447212875ea421a05747751efe76f8a068e"
+    prerequisites:
     - name: satellite_azimuth_angle
     - name: satellite_zenith_angle
     - name: solar_azimuth_angle

satpy/etc/composites/viirs.yaml

@@ -2,9 +2,10 @@ sensor_name: visir/viirs
 
 modifiers:
   rayleigh_corrected_crefl:
-    modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector
-    dem_filename: CMGDEM.hdf
-    optional_prerequisites:
+    modifier: !!python/name:satpy.modifiers.atmosphere.ReflectanceCorrector
+    url: "https://www.ssec.wisc.edu/~davidh/polar2grid/viirs_crefl/CMGDEM.hdf"
+    known_hash: "sha256:f33f1f867d79fff4fafe128f61c154236dd74fcc97bf418ea1437977a38d0604"
+    prerequisites:
     - name: satellite_azimuth_angle
       resolution: 742
     - name: satellite_zenith_angle
@@ -15,9 +16,10 @@ modifiers:
       resolution: 742
 
   rayleigh_corrected_crefl_iband:
-    modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector
-    dem_filename: CMGDEM.hdf
-    optional_prerequisites:
+    modifier: !!python/name:satpy.modifiers.atmosphere.ReflectanceCorrector
+    url: "https://www.ssec.wisc.edu/~davidh/polar2grid/viirs_crefl/CMGDEM.hdf"
+    known_hash: "sha256:f33f1f867d79fff4fafe128f61c154236dd74fcc97bf418ea1437977a38d0604"
+    prerequisites:
     - name: satellite_azimuth_angle
       resolution: 371
     - name: satellite_zenith_angle

satpy/modifiers/_crefl.py

@@ -0,0 +1,164 @@
+#!/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/>.
+"""Classes related to the CREFL (corrected reflectance) modifier."""
+
+import logging
+import warnings
+
+import numpy as np
+import xarray as xr
+from dask import array as da
+from satpy.aux_download import DataDownloadMixin, retrieve
+from satpy.modifiers import ModifierBase
+from satpy.utils import get_satpos
+
+LOG = logging.getLogger(__name__)
+
+
+class ReflectanceCorrector(ModifierBase, DataDownloadMixin):
+    """Corrected Reflectance (crefl) modifier.
+
+    Uses a python rewrite of the C CREFL code written for VIIRS and MODIS.
+    """
+
+    def __init__(self, *args, **kwargs):
+        """Initialize the compositor with values from the user or from the configuration file.
+
+        If `dem_filename` can't be found or opened then correction is done
+        assuming TOA or sealevel options.
+
+        Args:
+            dem_filename (str): DEPRECATED
+            url (str): URL or local path to the Digital Elevation Model (DEM)
+                HDF4 file. If unset (None or empty string), then elevation
+                is assumed to be 0 everywhere.
+            known_hash (str): Optional SHA256 checksum to verify the download
+                of ``url``.
+            dem_sds (str): Name of the variable in the elevation file to load.
+
+        """
+        dem_filename = kwargs.pop("dem_filename", None)
+        if dem_filename is not None:
+            warnings.warn("'dem_filename' for 'ReflectanceCorrector' is "
+                          "deprecated. Use 'url' instead.", DeprecationWarning)
+
+        self.dem_sds = kwargs.pop("dem_sds", "averaged elevation")
+        self.url = kwargs.pop('url', None)
+        self.known_hash = kwargs.pop('known_hash', None)
+        super(ReflectanceCorrector, self).__init__(*args, **kwargs)
+        self.dem_cache_key = None
+        if self.url:
+            reg_files = self.register_data_files([{
+                'url': self.url, 'known_hash': self.known_hash}
+            ])
+            self.dem_cache_key = reg_files[0]
+
+    def __call__(self, datasets, optional_datasets, **info):
+        """Create modified DataArray object by applying the crefl algorithm."""
+        from satpy.composites.crefl_utils import run_crefl, get_coefficients
+
+        datasets = self._get_data_and_angles(datasets, optional_datasets)
+        refl_data, sensor_aa, sensor_za, solar_aa, solar_za = datasets
+        avg_elevation = self._get_average_elevation()
+        is_percent = refl_data.attrs["units"] == "%"
+        coefficients = get_coefficients(refl_data.attrs["sensor"],
+                                        refl_data.attrs["wavelength"],
+                                        refl_data.attrs["resolution"])
+        use_abi = refl_data.attrs['sensor'] == 'abi'
+        lons, lats = refl_data.attrs['area'].get_lonlats(chunks=refl_data.chunks)
+        results = run_crefl(refl_data,
+                            coefficients,
+                            lons,
+                            lats,
+                            sensor_aa,
+                            sensor_za,
+                            solar_aa,
+                            solar_za,
+                            avg_elevation=avg_elevation,
+                            percent=is_percent,
+                            use_abi=use_abi)
+        info.update(refl_data.attrs)
+        info["rayleigh_corrected"] = True
+        factor = 100. if is_percent else 1.
+        results = results * factor
+        results.attrs = info
+        self.apply_modifier_info(refl_data, results)
+        return results
+
+    def _get_average_elevation(self):
+        if self.dem_cache_key is None:
+            return
+
+        LOG.debug("Loading CREFL averaged elevation information from: %s",
+                  self.dem_cache_key)
+        local_filename = retrieve(self.dem_cache_key)
+        from netCDF4 import Dataset as NCDataset
+        # HDF4 file, NetCDF library needs to be compiled with HDF4 support
+        nc = NCDataset(local_filename, "r")
+        # average elevation is stored as a 16-bit signed integer but with
+        # scale factor 1 and offset 0, convert it to float here
+        avg_elevation = nc.variables[self.dem_sds][:].astype(np.float64)
+        if isinstance(avg_elevation, np.ma.MaskedArray):
+            avg_elevation = avg_elevation.filled(np.nan)
+        return avg_elevation
+
+    def _get_data_and_angles(self, datasets, optional_datasets):
+        all_datasets = datasets + optional_datasets
+        if len(all_datasets) == 1:
+            vis = self.match_data_arrays(datasets)[0]
+            sensor_aa, sensor_za, solar_aa, solar_za = self.get_angles(vis)
+        elif len(all_datasets) == 5:
+            vis, sensor_aa, sensor_za, solar_aa, solar_za = self.match_data_arrays(
+                datasets + optional_datasets)
+            # get the dask array underneath
+            sensor_aa = sensor_aa.data
+            sensor_za = sensor_za.data
+            solar_aa = solar_aa.data
+            solar_za = solar_za.data
+        else:
+            raise ValueError("Not sure how to handle provided dependencies. "
+                             "Either all 4 angles must be provided or none of "
+                             "of them.")
+
+        # angles must be xarrays
+        sensor_aa = xr.DataArray(sensor_aa, dims=['y', 'x'])
+        sensor_za = xr.DataArray(sensor_za, dims=['y', 'x'])
+        solar_aa = xr.DataArray(solar_aa, dims=['y', 'x'])
+        solar_za = xr.DataArray(solar_za, dims=['y', 'x'])
+        refl_data = datasets[0]
+        return refl_data, sensor_aa, sensor_za, solar_aa, solar_za
+
+    def get_angles(self, vis):
+        """Get sun and satellite angles to use in crefl calculations."""
+        from pyorbital.astronomy import get_alt_az, sun_zenith_angle
+        from pyorbital.orbital import get_observer_look
+        lons, lats = vis.attrs['area'].get_lonlats(chunks=vis.data.chunks)
+        lons = da.where(lons >= 1e30, np.nan, lons)
+        lats = da.where(lats >= 1e30, np.nan, lats)
+        suna = get_alt_az(vis.attrs['start_time'], lons, lats)[1]
+        suna = np.rad2deg(suna)
+        sunz = sun_zenith_angle(vis.attrs['start_time'], lons, lats)
+        sat_lon, sat_lat, sat_alt = get_satpos(vis)
+        sata, satel = get_observer_look(
+            sat_lon,
+            sat_lat,
+            sat_alt / 1000.0,  # km
+            vis.attrs['start_time'],
+            lons, lats, 0)
+        satz = 90 - satel
+        return sata, satz, suna, sunz

satpy/modifiers/atmosphere.py

@@ -25,6 +25,7 @@
 import xarray as xr
 
 from satpy.modifiers import ModifierBase
+from satpy.modifiers._crefl import ReflectanceCorrector  # noqa
 from satpy.utils import get_satpos
 
 logger = logging.getLogger(__name__)