Fix generation of solar and satellite angles when lon/lats are invalid

satpy/composites/__init__.py

@@ -513,6 +513,8 @@ def get_angles(self, vis):
         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)
         sunalt, suna = get_alt_az(vis.attrs['start_time'], lons, lats)
         suna = np.rad2deg(suna)
         sunz = sun_zenith_angle(vis.attrs['start_time'], lons, lats)

satpy/composites/viirs.py

@@ -15,8 +15,7 @@
 #
 # You should have received a copy of the GNU General Public License along with
 # satpy.  If not, see <http://www.gnu.org/licenses/>.
-"""Composite classes for the VIIRS instrument.
-"""
+"""Composite classes for the VIIRS instrument."""
 
 import logging
 import os
@@ -35,9 +34,10 @@
 
 
 class VIIRSFog(CompositeBase):
+    """A simple temperature difference composite for showing fog."""
 
     def __call__(self, projectables, nonprojectables=None, **info):
-
+        """Create the temperature difference DataArray."""
         import warnings
         warnings.warn("VIIRSFog compositor is deprecated, use DifferenceCompositor "
                       "instead.", DeprecationWarning)
@@ -59,8 +59,7 @@ def __call__(self, projectables, nonprojectables=None, **info):
 
 
 class ReflectanceCorrector(CompositeBase):
-
-    """CREFL modifier
+    """Corrected Reflectance (crefl) modifier.
 
     Uses a python rewrite of the C CREFL code written for VIIRS and MODIS.
     """
@@ -87,6 +86,7 @@ def __init__(self, *args, **kwargs):
         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)
@@ -128,7 +128,7 @@ def __call__(self, datasets, optional_datasets, **info):
                                         refl_data.attrs["wavelength"],
                                         refl_data.attrs["resolution"])
         use_abi = vis.attrs['sensor'] == 'abi'
-        lons, lats = vis.attrs['area'].get_lonlats_dask(chunks=vis.chunks)
+        lons, lats = vis.attrs['area'].get_lonlats(chunks=vis.chunks)
         results = run_crefl(refl_data,
                             coefficients,
                             lons,
@@ -149,11 +149,12 @@ def __call__(self, datasets, optional_datasets, **info):
         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_dask(
-            chunks=vis.data.chunks)
+        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)
@@ -395,6 +396,7 @@ class ERFDNB(CompositeBase):
     """
 
     def __init__(self, *args, **kwargs):
+        """Initialize ERFDNB specific keyword arguments."""
         self.saturation_correction = kwargs.pop("saturation_correction",
                                                 False)
         super(ERFDNB, self).__init__(*args, **kwargs)
@@ -417,6 +419,7 @@ def _saturation_correction(self, dnb_data, unit_factor, min_val,
         return np.sqrt(inner_sqrt)
 
     def __call__(self, datasets, **info):
+        """Create the composite DataArray object for ERFDNB."""
         if len(datasets) != 4:
             raise ValueError("Expected 4 datasets, got %d" % (len(datasets), ))
 
@@ -492,17 +495,17 @@ def make_day_night_masks(solarZenithAngle,
                          highAngleCutoff,
                          lowAngleCutoff,
                          stepsDegrees=None):
-    """
-    given information on the solarZenithAngle for each point,
-    generate masks defining where the day, night, and mixed regions are
+    """Generate masks for day, night, and twilight regions.
+
+    Masks are created from the provided solar zenith angle data.
 
-    optionally provide the highAngleCutoff and lowAngleCutoff that define
+    Optionally provide the highAngleCutoff and lowAngleCutoff that define
     the limits of the terminator region (if no cutoffs are given the
-    DEFAULT_HIGH_ANGLE and DEFAULT_LOW_ANGLE will be used)
+    DEFAULT_HIGH_ANGLE and DEFAULT_LOW_ANGLE will be used).
 
-    optionally provide the stepsDegrees that define how many degrees each
+    Optionally provide the stepsDegrees that define how many degrees each
     "mixed" mask in the terminator region should be (if no stepsDegrees is
-    given, the whole terminator region will be one mask)
+    given, the whole terminator region will be one mask).
     """
     # if the caller passes None, we're only doing one step
     stepsDegrees = highAngleCutoff - lowAngleCutoff if stepsDegrees is None else stepsDegrees
@@ -545,8 +548,9 @@ def histogram_equalization(
         log_offset=None,
         local_radius_px=None,
         out=None):
-    """
-    Perform a histogram equalization on the data selected by mask_to_equalize.
+    """Perform a histogram equalization on the data.
+
+    Data is selected by the mask_to_equalize mask.
     The data will be separated into number_of_bins levels for equalization and
     outliers beyond +/- std_mult_cutoff*std will be ignored.
 
@@ -556,7 +560,6 @@ def histogram_equalization(
 
     Note: the data will be changed in place.
     """
-
     out = out if out is not None else data.copy()
     mask_to_use = mask_to_equalize if valid_data_mask is None else valid_data_mask
 
@@ -603,22 +606,18 @@ def local_histogram_equalization(data, mask_to_equalize, valid_data_mask=None, n
                                  ):
     """Equalize the provided data (in the mask_to_equalize) using adaptive histogram equalization.
 
-    tiles of width/height (2 * local_radius_px + 1) will be calculated and results for each pixel will be bilinerarly
+    Tiles of width/height (2 * local_radius_px + 1) will be calculated and results for each pixel will be bilinearly
     interpolated from the nearest 4 tiles when pixels fall near the edge of the image (there is no adjacent tile) the
-    resultant interpolated sum from the available tiles will be multipled to account for the weight of any missing
+    resultant interpolated sum from the available tiles will be multiplied to account for the weight of any missing
     tiles::
 
         pixel total interpolated value = pixel available interpolated value / (1 - missing interpolation weight)
 
-    if ``do_zerotoone_normalization`` is True the data will be scaled so that all data in the mask_to_equalize falls
-    between 0 and 1; otherwise the data in mask_to_equalize will all fall between 0 and number_of_bins
-
-    Returns:
-
-        The equalized data
+    If ``do_zerotoone_normalization`` is True the data will be scaled so that all data in the mask_to_equalize falls
+    between 0 and 1; otherwise the data in mask_to_equalize will all fall between 0 and number_of_bins.
 
+    Returns: The equalized data
     """
-
     out = out if out is not None else np.zeros_like(data)
     # if we don't have a valid mask, use the mask of what we should be
     # equalizing
@@ -823,7 +822,6 @@ def _histogram_equalization_helper(valid_data, number_of_bins, clip_limit=None,
         cumulative distribution function and bin information
 
     """
-
     # bucket all the selected data using np's histogram function
     temp_histogram, temp_bins = np.histogram(valid_data, number_of_bins)
 
@@ -873,17 +871,17 @@ def _histogram_equalization_helper(valid_data, number_of_bins, clip_limit=None,
 
 
 def _calculate_weights(tile_size):
-    """
-    calculate a weight array that will be used to quickly bilinearly-interpolate the histogram equalizations
-    tile size should be the width and height of a tile in pixels
+    """Calculate a weight array for bilinear interpolation of histogram tiles.
 
-    returns a 4D weight array, where the first 2 dimensions correspond to the grid of where the tiles are
-    relative to the tile being interpolated
-    """
+    The weight array will be used to quickly bilinearly-interpolate the
+    histogram equalizations tile size should be the width and height of a tile
+    in pixels.
 
+    Returns: 4D weight array where the first 2 dimensions correspond to the
+        grid of where the tiles are relative to the tile being interpolated.
+    """
     # we are essentially making a set of weight masks for an ideal center tile
     # that has all 8 surrounding tiles available
-
     # create our empty template tiles
     template_tile = np.zeros((3, 3, tile_size, tile_size), dtype=np.float32)
     """
@@ -998,7 +996,6 @@ def _linear_normalization_from_0to1(
     the correct theoretical current max and min so it can scale the data
     accordingly.
     """
-
     LOG.debug(message)
     if theoretical_min != 0:
         data[mask] = data[mask] - theoretical_min
@@ -1010,14 +1007,16 @@ class NCCZinke(CompositeBase):
     """Equalized DNB composite using the Zinke algorithm [#ncc1]_.
 
     References:
-
         .. [#ncc1] Stephan Zinke (2017),
-               A simplified high and near-constant contrast approach for the display of VIIRS day/night band imagery
+
+               A simplified high and near-constant contrast approach for the
+               display of VIIRS day/night band imagery
                :doi:`10.1080/01431161.2017.1338838`
 
     """
 
     def __call__(self, datasets, **info):
+        """Create HNCC DNB composite."""
         if len(datasets) != 4:
             raise ValueError("Expected 4 datasets, got %d" % (len(datasets),))
 
@@ -1060,8 +1059,8 @@ def __call__(self, datasets, **info):
         return dnb_data
 
     def gain_factor(self, theta):
-        return theta.map_blocks(self._gain_factor,
-                                dtype=theta.dtype)
+        """Compute gain factor in a dask-friendly manner."""
+        return theta.map_blocks(self._gain_factor, dtype=theta.dtype)
 
     @staticmethod
     def _gain_factor(theta):

satpy/tests/compositor_tests/__init__.py

@@ -1067,6 +1067,8 @@ def test_multiple_sensors(self):
 
 
 class TestPSPAtmosphericalCorrection(unittest.TestCase):
+    """Test the pyspectral-based atmospheric correction modifier."""
+
     def setUp(self):
         """Patch in-class imports."""
         self.orbital = mock.MagicMock()
@@ -1109,6 +1111,8 @@ def test_call(self, get_satpos, *mocks):
 
 
 class TestPSPRayleighReflectance(unittest.TestCase):
+    """Test the pyspectral-based rayleigh correction modifier."""
+
     def setUp(self):
         """Patch in-class imports."""
         self.astronomy = mock.MagicMock()
@@ -1134,7 +1138,13 @@ def test_get_angles(self, get_satpos):
         self.orbital.get_observer_look.return_value = 0, 0
         self.astronomy.get_alt_az.return_value = 0, 0
         area = mock.MagicMock()
-        area.get_lonlats.return_value = 'lons', 'lats'
+        lons = np.zeros((5, 5))
+        lons[1, 1] = np.inf
+        lons = da.from_array(lons, chunks=5)
+        lats = np.zeros((5, 5))
+        lats[1, 1] = np.inf
+        lats = da.from_array(lats, chunks=5)
+        area.get_lonlats.return_value = (lons, lats)
         vis = mock.MagicMock(attrs={'area': area,
                                     'start_time': 'start_time'})
 
@@ -1144,8 +1154,12 @@ def test_get_angles(self, get_satpos):
 
         # Check arguments of get_orbserver_look() call, especially the altitude
         # unit conversion from meters to kilometers
-        self.orbital.get_observer_look.assert_called_with(
-            'sat_lon', 'sat_lat', 12345.678, 'start_time', 'lons', 'lats', 0)
+        self.orbital.get_observer_look.assert_called_once()
+        args = self.orbital.get_observer_look.call_args[0]
+        self.assertEqual(args[:4], ('sat_lon', 'sat_lat', 12345.678, 'start_time'))
+        self.assertIsInstance(args[4], da.Array)
+        self.assertIsInstance(args[5], da.Array)
+        self.assertEqual(args[6], 0)
 
 
 def suite():

satpy/tests/compositor_tests/test_viirs.py

@@ -17,12 +17,7 @@
 # satpy.  If not, see <http://www.gnu.org/licenses/>.
 """Tests for VIIRS compositors."""
 
-import sys
-
-if sys.version_info < (2, 7):
-    import unittest2 as unittest
-else:
-    import unittest
+import unittest
 try:
     from unittest import mock
 except ImportError:
@@ -283,6 +278,7 @@ def test_reflectance_corrector_abi(self):
             DatasetID(name='solar_zenith_angle')])
 
         area, dnb = self.data_area_ref_corrector()
+        print(dnb.compute())
         c01 = xr.DataArray(dnb,
                            dims=('y', 'x'),
                            attrs={'satellite_longitude': -89.5, 'satellite_latitude': 0.0,
@@ -311,21 +307,21 @@ def test_reflectance_corrector_abi(self):
         self.assertEqual(res.attrs['area'], area)
         self.assertEqual(res.attrs['ancillary_variables'], [])
         data = res.values
-        self.assertLess(abs(np.mean(data) - 29.907390988422513), 1e-10)
+        self.assertLess(abs(np.nanmean(data) - 26.00760944144745), 1e-10)
         self.assertEqual(data.shape, (5, 10))
-        unique = np.unique(data)
+        unique = np.unique(data[~np.isnan(data)])
         np.testing.assert_allclose(unique, [-1.0, 4.210745457958135, 6.7833906076177595, 8.730371329824473,
                                             10.286627569545209, 11.744159436709374, 12.20226097829902,
                                             13.501444598985305, 15.344399223932212, 17.173329483996515,
                                             17.28798660754271, 18.29594550575925, 19.076835059905125,
                                             19.288331720959864, 19.77043407084455, 19.887082168377006,
                                             20.091028778326375, 20.230341149334617, 20.457671064690196,
                                             20.82686905639114, 21.021094816441195, 21.129963777952124,
-                                            21.94957397026227, 41.601857910095575, 43.963919057675504,
+                                            41.601857910095575, 43.963919057675504,
                                             46.21672174361075, 46.972099490462085, 47.497072794632835,
                                             47.80393007974336, 47.956765988770385, 48.043025685032106,
                                             51.909142813383916, 58.8234273736508, 68.84706145641482, 69.91085190887961,
-                                            71.10179768327806, 71.33161009169649, 78.81291424983952])
+                                            71.10179768327806, 71.33161009169649])
 
     def test_reflectance_corrector_viirs(self):
         """Test ReflectanceCorrector modifier with VIIRS data."""
@@ -479,6 +475,8 @@ def make_xarray(self, name, calibration, wavelength=None, modifiers=None, resolu
 
 
 class ViirsReflectanceCorrectorTest(unittest.TestCase):
+    """Tests for the VIIRS/MODIS Corrected Reflectance modifier."""
+
     def setUp(self):
         """Patch in-class imports."""
         self.astronomy = mock.MagicMock()
@@ -497,14 +495,22 @@ def tearDown(self):
     @mock.patch('satpy.composites.viirs.get_satpos')
     def test_get_angles(self, get_satpos):
         """Test sun and satellite angle calculation."""
+        import numpy as np
+        import dask.array as da
         from satpy.composites.viirs import ReflectanceCorrector
 
         # Patch methods
         get_satpos.return_value = 'sat_lon', 'sat_lat', 12345678
         self.orbital.get_observer_look.return_value = 0, 0
         self.astronomy.get_alt_az.return_value = 0, 0
         area = mock.MagicMock()
-        area.get_lonlats_dask.return_value = 'lons', 'lats'
+        lons = np.zeros((5, 5))
+        lons[1, 1] = np.inf
+        lons = da.from_array(lons, chunks=5)
+        lats = np.zeros((5, 5))
+        lats[1, 1] = np.inf
+        lats = da.from_array(lats, chunks=5)
+        area.get_lonlats.return_value = (lons, lats)
         vis = mock.MagicMock(attrs={'area': area,
                                     'start_time': 'start_time'})
 
@@ -514,8 +520,12 @@ def test_get_angles(self, get_satpos):
 
         # Check arguments of get_orbserver_look() call, especially the altitude
         # unit conversion from meters to kilometers
-        self.orbital.get_observer_look.assert_called_with(
-            'sat_lon', 'sat_lat', 12345.678, 'start_time', 'lons', 'lats', 0)
+        self.orbital.get_observer_look.assert_called_once()
+        args = self.orbital.get_observer_look.call_args[0]
+        self.assertEqual(args[:4], ('sat_lon', 'sat_lat', 12345.678, 'start_time'))
+        self.assertIsInstance(args[4], da.Array)
+        self.assertIsInstance(args[5], da.Array)
+        self.assertEqual(args[6], 0)
 
 
 def suite():