diff --git a/satpy/composites/__init__.py b/satpy/composites/__init__.py index 2bdcabe27e..ced0d00b8d 100644 --- a/satpy/composites/__init__.py +++ b/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) diff --git a/satpy/composites/viirs.py b/satpy/composites/viirs.py index c1527ca38a..beb11b4688 100644 --- a/satpy/composites/viirs.py +++ b/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 . -"""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): diff --git a/satpy/tests/compositor_tests/__init__.py b/satpy/tests/compositor_tests/__init__.py index 042377fcaf..5c0f72986f 100644 --- a/satpy/tests/compositor_tests/__init__.py +++ b/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(): diff --git a/satpy/tests/compositor_tests/test_viirs.py b/satpy/tests/compositor_tests/test_viirs.py index 7a1a468faf..44f5644398 100644 --- a/satpy/tests/compositor_tests/test_viirs.py +++ b/satpy/tests/compositor_tests/test_viirs.py @@ -17,12 +17,7 @@ # satpy. If not, see . """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,9 +307,9 @@ 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, @@ -321,11 +317,11 @@ def test_reflectance_corrector_abi(self): 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,6 +495,8 @@ 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 @@ -504,7 +504,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_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():