diff --git a/satpy/composites/__init__.py b/satpy/composites/__init__.py index fd50cced99..8e50121c88 100644 --- a/satpy/composites/__init__.py +++ b/satpy/composites/__init__.py @@ -313,6 +313,18 @@ class SunZenithCorrectorBase(CompositeBase): coszen = WeakValueDictionary() + def __init__(self, max_sza=95.0, **kwargs): + """Collect custom configuration values. + + Args: + max_sza (float): Maximum solar zenith angle in degrees that is + considered valid and correctable. Default 95.0. + + """ + self.max_sza = max_sza + self.max_sza_cos = np.cos(np.deg2rad(max_sza)) if max_sza is not None else None + super(SunZenithCorrectorBase, self).__init__(**kwargs) + def __call__(self, projectables, **info): projectables = self.check_areas(projectables) vis = projectables[0] @@ -320,10 +332,7 @@ def __call__(self, projectables, **info): LOG.debug("Sun zen correction already applied") return vis - if hasattr(vis.attrs["area"], 'name'): - area_name = vis.attrs["area"].name - else: - area_name = 'swath' + str(vis.shape) + area_name = hash(vis.attrs['area']) key = (vis.attrs["start_time"], area_name) tic = time.time() LOG.debug("Applying sun zen correction") @@ -334,22 +343,19 @@ def __call__(self, projectables, **info): LOG.debug("Computing sun zenith angles.") lons, lats = vis.attrs["area"].get_lonlats_dask(CHUNK_SIZE) - coszen = xr.DataArray(cos_zen(vis.attrs["start_time"], - lons, lats), - dims=['y', 'x'], - coords=[vis['y'], vis['x']]) - coszen = coszen.where((coszen > 0.035) & (coszen < 1)) + coszen = xr.DataArray(cos_zen(vis.attrs["start_time"], lons, lats), + dims=['y', 'x'], coords=[vis['y'], vis['x']]) + if self.max_sza is not None: + coszen = coszen.where(coszen >= self.max_sza_cos) self.coszen[key] = coszen else: - coszen = xu.cos(xu.deg2rad(projectables[1])) + coszen = np.cos(np.deg2rad(projectables[1])) self.coszen[key] = coszen proj = self._apply_correction(vis, coszen) proj.attrs = vis.attrs.copy() self.apply_modifier_info(vis, proj) - LOG.debug( - "Sun-zenith correction applied. Computation time: %5.1f (sec)", - time.time() - tic) + LOG.debug("Sun-zenith correction applied. Computation time: %5.1f (sec)", time.time() - tic) return proj def _apply_correction(self, proj, coszen): @@ -357,11 +363,44 @@ def _apply_correction(self, proj, coszen): class SunZenithCorrector(SunZenithCorrectorBase): - """Standard sun zenith correction, 1/cos(sunz).""" + """Standard sun zenith correction using ``1 / cos(sunz)``. + + In addition to adjusting the provided reflectances by the cosine of the + solar zenith angle, this modifier forces all reflectances beyond a + solar zenith angle of ``max_sza`` to 0. It also gradually reduces the + amount of correction done between ``correction_limit`` and ``max_sza``. If + ``max_sza`` is ``None`` then a constant correction is applied to zenith + angles beyond ``correction_limit``. + + To set ``max_sza`` to ``None`` in a YAML configuration file use: + + .. code-block:: yaml + + sunz_corrected: + compositor: !!python/name:satpy.composites.SunZenithCorrector + max_sza: !!null + optional_prerequisites: + - solar_zenith_angle + + """ + + def __init__(self, correction_limit=88., **kwargs): + """Collect custom configuration values. + + Args: + correction_limit (float): Maximum solar zenith angle to apply the + correction in degrees. Pixels beyond this limit have a + constant correction applied. Default 88. + max_sza (float): Maximum solar zenith angle in degrees that is + considered valid and correctable. Default 95.0. + + """ + self.correction_limit = correction_limit + super(SunZenithCorrector, self).__init__(**kwargs) def _apply_correction(self, proj, coszen): LOG.debug("Apply the standard sun-zenith correction [1/cos(sunz)]") - return sunzen_corr_cos(proj, coszen) + return sunzen_corr_cos(proj, coszen, limit=self.correction_limit, max_sza=self.max_sza) class EffectiveSolarPathLengthCorrector(SunZenithCorrectorBase): @@ -369,11 +408,43 @@ class EffectiveSolarPathLengthCorrector(SunZenithCorrectorBase): (2006): https://doi.org/10.1175/JAS3682.1 + In addition to adjusting the provided reflectances by the cosine of the + solar zenith angle, this modifier forces all reflectances beyond a + solar zenith angle of `max_sza` to 0 to reduce noise in the final data. + It also gradually reduces the amount of correction done between + ``correction_limit`` and ``max_sza``. If ``max_sza`` is ``None`` then a + constant correction is applied to zenith angles beyond + ``correction_limit``. + + To set ``max_sza`` to ``None`` in a YAML configuration file use: + + .. code-block:: yaml + + effective_solar_pathlength_corrected: + compositor: !!python/name:satpy.composites.EffectiveSolarPathLengthCorrector + max_sza: !!null + optional_prerequisites: + - solar_zenith_angle + """ + def __init__(self, correction_limit=88., **kwargs): + """Collect custom configuration values. + + Args: + correction_limit (float): Maximum solar zenith angle to apply the + correction in degrees. Pixels beyond this limit have a + constant correction applied. Default 88. + max_sza (float): Maximum solar zenith angle in degrees that is + considered valid and correctable. Default 95.0. + + """ + self.correction_limit = correction_limit + super(EffectiveSolarPathLengthCorrector, self).__init__(**kwargs) + def _apply_correction(self, proj, coszen): LOG.debug("Apply the effective solar atmospheric path length correction method by Li and Shibata") - return atmospheric_path_length_correction(proj, coszen) + return atmospheric_path_length_correction(proj, coszen, limit=self.correction_limit, max_sza=self.max_sza) class PSPRayleighReflectance(CompositeBase): @@ -829,7 +900,7 @@ def __call__(self, projectables, **kwargs): lim_low = np.cos(np.deg2rad(self.lim_low)) lim_high = np.cos(np.deg2rad(self.lim_high)) try: - coszen = xu.cos(xu.deg2rad(projectables[2])) + coszen = np.cos(np.deg2rad(projectables[2])) except IndexError: from pyorbital.astronomy import cos_zen LOG.debug("Computing sun zenith angles.") @@ -1163,7 +1234,7 @@ def __call__(self, datasets, optional_datasets=None, **info): r, g, b = datasets[:3] # combine the masks - mask = ~(da.isnull(r.data) | da.isnull(g.data) | da.isnull(b.data)) + mask = ~(r.isnull() | g.isnull() | b.isnull()) r = r.where(mask) g = g.where(mask) b = b.where(mask) @@ -1202,8 +1273,6 @@ def four_element_average_dask(d): """Average every 4 elements (2x2) in a 2D array""" def _mean4(data, offset=(0, 0), block_id=None): rows, cols = data.shape - rows2, cols2 = data.shape - pad = [] # we assume that the chunks except the first ones are aligned if block_id[0] == 0: row_offset = offset[0] % 2 diff --git a/satpy/composites/viirs.py b/satpy/composites/viirs.py index 0a163dd681..3bb233a7a9 100644 --- a/satpy/composites/viirs.py +++ b/satpy/composites/viirs.py @@ -160,7 +160,7 @@ def get_angles(self, vis): lons, lats = vis.attrs['area'].get_lonlats_dask( chunks=vis.data.chunks) suna = get_alt_az(vis.attrs['start_time'], lons, lats)[1] - suna = xu.rad2deg(suna) + suna = np.rad2deg(suna) sunz = sun_zenith_angle(vis.attrs['start_time'], lons, lats) sata, satel = get_observer_look( vis.attrs['satellite_longitude'], diff --git a/satpy/etc/composites/visir.yaml b/satpy/etc/composites/visir.yaml index 237b8cfb3b..b7608986f4 100644 --- a/satpy/etc/composites/visir.yaml +++ b/satpy/etc/composites/visir.yaml @@ -4,10 +4,12 @@ modifiers: sunz_corrected: compositor: !!python/name:satpy.composites.SunZenithCorrector optional_prerequisites: - - satellite_zenith_angle + - solar_zenith_angle effective_solar_pathlength_corrected: compositor: !!python/name:satpy.composites.EffectiveSolarPathLengthCorrector + optional_prerequisites: + - solar_zenith_angle co2_corrected: compositor: !!python/name:satpy.composites.CO2Corrector diff --git a/satpy/tests/compositor_tests/__init__.py b/satpy/tests/compositor_tests/__init__.py index b509824a68..b00dbe71b1 100644 --- a/satpy/tests/compositor_tests/__init__.py +++ b/satpy/tests/compositor_tests/__init__.py @@ -223,6 +223,58 @@ def test_self_sharpened_basic(self): np.testing.assert_allclose(res[2], np.array([[4, 4], [4, 4]], dtype=np.float64)) +class TestSunZenithCorrector(unittest.TestCase): + def setUp(self): + """Create test data.""" + from pyresample.geometry import AreaDefinition + area = AreaDefinition('test', 'test', 'test', + {'proj': 'merc'}, 2, 2, + (-2000, -2000, 2000, 2000)) + attrs = {'area': area, + 'start_time': datetime(2018, 1, 1, 18), + 'modifiers': tuple(), + 'name': 'test_vis'} + ds1 = xr.DataArray(da.ones((2, 2), chunks=2, dtype=np.float64), + attrs=attrs, dims=('y', 'x'), + coords={'y': [0, 1], 'x': [0, 1]}) + self.ds1 = ds1 + self.sza = xr.DataArray( + np.rad2deg(np.arccos(da.from_array([[0.0149581333, 0.0146694376], [0.0150812684, 0.0147925727]], + chunks=2))), + attrs={'area': area}, + dims=('y', 'x'), + coords={'y': [0, 1], 'x': [0, 1]}, + ) + + def test_basic_default_not_provided(self): + """Test default limits when SZA isn't provided.""" + from satpy.composites import SunZenithCorrector + comp = SunZenithCorrector(name='sza_test', modifiers=tuple()) + res = comp((self.ds1,), test_attr='test') + np.testing.assert_allclose(res.values, np.array([[22.401667, 22.31777], [22.437503, 22.353533]])) + + def test_basic_lims_not_provided(self): + """Test custom limits when SZA isn't provided.""" + from satpy.composites import SunZenithCorrector + comp = SunZenithCorrector(name='sza_test', modifiers=tuple(), correction_limit=90) + res = comp((self.ds1,), test_attr='test') + np.testing.assert_allclose(res.values, np.array([[66.853262, 68.168939], [66.30742, 67.601493]])) + + def test_basic_default_provided(self): + """Test default limits when SZA is provided.""" + from satpy.composites import SunZenithCorrector + comp = SunZenithCorrector(name='sza_test', modifiers=tuple()) + res = comp((self.ds1, self.sza), test_attr='test') + np.testing.assert_allclose(res.values, np.array([[22.401667, 22.31777], [22.437503, 22.353533]])) + + def test_basic_lims_provided(self): + """Test custom limits when SZA is provided.""" + from satpy.composites import SunZenithCorrector + comp = SunZenithCorrector(name='sza_test', modifiers=tuple(), correction_limit=90) + res = comp((self.ds1, self.sza), test_attr='test') + np.testing.assert_allclose(res.values, np.array([[66.853262, 68.168939], [66.30742, 67.601493]])) + + class TestDayNightCompositor(unittest.TestCase): """Test DayNightCompositor.""" @@ -464,6 +516,7 @@ def suite(): mysuite.addTests(test_viirs.suite()) mysuite.addTest(loader.loadTestsFromTestCase(TestCheckArea)) mysuite.addTest(loader.loadTestsFromTestCase(TestRatioSharpenedCompositors)) + mysuite.addTest(loader.loadTestsFromTestCase(TestSunZenithCorrector)) mysuite.addTest(loader.loadTestsFromTestCase(TestDayNightCompositor)) mysuite.addTest(loader.loadTestsFromTestCase(TestFillingCompositor)) mysuite.addTest(loader.loadTestsFromTestCase(TestSandwichCompositor)) diff --git a/satpy/utils.py b/satpy/utils.py index 8cf063a113..46cbed9aef 100644 --- a/satpy/utils.py +++ b/satpy/utils.py @@ -31,6 +31,7 @@ import os import re +import numpy as np import xarray.ufuncs as xu try: @@ -208,56 +209,82 @@ def proj_units_to_meters(proj_str): def _get_sunz_corr_li_and_shibata(cos_zen): + return 24.35 / (2. * cos_zen + np.sqrt(498.5225 * cos_zen**2 + 1)) - return 24.35 / (2. * cos_zen + - xu.sqrt(498.5225 * cos_zen**2 + 1)) - -def sunzen_corr_cos(data, cos_zen, limit=88.): +def sunzen_corr_cos(data, cos_zen, limit=88., max_sza=95.): """Perform Sun zenith angle correction. The correction is based on the provided cosine of the zenith - angle (*cos_zen*). The correction is limited - to *limit* degrees (default: 88.0 degrees). For larger zenith - angles, the correction is the same as at the *limit*. Both *data* - and *cos_zen* are given as 2-dimensional Numpy arrays or Numpy - MaskedArrays, and they should have equal shapes. + angle (``cos_zen``). The correction is limited + to ``limit`` degrees (default: 88.0 degrees). For larger zenith + angles, the correction is the same as at the ``limit`` if ``max_sza`` + is `None`. The default behavior is to gradually reduce the correction + past ``limit`` degrees up to ``max_sza`` where the correction becomes + 0. Both ``data`` and ``cos_zen`` should be 2D arrays of the same shape. """ # Convert the zenith angle limit to cosine of zenith angle - limit = xu.cos(xu.deg2rad(limit)) + limit_rad = np.deg2rad(limit) + limit_cos = np.cos(limit_rad) + max_sza_rad = np.deg2rad(max_sza) if max_sza is not None else max_sza # Cosine correction corr = 1. / cos_zen - # Use constant value (the limit) for larger zenith - # angles - corr = corr.where(cos_zen > limit).fillna(1 / limit) + if max_sza is not None: + # gradually fall off for larger zenith angle + grad_factor = (np.arccos(cos_zen) - limit_rad) / (max_sza_rad - limit_rad) + # invert the factor so maximum correction is done at `limit` and falls off later + grad_factor = 1. - np.log(grad_factor + 1) / np.log(2) + # make sure we don't make anything negative + grad_factor = grad_factor.clip(0.) + else: + # Use constant value (the limit) for larger zenith angles + grad_factor = 1. + corr = corr.where(cos_zen > limit_cos, grad_factor / limit_cos) + # Force "night" pixels to 0 (where SZA is invalid) + corr = corr.where(cos_zen.notnull(), 0) return data * corr -def atmospheric_path_length_correction(data, cos_zen, limit=88.): +def atmospheric_path_length_correction(data, cos_zen, limit=88., max_sza=95.): """Perform Sun zenith angle correction. This function uses the correction method proposed by Li and Shibata (2006): https://doi.org/10.1175/JAS3682.1 - The correction is limited to *limit* degrees (default: 88.0 degrees). For - larger zenith angles, the correction is the same as at the *limit*. Both - *data* and *cos_zen* are given as 2-dimensional Numpy arrays or Numpy - MaskedArrays, and they should have equal shapes. + The correction is limited to ``limit`` degrees (default: 88.0 degrees). For + larger zenith angles, the correction is the same as at the ``limit`` if + ``max_sza`` is `None`. The default behavior is to gradually reduce the + correction past ``limit`` degrees up to ``max_sza`` where the correction + becomes 0. Both ``data`` and ``cos_zen`` should be 2D arrays of the same + shape. """ - # Convert the zenith angle limit to cosine of zenith angle - limit = xu.cos(xu.radians(limit)) + limit_rad = np.deg2rad(limit) + limit_cos = np.cos(limit_rad) + max_sza_rad = np.deg2rad(max_sza) if max_sza is not None else max_sza # Cosine correction corr = _get_sunz_corr_li_and_shibata(cos_zen) - # Use constant value (the limit) for larger zenith - # angles - corr_lim = _get_sunz_corr_li_and_shibata(limit) - corr = corr.where(cos_zen > limit).fillna(corr_lim) + # Use constant value (the limit) for larger zenith angles + corr_lim = _get_sunz_corr_li_and_shibata(limit_cos) + + if max_sza is not None: + # gradually fall off for larger zenith angle + grad_factor = (np.arccos(cos_zen) - limit_rad) / (max_sza_rad - limit_rad) + # invert the factor so maximum correction is done at `limit` and falls off later + grad_factor = 1. - np.log(grad_factor + 1) / np.log(2) + # make sure we don't make anything negative + grad_factor = grad_factor.clip(0.) + else: + # Use constant value (the limit) for larger zenith angles + grad_factor = 1. + corr = corr.where(cos_zen > limit_cos, grad_factor * corr_lim) + # Force "night" pixels to 0 (where SZA is invalid) + corr = corr.where(cos_zen.notnull(), 0) return data * corr diff --git a/setup.py b/setup.py index 83542db4a5..f49d51b91c 100644 --- a/setup.py +++ b/setup.py @@ -31,7 +31,7 @@ from setuptools import find_packages, setup -requires = ['numpy >=1.12', 'pillow', 'pyresample >=1.10.3', 'trollsift', +requires = ['numpy >=1.13', 'pillow', 'pyresample >=1.10.3', 'trollsift', 'trollimage >=1.5.1', 'pykdtree', 'six', 'pyyaml', 'xarray >=0.10.1', 'dask[array] >=0.17.1']