From 937950dcb375a719d5a542c987bee723176ed9a1 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 24 Mar 2021 14:15:13 -0500 Subject: [PATCH 1/7] Switch VIIRS CREFL ReflectanceCorrector to download DEM file --- satpy/composites/viirs.py | 55 +-- satpy/etc/composites/modis.yaml | 3 +- satpy/etc/composites/viirs.yaml | 6 +- satpy/tests/compositor_tests/test_viirs.py | 396 +++++++++++---------- 4 files changed, 244 insertions(+), 216 deletions(-) diff --git a/satpy/composites/viirs.py b/satpy/composites/viirs.py index 7f6a9c03d1..ae669e7918 100644 --- a/satpy/composites/viirs.py +++ b/satpy/composites/viirs.py @@ -18,22 +18,23 @@ """Composite classes for the VIIRS instrument.""" import logging -import os +import warnings 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.modifiers import ModifierBase from satpy.dataset import combine_metadata from satpy.utils import get_satpos +from satpy.aux_download import DataDownloadMixin, retrieve LOG = logging.getLogger(__name__) -class ReflectanceCorrector(CompositeBase): +class ReflectanceCorrector(ModifierBase, DataDownloadMixin): """Corrected Reflectance (crefl) modifier. Uses a python rewrite of the C CREFL code written for VIIRS and MODIS. @@ -45,21 +46,31 @@ def __init__(self, *args, **kwargs): 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(, ) - :param dem_sds: variable name to load from the ancillary file + 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", - 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) + 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.""" @@ -82,12 +93,13 @@ def __call__(self, datasets, optional_datasets, **info): refl_data = datasets[0] if refl_data.attrs.get("rayleigh_corrected"): return refl_data - if os.path.isfile(self.dem_file): + if self.dem_cache_key is not None: LOG.debug("Loading CREFL averaged elevation information from: %s", - self.dem_file) + 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(self.dem_file, "r") + 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) @@ -98,8 +110,7 @@ def __call__(self, datasets, optional_datasets, **info): from satpy.composites.crefl_utils import run_crefl, get_coefficients - percent = refl_data.attrs["units"] == "%" - + is_percent = refl_data.attrs["units"] == "%" coefficients = get_coefficients(refl_data.attrs["sensor"], refl_data.attrs["wavelength"], refl_data.attrs["resolution"]) @@ -114,11 +125,11 @@ def __call__(self, datasets, optional_datasets, **info): solar_aa, solar_za, avg_elevation=avg_elevation, - percent=percent, + percent=is_percent, use_abi=use_abi) info.update(refl_data.attrs) info["rayleigh_corrected"] = True - factor = 100. if percent else 1. + factor = 100. if is_percent else 1. results = results * factor results.attrs = info self.apply_modifier_info(refl_data, results) diff --git a/satpy/etc/composites/modis.yaml b/satpy/etc/composites/modis.yaml index 3778646b11..5afe057de2 100644 --- a/satpy/etc/composites/modis.yaml +++ b/satpy/etc/composites/modis.yaml @@ -2,7 +2,8 @@ sensor_name: visir/modis modifiers: rayleigh_corrected_crefl: modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector - dem_filename: CMGDEM.hdf + url: "https://www.ssec.wisc.edu/~davidh/polar2grid/viirs_crefl/CMGDEM.hdf" + known_hash: "sha256:f33f1f867d79fff4fafe128f61c154236dd74fcc97bf418ea1437977a38d0604" optional_prerequisites: - name: satellite_azimuth_angle - name: satellite_zenith_angle diff --git a/satpy/etc/composites/viirs.yaml b/satpy/etc/composites/viirs.yaml index 7c14e6c2e1..bcc183705d 100644 --- a/satpy/etc/composites/viirs.yaml +++ b/satpy/etc/composites/viirs.yaml @@ -3,7 +3,8 @@ sensor_name: visir/viirs modifiers: rayleigh_corrected_crefl: modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector - dem_filename: CMGDEM.hdf + url: "https://www.ssec.wisc.edu/~davidh/polar2grid/viirs_crefl/CMGDEM.hdf" + known_hash: "sha256:f33f1f867d79fff4fafe128f61c154236dd74fcc97bf418ea1437977a38d0604" optional_prerequisites: - name: satellite_azimuth_angle resolution: 742 @@ -16,7 +17,8 @@ modifiers: rayleigh_corrected_crefl_iband: modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector - dem_filename: CMGDEM.hdf + url: "https://www.ssec.wisc.edu/~davidh/polar2grid/viirs_crefl/CMGDEM.hdf" + known_hash: "sha256:f33f1f867d79fff4fafe128f61c154236dd74fcc97bf418ea1437977a38d0604" optional_prerequisites: - name: satellite_azimuth_angle resolution: 371 diff --git a/satpy/tests/compositor_tests/test_viirs.py b/satpy/tests/compositor_tests/test_viirs.py index 1cce2356a8..00da0359e9 100644 --- a/satpy/tests/compositor_tests/test_viirs.py +++ b/satpy/tests/compositor_tests/test_viirs.py @@ -20,30 +20,14 @@ import unittest from unittest import mock +import dask.array as da +import numpy as np +from pyresample.geometry import AreaDefinition + class TestVIIRSComposites(unittest.TestCase): """Test VIIRS-specific composites.""" - def data_area_ref_corrector(self): - """Create test area definition and data.""" - import dask.array as da - import numpy as np - from pyresample.geometry import AreaDefinition - rows = 5 - cols = 10 - area = AreaDefinition( - 'some_area_name', 'On-the-fly area', 'geosabii', - {'a': '6378137.0', 'b': '6356752.31414', 'h': '35786023.0', 'lon_0': '-89.5', 'proj': 'geos', 'sweep': 'x', - 'units': 'm'}, - cols, rows, - (-5434894.954752679, -5434894.964451744, 5434894.964451744, 5434894.954752679)) - - dnb = np.zeros((rows, cols)) + 25 - dnb[3, :] += 25 - dnb[4:, :] += 50 - dnb = da.from_array(dnb, chunks=100) - return area, dnb - def test_load_composite_yaml(self): """Test loading the yaml for this sensor.""" from satpy.composites.config_loader import CompositorLoader @@ -246,6 +230,82 @@ def test_hncc_dnb(self): 9.50784532e-03, 1.42617433e-02, 1.50001560e+03, 3.00001560e+03, 4.50001560e+03]) + +class TestViirsReflectanceCorrectorAnglesTest(unittest.TestCase): + """Tests for the VIIRS/MODIS Corrected Reflectance modifier handling angles.""" + + def setUp(self): + """Patch in-class imports.""" + self.astronomy = mock.MagicMock() + self.orbital = mock.MagicMock() + modules = { + 'pyorbital.astronomy': self.astronomy, + 'pyorbital.orbital': self.orbital, + } + self.module_patcher = mock.patch.dict('sys.modules', modules) + self.module_patcher.start() + + def tearDown(self): + """Unpatch in-class imports.""" + self.module_patcher.stop() + + @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() + 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'}) + + # Compute angles + psp = ReflectanceCorrector(name='dummy') + psp.get_angles(vis) + + # Check arguments of get_orbserver_look() call, especially the altitude + # unit conversion from meters to kilometers + 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) + + +class TestReflectanceCorrectorModifier: + """Test the CREFL modifier.""" + + @staticmethod + def data_area_ref_corrector(): + """Create test area definition and data.""" + rows = 5 + cols = 10 + area = AreaDefinition( + 'some_area_name', 'On-the-fly area', 'geosabii', + {'a': '6378137.0', 'b': '6356752.31414', 'h': '35786023.0', 'lon_0': '-89.5', 'proj': 'geos', 'sweep': 'x', + 'units': 'm'}, + cols, rows, + (-5434894.954752679, -5434894.964451744, 5434894.964451744, 5434894.954752679)) + + dnb = np.zeros((rows, cols)) + 25 + dnb[3, :] += 25 + dnb[4:, :] += 50 + dnb = da.from_array(dnb, chunks=100) + return area, dnb + def test_reflectance_corrector_abi(self): """Test ReflectanceCorrector modifier with ABI data.""" import xarray as xr @@ -261,50 +321,52 @@ def test_reflectance_corrector_abi(self): wavelength=(0.45, 0.47, 0.49), resolution=1000, calibration='reflectance', modifiers=('sunz_corrected', 'rayleigh_corrected_crefl',), sensor='abi') - self.assertEqual(ref_cor.attrs['modifiers'], ('sunz_corrected', 'rayleigh_corrected_crefl',)) - self.assertEqual(ref_cor.attrs['calibration'], 'reflectance') - self.assertEqual(ref_cor.attrs['wavelength'], (0.45, 0.47, 0.49)) - self.assertEqual(ref_cor.attrs['name'], 'C01') - self.assertEqual(ref_cor.attrs['resolution'], 1000) - self.assertEqual(ref_cor.attrs['sensor'], 'abi') - self.assertEqual(ref_cor.attrs['prerequisites'], []) - self.assertEqual(ref_cor.attrs['optional_prerequisites'], [ + assert ref_cor.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl') + assert ref_cor.attrs['calibration'] == 'reflectance' + assert ref_cor.attrs['wavelength'] == (0.45, 0.47, 0.49) + assert ref_cor.attrs['name'] == 'C01' + assert ref_cor.attrs['resolution'] == 1000 + assert ref_cor.attrs['sensor'] == 'abi' + assert ref_cor.attrs['prerequisites'] == [] + assert ref_cor.attrs['optional_prerequisites'] == [ make_dsq(name='satellite_azimuth_angle'), make_dsq(name='satellite_zenith_angle'), make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle')]) + make_dsq(name='solar_zenith_angle')] area, dnb = self.data_area_ref_corrector() c01 = xr.DataArray(dnb, dims=('y', 'x'), - attrs={'satellite_longitude': -89.5, 'satellite_latitude': 0.0, - 'satellite_altitude': 35786023.4375, 'platform_name': 'GOES-16', - 'calibration': 'reflectance', 'units': '%', 'wavelength': (0.45, 0.47, 0.49), - 'name': 'C01', 'resolution': 1000, 'sensor': 'abi', - 'start_time': '2017-09-20 17:30:40.800000', 'end_time': '2017-09-20 17:41:17.500000', - 'area': area, 'ancillary_variables': []}) + attrs={ + 'satellite_longitude': -89.5, 'satellite_latitude': 0.0, + 'satellite_altitude': 35786023.4375, 'platform_name': 'GOES-16', + 'calibration': 'reflectance', 'units': '%', 'wavelength': (0.45, 0.47, 0.49), + 'name': 'C01', 'resolution': 1000, 'sensor': 'abi', + 'start_time': '2017-09-20 17:30:40.800000', 'end_time': '2017-09-20 17:41:17.500000', + 'area': area, 'ancillary_variables': [] + }) res = ref_cor([c01], []) - self.assertIsInstance(res, xr.DataArray) - self.assertIsInstance(res.data, da.Array) - self.assertEqual(res.attrs['satellite_longitude'], -89.5) - self.assertEqual(res.attrs['satellite_latitude'], 0.0) - self.assertEqual(res.attrs['satellite_altitude'], 35786023.4375) - self.assertEqual(res.attrs['modifiers'], ('sunz_corrected', 'rayleigh_corrected_crefl',)) - self.assertEqual(res.attrs['platform_name'], 'GOES-16') - self.assertEqual(res.attrs['calibration'], 'reflectance') - self.assertEqual(res.attrs['units'], '%') - self.assertEqual(res.attrs['wavelength'], (0.45, 0.47, 0.49)) - self.assertEqual(res.attrs['name'], 'C01') - self.assertEqual(res.attrs['resolution'], 1000) - self.assertEqual(res.attrs['sensor'], 'abi') - self.assertEqual(res.attrs['start_time'], '2017-09-20 17:30:40.800000') - self.assertEqual(res.attrs['end_time'], '2017-09-20 17:41:17.500000') - self.assertEqual(res.attrs['area'], area) - self.assertEqual(res.attrs['ancillary_variables'], []) + assert isinstance(res, xr.DataArray) + assert isinstance(res.data, da.Array) + assert res.attrs['satellite_longitude'] == -89.5 + assert res.attrs['satellite_latitude'] == 0.0 + assert res.attrs['satellite_altitude'] == 35786023.4375 + assert res.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl') + assert res.attrs['platform_name'] == 'GOES-16' + assert res.attrs['calibration'] == 'reflectance' + assert res.attrs['units'] == '%' + assert res.attrs['wavelength'] == (0.45, 0.47, 0.49) + assert res.attrs['name'] == 'C01' + assert res.attrs['resolution'] == 1000 + assert res.attrs['sensor'] == 'abi' + assert res.attrs['start_time'] == '2017-09-20 17:30:40.800000' + assert res.attrs['end_time'] == '2017-09-20 17:41:17.500000' + assert res.attrs['area'] == area + assert res.attrs['ancillary_variables'] == [] data = res.values - self.assertLess(abs(np.nanmean(data) - 26.00760944144745), 1e-10) - self.assertEqual(data.shape, (5, 10)) + assert abs(np.nanmean(data) - 26.00760944144745) < 1e-10 + assert data.shape == (5, 10) 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, @@ -328,72 +390,75 @@ def test_reflectance_corrector_viirs(self): from satpy.composites.viirs import ReflectanceCorrector from satpy.tests.utils import make_dsq ref_cor = ReflectanceCorrector(dem_filename='_fake.hdf', optional_prerequisites=[ - make_dsq(name='satellite_azimuth_angle'), - make_dsq(name='satellite_zenith_angle'), - make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle')], + make_dsq(name='satellite_azimuth_angle'), + make_dsq(name='satellite_zenith_angle'), + make_dsq(name='solar_azimuth_angle'), + make_dsq(name='solar_zenith_angle')], name='I01', prerequisites=[], wavelength=(0.6, 0.64, 0.68), resolution=371, calibration='reflectance', modifiers=('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband'), sensor='viirs') - self.assertEqual(ref_cor.attrs['modifiers'], ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband')) - self.assertEqual(ref_cor.attrs['calibration'], 'reflectance') - self.assertEqual(ref_cor.attrs['wavelength'], (0.6, 0.64, 0.68)) - self.assertEqual(ref_cor.attrs['name'], 'I01') - self.assertEqual(ref_cor.attrs['resolution'], 371) - self.assertEqual(ref_cor.attrs['sensor'], 'viirs') - self.assertEqual(ref_cor.attrs['prerequisites'], []) - self.assertEqual(ref_cor.attrs['optional_prerequisites'], [ + assert ref_cor.attrs['modifiers'] == ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband') + assert ref_cor.attrs['calibration'] == 'reflectance' + assert ref_cor.attrs['wavelength'] == (0.6, 0.64, 0.68) + assert ref_cor.attrs['name'] == 'I01' + assert ref_cor.attrs['resolution'] == 371 + assert ref_cor.attrs['sensor'] == 'viirs' + assert ref_cor.attrs['prerequisites'] == [] + assert ref_cor.attrs['optional_prerequisites'] == [ make_dsq(name='satellite_azimuth_angle'), make_dsq(name='satellite_zenith_angle'), make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle')]) + make_dsq(name='solar_zenith_angle')] area, dnb = self.data_area_ref_corrector() - def make_xarray(self, file_key, name, standard_name, wavelength=None, units='degrees', calibration=None, + def make_xarray(file_key, name, standard_name, wavelength=None, units='degrees', calibration=None, file_type=('gitco', 'gimgo')): return xr.DataArray(dnb, dims=('y', 'x'), - attrs={'start_orbit': 1708, 'end_orbit': 1708, 'wavelength': wavelength, 'level': None, - 'modifiers': None, 'calibration': calibration, 'file_key': file_key, - 'resolution': 371, 'file_type': file_type, 'name': name, - 'standard_name': standard_name, 'platform_name': 'Suomi-NPP', - 'polarization': None, 'sensor': 'viirs', 'units': units, - 'start_time': datetime.datetime(2012, 2, 25, 18, 1, 24, 570942), - 'end_time': datetime.datetime(2012, 2, 25, 18, 11, 21, 175760), 'area': area, - 'ancillary_variables': []}) - c01 = make_xarray(self, None, 'I01', 'toa_bidirectional_reflectance', wavelength=(0.6, 0.64, 0.68), units='%', + attrs={ + 'start_orbit': 1708, 'end_orbit': 1708, 'wavelength': wavelength, 'level': None, + 'modifiers': None, 'calibration': calibration, 'file_key': file_key, + 'resolution': 371, 'file_type': file_type, 'name': name, + 'standard_name': standard_name, 'platform_name': 'Suomi-NPP', + 'polarization': None, 'sensor': 'viirs', 'units': units, + 'start_time': datetime.datetime(2012, 2, 25, 18, 1, 24, 570942), + 'end_time': datetime.datetime(2012, 2, 25, 18, 11, 21, 175760), 'area': area, + 'ancillary_variables': [] + }) + + c01 = make_xarray(None, 'I01', 'toa_bidirectional_reflectance', wavelength=(0.6, 0.64, 0.68), units='%', calibration='reflectance', file_type='svi01') - c02 = make_xarray(self, 'All_Data/{file_group}_All/SatelliteAzimuthAngle', 'satellite_azimuth_angle', + c02 = make_xarray('All_Data/{file_group}_All/SatelliteAzimuthAngle', 'satellite_azimuth_angle', 'sensor_azimuth_angle') - c03 = make_xarray(self, 'All_Data/{file_group}_All/SatelliteZenithAngle', 'satellite_zenith_angle', + c03 = make_xarray('All_Data/{file_group}_All/SatelliteZenithAngle', 'satellite_zenith_angle', 'sensor_zenith_angle') - c04 = make_xarray(self, 'All_Data/{file_group}_All/SolarAzimuthAngle', 'solar_azimuth_angle', + c04 = make_xarray('All_Data/{file_group}_All/SolarAzimuthAngle', 'solar_azimuth_angle', 'solar_azimuth_angle') - c05 = make_xarray(self, 'All_Data/{file_group}_All/SolarZenithAngle', 'solar_zenith_angle', + c05 = make_xarray('All_Data/{file_group}_All/SolarZenithAngle', 'solar_zenith_angle', 'solar_zenith_angle') res = ref_cor([c01], [c02, c03, c04, c05]) - self.assertIsInstance(res, xr.DataArray) - self.assertIsInstance(res.data, da.Array) - self.assertEqual(res.attrs['wavelength'], (0.6, 0.64, 0.68)) - self.assertEqual(res.attrs['modifiers'], ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband')) - self.assertEqual(res.attrs['calibration'], 'reflectance') - self.assertEqual(res.attrs['resolution'], 371) - self.assertEqual(res.attrs['file_type'], 'svi01') - self.assertEqual(res.attrs['name'], 'I01') - self.assertEqual(res.attrs['standard_name'], 'toa_bidirectional_reflectance') - self.assertEqual(res.attrs['platform_name'], 'Suomi-NPP') - self.assertEqual(res.attrs['sensor'], 'viirs') - self.assertEqual(res.attrs['units'], '%') - self.assertEqual(res.attrs['start_time'], datetime.datetime(2012, 2, 25, 18, 1, 24, 570942)) - self.assertEqual(res.attrs['end_time'], datetime.datetime(2012, 2, 25, 18, 11, 21, 175760)) - self.assertEqual(res.attrs['area'], area) - self.assertEqual(res.attrs['ancillary_variables'], []) + assert isinstance(res, xr.DataArray) + assert isinstance(res.data, da.Array) + assert res.attrs['wavelength'] == (0.6, 0.64, 0.68) + assert res.attrs['modifiers'] == ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband') + assert res.attrs['calibration'] == 'reflectance' + assert res.attrs['resolution'] == 371 + assert res.attrs['file_type'] == 'svi01' + assert res.attrs['name'] == 'I01' + assert res.attrs['standard_name'] == 'toa_bidirectional_reflectance' + assert res.attrs['platform_name'] == 'Suomi-NPP' + assert res.attrs['sensor'] == 'viirs' + assert res.attrs['units'] == '%' + assert res.attrs['start_time'] == datetime.datetime(2012, 2, 25, 18, 1, 24, 570942) + assert res.attrs['end_time'] == datetime.datetime(2012, 2, 25, 18, 11, 21, 175760) + assert res.attrs['area'] == area + assert res.attrs['ancillary_variables'] == [] data = res.values - self.assertLess(abs(np.mean(data) - 40.7578684169142), 1e-10) - self.assertEqual(data.shape, (5, 10)) + assert abs(np.mean(data) - 40.7578684169142) < 1e-10 + assert data.shape == (5, 10) unique = np.unique(data) np.testing.assert_allclose(unique, [25.20341702519979, 52.38819447051263, 75.79089653845898]) @@ -413,112 +478,61 @@ def test_reflectance_corrector_modis(self): dem_filename='_fake.hdf', optional_prerequisites=[sataa_did, satza_did, solaa_did, solza_did], name='1', prerequisites=[], wavelength=(0.62, 0.645, 0.67), resolution=250, calibration='reflectance', modifiers=('sunz_corrected', 'rayleigh_corrected_crefl'), sensor='modis') - self.assertEqual(ref_cor.attrs['modifiers'], ('sunz_corrected', 'rayleigh_corrected_crefl')) - self.assertEqual(ref_cor.attrs['calibration'], 'reflectance') - self.assertEqual(ref_cor.attrs['wavelength'], (0.62, 0.645, 0.67)) - self.assertEqual(ref_cor.attrs['name'], '1') - self.assertEqual(ref_cor.attrs['resolution'], 250) - self.assertEqual(ref_cor.attrs['sensor'], 'modis') - self.assertEqual(ref_cor.attrs['prerequisites'], []) - self.assertEqual(ref_cor.attrs['optional_prerequisites'], [ + assert ref_cor.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl') + assert ref_cor.attrs['calibration'] == 'reflectance' + assert ref_cor.attrs['wavelength'] == (0.62, 0.645, 0.67) + assert ref_cor.attrs['name'] == '1' + assert ref_cor.attrs['resolution'] == 250 + assert ref_cor.attrs['sensor'] == 'modis' + assert ref_cor.attrs['prerequisites'] == [] + assert ref_cor.attrs['optional_prerequisites'] == [ make_dsq(name='satellite_azimuth_angle'), make_dsq(name='satellite_zenith_angle'), make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle')]) + make_dsq(name='solar_zenith_angle')] area, dnb = self.data_area_ref_corrector() - def make_xarray(self, name, calibration, wavelength=None, modifiers=None, resolution=1000, + def make_xarray(name, calibration, wavelength=None, modifiers=None, resolution=1000, file_type='hdf_eos_geo'): return xr.DataArray(dnb, dims=('y', 'x'), - attrs={'wavelength': wavelength, 'level': None, 'modifiers': modifiers, - 'calibration': calibration, 'resolution': resolution, 'file_type': file_type, - 'name': name, 'coordinates': ['longitude', 'latitude'], - 'platform_name': 'EOS-Aqua', 'polarization': None, 'sensor': 'modis', - 'units': '%', 'start_time': datetime.datetime(2012, 8, 13, 18, 46, 1, 439838), - 'end_time': datetime.datetime(2012, 8, 13, 18, 57, 47, 746296), 'area': area, - 'ancillary_variables': []}) - c01 = make_xarray(self, '1', 'reflectance', wavelength=(0.62, 0.645, 0.67), modifiers='sunz_corrected', + attrs={ + 'wavelength': wavelength, 'level': None, 'modifiers': modifiers, + 'calibration': calibration, 'resolution': resolution, 'file_type': file_type, + 'name': name, 'coordinates': ['longitude', 'latitude'], + 'platform_name': 'EOS-Aqua', 'polarization': None, 'sensor': 'modis', + 'units': '%', 'start_time': datetime.datetime(2012, 8, 13, 18, 46, 1, 439838), + 'end_time': datetime.datetime(2012, 8, 13, 18, 57, 47, 746296), 'area': area, + 'ancillary_variables': [] + }) + + c01 = make_xarray('1', 'reflectance', wavelength=(0.62, 0.645, 0.67), modifiers='sunz_corrected', resolution=500, file_type='hdf_eos_data_500m') - c02 = make_xarray(self, 'satellite_azimuth_angle', None) - c03 = make_xarray(self, 'satellite_zenith_angle', None) - c04 = make_xarray(self, 'solar_azimuth_angle', None) - c05 = make_xarray(self, 'solar_zenith_angle', None) + c02 = make_xarray('satellite_azimuth_angle', None) + c03 = make_xarray('satellite_zenith_angle', None) + c04 = make_xarray('solar_azimuth_angle', None) + c05 = make_xarray('solar_zenith_angle', None) res = ref_cor([c01], [c02, c03, c04, c05]) - self.assertIsInstance(res, xr.DataArray) - self.assertIsInstance(res.data, da.Array) - self.assertEqual(res.attrs['wavelength'], (0.62, 0.645, 0.67)) - self.assertEqual(res.attrs['modifiers'], ('sunz_corrected', 'rayleigh_corrected_crefl',)) - self.assertEqual(res.attrs['calibration'], 'reflectance') - self.assertEqual(res.attrs['resolution'], 500) - self.assertEqual(res.attrs['file_type'], 'hdf_eos_data_500m') - self.assertEqual(res.attrs['name'], '1') - self.assertEqual(res.attrs['platform_name'], 'EOS-Aqua') - self.assertEqual(res.attrs['sensor'], 'modis') - self.assertEqual(res.attrs['units'], '%') - self.assertEqual(res.attrs['start_time'], datetime.datetime(2012, 8, 13, 18, 46, 1, 439838)) - self.assertEqual(res.attrs['end_time'], datetime.datetime(2012, 8, 13, 18, 57, 47, 746296)) - self.assertEqual(res.attrs['area'], area) - self.assertEqual(res.attrs['ancillary_variables'], []) + assert isinstance(res, xr.DataArray) + assert isinstance(res.data, da.Array) + assert res.attrs['wavelength'] == (0.62, 0.645, 0.67) + assert res.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl',) + assert res.attrs['calibration'] == 'reflectance' + assert res.attrs['resolution'] == 500 + assert res.attrs['file_type'] == 'hdf_eos_data_500m' + assert res.attrs['name'] == '1' + assert res.attrs['platform_name'] == 'EOS-Aqua' + assert res.attrs['sensor'] == 'modis' + assert res.attrs['units'] == '%' + assert res.attrs['start_time'] == datetime.datetime(2012, 8, 13, 18, 46, 1, 439838) + assert res.attrs['end_time'] == datetime.datetime(2012, 8, 13, 18, 57, 47, 746296) + assert res.attrs['area'] == area + assert res.attrs['ancillary_variables'] == [] data = res.values if abs(np.mean(data) - 38.734365117099145) >= 1e-10: raise AssertionError('{} is not within {} of {}'.format(np.mean(data), 1e-10, 38.734365117099145)) - self.assertEqual(data.shape, (5, 10)) + assert data.shape == (5, 10) unique = np.unique(data) np.testing.assert_allclose(unique, [24.641586, 50.431692, 69.315375]) - - -class ViirsReflectanceCorrectorTest(unittest.TestCase): - """Tests for the VIIRS/MODIS Corrected Reflectance modifier.""" - - def setUp(self): - """Patch in-class imports.""" - self.astronomy = mock.MagicMock() - self.orbital = mock.MagicMock() - modules = { - 'pyorbital.astronomy': self.astronomy, - 'pyorbital.orbital': self.orbital, - } - self.module_patcher = mock.patch.dict('sys.modules', modules) - self.module_patcher.start() - - def tearDown(self): - """Unpatch in-class imports.""" - self.module_patcher.stop() - - @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() - 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'}) - - # Compute angles - psp = ReflectanceCorrector(name='dummy') - psp.get_angles(vis) - - # Check arguments of get_orbserver_look() call, especially the altitude - # unit conversion from meters to kilometers - 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) From 779d7a336b7d4b4eeb839243e846622059b59a75 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 24 Mar 2021 19:45:41 -0500 Subject: [PATCH 2/7] Refactor crefl reflectance corrector modifier --- satpy/composites/viirs.py | 84 ++++++++++++++++++++++----------------- 1 file changed, 48 insertions(+), 36 deletions(-) diff --git a/satpy/composites/viirs.py b/satpy/composites/viirs.py index ae669e7918..c467e1e454 100644 --- a/satpy/composites/viirs.py +++ b/satpy/composites/viirs.py @@ -74,48 +74,17 @@ def __init__(self, *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 self.dem_cache_key is not None: - 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) - else: - avg_elevation = None - 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 = vis.attrs['sensor'] == 'abi' - lons, lats = vis.attrs['area'].get_lonlats(chunks=vis.chunks) + 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, @@ -135,6 +104,49 @@ def __call__(self, datasets, optional_datasets, **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 aa0e556889cc18d8e9e59514b48a07dc6ad10eb5 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 24 Mar 2021 20:06:04 -0500 Subject: [PATCH 3/7] Change viirs and modis rayleigh correction to use required prereqs --- satpy/etc/composites/modis.yaml | 6 +++--- satpy/etc/composites/viirs.yaml | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/satpy/etc/composites/modis.yaml b/satpy/etc/composites/modis.yaml index 5afe057de2..b4f775aa62 100644 --- a/satpy/etc/composites/modis.yaml +++ b/satpy/etc/composites/modis.yaml @@ -2,9 +2,9 @@ sensor_name: visir/modis modifiers: rayleigh_corrected_crefl: modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector - url: "https://www.ssec.wisc.edu/~davidh/polar2grid/viirs_crefl/CMGDEM.hdf" - known_hash: "sha256:f33f1f867d79fff4fafe128f61c154236dd74fcc97bf418ea1437977a38d0604" - optional_prerequisites: + 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 diff --git a/satpy/etc/composites/viirs.yaml b/satpy/etc/composites/viirs.yaml index bcc183705d..44fde79594 100644 --- a/satpy/etc/composites/viirs.yaml +++ b/satpy/etc/composites/viirs.yaml @@ -5,7 +5,7 @@ modifiers: modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector url: "https://www.ssec.wisc.edu/~davidh/polar2grid/viirs_crefl/CMGDEM.hdf" known_hash: "sha256:f33f1f867d79fff4fafe128f61c154236dd74fcc97bf418ea1437977a38d0604" - optional_prerequisites: + prerequisites: - name: satellite_azimuth_angle resolution: 742 - name: satellite_zenith_angle @@ -19,7 +19,7 @@ modifiers: modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector url: "https://www.ssec.wisc.edu/~davidh/polar2grid/viirs_crefl/CMGDEM.hdf" known_hash: "sha256:f33f1f867d79fff4fafe128f61c154236dd74fcc97bf418ea1437977a38d0604" - optional_prerequisites: + prerequisites: - name: satellite_azimuth_angle resolution: 371 - name: satellite_zenith_angle From 1a67341a9bc1074cc6149e49f62fba78ea3c2e23 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 25 Mar 2021 06:59:00 -0500 Subject: [PATCH 4/7] Add CREFL test that uses a DEM file --- satpy/tests/compositor_tests/test_viirs.py | 95 +++++++++++++++++++++- 1 file changed, 92 insertions(+), 3 deletions(-) diff --git a/satpy/tests/compositor_tests/test_viirs.py b/satpy/tests/compositor_tests/test_viirs.py index 00da0359e9..167066a163 100644 --- a/satpy/tests/compositor_tests/test_viirs.py +++ b/satpy/tests/compositor_tests/test_viirs.py @@ -20,6 +20,7 @@ import unittest from unittest import mock +import pytest import dask.array as da import numpy as np from pyresample.geometry import AreaDefinition @@ -313,7 +314,7 @@ def test_reflectance_corrector_abi(self): import numpy as np from satpy.composites.viirs import ReflectanceCorrector from satpy.tests.utils import make_dsq - ref_cor = ReflectanceCorrector(dem_filename='_fake.hdf', optional_prerequisites=[ + ref_cor = ReflectanceCorrector(optional_prerequisites=[ make_dsq(name='satellite_azimuth_angle'), make_dsq(name='satellite_zenith_angle'), make_dsq(name='solar_azimuth_angle'), @@ -389,7 +390,7 @@ def test_reflectance_corrector_viirs(self): import datetime from satpy.composites.viirs import ReflectanceCorrector from satpy.tests.utils import make_dsq - ref_cor = ReflectanceCorrector(dem_filename='_fake.hdf', optional_prerequisites=[ + ref_cor = ReflectanceCorrector(optional_prerequisites=[ make_dsq(name='satellite_azimuth_angle'), make_dsq(name='satellite_zenith_angle'), make_dsq(name='solar_azimuth_angle'), @@ -475,7 +476,7 @@ def test_reflectance_corrector_modis(self): solaa_did = make_dsq(name='solar_azimuth_angle') solza_did = make_dsq(name='solar_zenith_angle') ref_cor = ReflectanceCorrector( - dem_filename='_fake.hdf', optional_prerequisites=[sataa_did, satza_did, solaa_did, solza_did], name='1', + optional_prerequisites=[sataa_did, satza_did, solaa_did, solza_did], name='1', prerequisites=[], wavelength=(0.62, 0.645, 0.67), resolution=250, calibration='reflectance', modifiers=('sunz_corrected', 'rayleigh_corrected_crefl'), sensor='modis') assert ref_cor.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl') @@ -536,3 +537,91 @@ def make_xarray(name, calibration, wavelength=None, modifiers=None, resolution=1 assert data.shape == (5, 10) unique = np.unique(data) np.testing.assert_allclose(unique, [24.641586, 50.431692, 69.315375]) + + def test_reflectance_corrector_bad_prereqs(self): + """Test ReflectanceCorrector modifier with wrong number of inputs.""" + from satpy.composites.viirs import ReflectanceCorrector + ref_cor = ReflectanceCorrector("test") + pytest.raises(ValueError, ref_cor, [1], [2, 3, 4]) + pytest.raises(ValueError, ref_cor, [1, 2, 3, 4], []) + pytest.raises(ValueError, ref_cor, [], [1, 2, 3, 4]) + + def test_crefl_viirs_dem(self, tmpdir): + """Test using the modifier with a DEM file.""" + import xarray as xr + import dask.array as da + import numpy as np + import datetime + from satpy.composites.viirs import ReflectanceCorrector + from satpy.tests.utils import make_dsq + ref_cor = ReflectanceCorrector( + url='CMGDEM.hdf', + optional_prerequisites=[ + make_dsq(name='satellite_azimuth_angle'), + make_dsq(name='satellite_zenith_angle'), + make_dsq(name='solar_azimuth_angle'), + make_dsq(name='solar_zenith_angle') + ], + name='I01', prerequisites=[], wavelength=(0.6, 0.64, 0.68), + resolution=371, calibration='reflectance', + modifiers=('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband'), + sensor='viirs') + area, refl = self.data_area_ref_corrector() + + def make_xarray(file_key, name, standard_name, wavelength=None, units='degrees', calibration=None, + file_type=('gitco', 'gimgo')): + return xr.DataArray(refl, dims=('y', 'x'), + attrs={ + 'start_orbit': 1708, 'end_orbit': 1708, 'wavelength': wavelength, 'level': None, + 'modifiers': None, 'calibration': calibration, 'file_key': file_key, + 'resolution': 371, 'file_type': file_type, 'name': name, + 'standard_name': standard_name, 'platform_name': 'Suomi-NPP', + 'polarization': None, 'sensor': 'viirs', 'units': units, + 'start_time': datetime.datetime(2012, 2, 25, 18, 1, 24, 570942), + 'end_time': datetime.datetime(2012, 2, 25, 18, 11, 21, 175760), 'area': area, + 'ancillary_variables': [] + }) + + c01 = make_xarray(None, 'I01', 'toa_bidirectional_reflectance', wavelength=(0.6, 0.64, 0.68), units='%', + calibration='reflectance', file_type='svi01') + c02 = make_xarray('All_Data/{file_group}_All/SatelliteAzimuthAngle', 'satellite_azimuth_angle', + 'sensor_azimuth_angle') + c03 = make_xarray('All_Data/{file_group}_All/SatelliteZenithAngle', 'satellite_zenith_angle', + 'sensor_zenith_angle') + c04 = make_xarray('All_Data/{file_group}_All/SolarAzimuthAngle', 'solar_azimuth_angle', + 'solar_azimuth_angle') + c05 = make_xarray('All_Data/{file_group}_All/SolarZenithAngle', 'solar_zenith_angle', + 'solar_zenith_angle') + with mock.patch('satpy.composites.viirs.retrieve') as rmock: + dem_fn = str(tmpdir.join('CMGDEM.hdf')) + rmock.return_value = dem_fn + from netCDF4 import Dataset + nc = Dataset(dem_fn, 'w') + nc.createDimension('y', 10) + nc.createDimension('x', 10) + dem_var = nc.createVariable('averaged elevation', np.float32, + ('y', 'x')) + dem_var[:] = 0 + res = ref_cor([c01], [c02, c03, c04, c05]) + + assert isinstance(res, xr.DataArray) + assert isinstance(res.data, da.Array) + assert res.attrs['wavelength'] == (0.6, 0.64, 0.68) + assert res.attrs['modifiers'] == ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband') + assert res.attrs['calibration'] == 'reflectance' + assert res.attrs['resolution'] == 371 + assert res.attrs['file_type'] == 'svi01' + assert res.attrs['name'] == 'I01' + assert res.attrs['standard_name'] == 'toa_bidirectional_reflectance' + assert res.attrs['platform_name'] == 'Suomi-NPP' + assert res.attrs['sensor'] == 'viirs' + assert res.attrs['units'] == '%' + assert res.attrs['start_time'] == datetime.datetime(2012, 2, 25, 18, 1, 24, 570942) + assert res.attrs['end_time'] == datetime.datetime(2012, 2, 25, 18, 11, 21, 175760) + assert res.attrs['area'] == area + assert res.attrs['ancillary_variables'] == [] + data = res.values + assert abs(np.mean(data) - 40.7578684169142) < 1e-10 + assert data.shape == (5, 10) + unique = np.unique(data) + np.testing.assert_allclose(unique, [25.20341702519979, 52.38819447051263, 75.79089653845898]) From 086127e817ae250d7c61f37ad1a9e0c13ad59c6a Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 25 Mar 2021 14:01:04 -0500 Subject: [PATCH 5/7] Refactor crefl tests to be a little cleaner --- satpy/tests/compositor_tests/test_viirs.py | 155 +++++++-------------- 1 file changed, 51 insertions(+), 104 deletions(-) diff --git a/satpy/tests/compositor_tests/test_viirs.py b/satpy/tests/compositor_tests/test_viirs.py index 167066a163..f2bb724aaf 100644 --- a/satpy/tests/compositor_tests/test_viirs.py +++ b/satpy/tests/compositor_tests/test_viirs.py @@ -382,7 +382,8 @@ def test_reflectance_corrector_abi(self): 51.909142813383916, 58.8234273736508, 68.84706145641482, 69.91085190887961, 71.10179768327806, 71.33161009169649]) - def test_reflectance_corrector_viirs(self): + @pytest.mark.parametrize('url', [None, 'CMGDEM.hdf']) + def test_reflectance_corrector_viirs(self, tmpdir, url): """Test ReflectanceCorrector modifier with VIIRS data.""" import xarray as xr import dask.array as da @@ -390,15 +391,22 @@ def test_reflectance_corrector_viirs(self): import datetime from satpy.composites.viirs import ReflectanceCorrector from satpy.tests.utils import make_dsq - ref_cor = ReflectanceCorrector(optional_prerequisites=[ - make_dsq(name='satellite_azimuth_angle'), - make_dsq(name='satellite_zenith_angle'), - make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle')], - name='I01', prerequisites=[], wavelength=(0.6, 0.64, 0.68), resolution=371, - calibration='reflectance', modifiers=('sunz_corrected_iband', - 'rayleigh_corrected_crefl_iband'), - sensor='viirs') + + ref_cor = ReflectanceCorrector( + optional_prerequisites=[ + make_dsq(name='satellite_azimuth_angle'), + make_dsq(name='satellite_zenith_angle'), + make_dsq(name='solar_azimuth_angle'), + make_dsq(name='solar_zenith_angle') + ], + name='I01', + prerequisites=[], + wavelength=(0.6, 0.64, 0.68), + resolution=371, + calibration='reflectance', + modifiers=('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband'), + sensor='viirs', + url=url) assert ref_cor.attrs['modifiers'] == ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband') assert ref_cor.attrs['calibration'] == 'reflectance' @@ -415,13 +423,12 @@ def test_reflectance_corrector_viirs(self): area, dnb = self.data_area_ref_corrector() - def make_xarray(file_key, name, standard_name, wavelength=None, units='degrees', calibration=None, - file_type=('gitco', 'gimgo')): + def make_xarray(name, standard_name, wavelength=None, units='degrees', calibration=None): return xr.DataArray(dnb, dims=('y', 'x'), attrs={ 'start_orbit': 1708, 'end_orbit': 1708, 'wavelength': wavelength, 'level': None, - 'modifiers': None, 'calibration': calibration, 'file_key': file_key, - 'resolution': 371, 'file_type': file_type, 'name': name, + 'modifiers': None, 'calibration': calibration, + 'resolution': 371, 'name': name, 'standard_name': standard_name, 'platform_name': 'Suomi-NPP', 'polarization': None, 'sensor': 'viirs', 'units': units, 'start_time': datetime.datetime(2012, 2, 25, 18, 1, 24, 570942), @@ -429,17 +436,17 @@ def make_xarray(file_key, name, standard_name, wavelength=None, units='degrees', 'ancillary_variables': [] }) - c01 = make_xarray(None, 'I01', 'toa_bidirectional_reflectance', wavelength=(0.6, 0.64, 0.68), units='%', - calibration='reflectance', file_type='svi01') - c02 = make_xarray('All_Data/{file_group}_All/SatelliteAzimuthAngle', 'satellite_azimuth_angle', - 'sensor_azimuth_angle') - c03 = make_xarray('All_Data/{file_group}_All/SatelliteZenithAngle', 'satellite_zenith_angle', - 'sensor_zenith_angle') - c04 = make_xarray('All_Data/{file_group}_All/SolarAzimuthAngle', 'solar_azimuth_angle', - 'solar_azimuth_angle') - c05 = make_xarray('All_Data/{file_group}_All/SolarZenithAngle', 'solar_zenith_angle', - 'solar_zenith_angle') + c01 = make_xarray('I01', 'toa_bidirectional_reflectance', + wavelength=(0.6, 0.64, 0.68), units='%', + calibration='reflectance') + c02 = make_xarray('satellite_azimuth_angle', 'sensor_azimuth_angle') + c03 = make_xarray('satellite_zenith_angle', 'sensor_zenith_angle') + c04 = make_xarray('solar_azimuth_angle', 'solar_azimuth_angle') + c05 = make_xarray('solar_zenith_angle', 'solar_zenith_angle') + + rmock_obj = self._start_dem_mock(tmpdir, url) res = ref_cor([c01], [c02, c03, c04, c05]) + self._stop_dem_mock(rmock_obj) assert isinstance(res, xr.DataArray) assert isinstance(res.data, da.Array) @@ -447,7 +454,6 @@ def make_xarray(file_key, name, standard_name, wavelength=None, units='degrees', assert res.attrs['modifiers'] == ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband') assert res.attrs['calibration'] == 'reflectance' assert res.attrs['resolution'] == 371 - assert res.attrs['file_type'] == 'svi01' assert res.attrs['name'] == 'I01' assert res.attrs['standard_name'] == 'toa_bidirectional_reflectance' assert res.attrs['platform_name'] == 'Suomi-NPP' @@ -546,82 +552,23 @@ def test_reflectance_corrector_bad_prereqs(self): pytest.raises(ValueError, ref_cor, [1, 2, 3, 4], []) pytest.raises(ValueError, ref_cor, [], [1, 2, 3, 4]) - def test_crefl_viirs_dem(self, tmpdir): - """Test using the modifier with a DEM file.""" - import xarray as xr - import dask.array as da - import numpy as np - import datetime - from satpy.composites.viirs import ReflectanceCorrector - from satpy.tests.utils import make_dsq - ref_cor = ReflectanceCorrector( - url='CMGDEM.hdf', - optional_prerequisites=[ - make_dsq(name='satellite_azimuth_angle'), - make_dsq(name='satellite_zenith_angle'), - make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle') - ], - name='I01', prerequisites=[], wavelength=(0.6, 0.64, 0.68), - resolution=371, calibration='reflectance', - modifiers=('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband'), - sensor='viirs') - area, refl = self.data_area_ref_corrector() - - def make_xarray(file_key, name, standard_name, wavelength=None, units='degrees', calibration=None, - file_type=('gitco', 'gimgo')): - return xr.DataArray(refl, dims=('y', 'x'), - attrs={ - 'start_orbit': 1708, 'end_orbit': 1708, 'wavelength': wavelength, 'level': None, - 'modifiers': None, 'calibration': calibration, 'file_key': file_key, - 'resolution': 371, 'file_type': file_type, 'name': name, - 'standard_name': standard_name, 'platform_name': 'Suomi-NPP', - 'polarization': None, 'sensor': 'viirs', 'units': units, - 'start_time': datetime.datetime(2012, 2, 25, 18, 1, 24, 570942), - 'end_time': datetime.datetime(2012, 2, 25, 18, 11, 21, 175760), 'area': area, - 'ancillary_variables': [] - }) - - c01 = make_xarray(None, 'I01', 'toa_bidirectional_reflectance', wavelength=(0.6, 0.64, 0.68), units='%', - calibration='reflectance', file_type='svi01') - c02 = make_xarray('All_Data/{file_group}_All/SatelliteAzimuthAngle', 'satellite_azimuth_angle', - 'sensor_azimuth_angle') - c03 = make_xarray('All_Data/{file_group}_All/SatelliteZenithAngle', 'satellite_zenith_angle', - 'sensor_zenith_angle') - c04 = make_xarray('All_Data/{file_group}_All/SolarAzimuthAngle', 'solar_azimuth_angle', - 'solar_azimuth_angle') - c05 = make_xarray('All_Data/{file_group}_All/SolarZenithAngle', 'solar_zenith_angle', - 'solar_zenith_angle') - with mock.patch('satpy.composites.viirs.retrieve') as rmock: - dem_fn = str(tmpdir.join('CMGDEM.hdf')) - rmock.return_value = dem_fn - from netCDF4 import Dataset - nc = Dataset(dem_fn, 'w') - nc.createDimension('y', 10) - nc.createDimension('x', 10) - dem_var = nc.createVariable('averaged elevation', np.float32, - ('y', 'x')) - dem_var[:] = 0 - res = ref_cor([c01], [c02, c03, c04, c05]) - - assert isinstance(res, xr.DataArray) - assert isinstance(res.data, da.Array) - assert res.attrs['wavelength'] == (0.6, 0.64, 0.68) - assert res.attrs['modifiers'] == ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband') - assert res.attrs['calibration'] == 'reflectance' - assert res.attrs['resolution'] == 371 - assert res.attrs['file_type'] == 'svi01' - assert res.attrs['name'] == 'I01' - assert res.attrs['standard_name'] == 'toa_bidirectional_reflectance' - assert res.attrs['platform_name'] == 'Suomi-NPP' - assert res.attrs['sensor'] == 'viirs' - assert res.attrs['units'] == '%' - assert res.attrs['start_time'] == datetime.datetime(2012, 2, 25, 18, 1, 24, 570942) - assert res.attrs['end_time'] == datetime.datetime(2012, 2, 25, 18, 11, 21, 175760) - assert res.attrs['area'] == area - assert res.attrs['ancillary_variables'] == [] - data = res.values - assert abs(np.mean(data) - 40.7578684169142) < 1e-10 - assert data.shape == (5, 10) - unique = np.unique(data) - np.testing.assert_allclose(unique, [25.20341702519979, 52.38819447051263, 75.79089653845898]) + def _start_dem_mock(self, tmpdir, url): + if not url: + return + rmock_obj = mock.patch('satpy.composites.viirs.retrieve') + rmock = rmock_obj.start() + dem_fn = str(tmpdir.join(url)) + rmock.return_value = dem_fn + from netCDF4 import Dataset + + nc = Dataset(dem_fn, 'w') + nc.createDimension('y', 10) + nc.createDimension('x', 10) + dem_var = nc.createVariable('averaged elevation', np.float32, + ('y', 'x')) + dem_var[:] = 0 + return rmock_obj + + def _stop_dem_mock(self, rmock_obj): + if rmock_obj: + rmock_obj.stop() From 71fb408f04bdfb959bf0a04b063d0da6e4ef6534 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 25 Mar 2021 14:49:19 -0500 Subject: [PATCH 6/7] Move ReflectanceCorrector to satpy.modifiers --- satpy/composites/viirs.py | 138 -------- satpy/etc/composites/abi.yaml | 2 +- satpy/etc/composites/modis.yaml | 2 +- satpy/etc/composites/viirs.yaml | 4 +- satpy/modifiers/_crefl.py | 164 ++++++++++ satpy/modifiers/atmosphere.py | 1 + satpy/tests/compositor_tests/test_viirs.py | 348 -------------------- satpy/tests/modifier_tests/__init__.py | 15 + satpy/tests/modifier_tests/test_crefl.py | 364 +++++++++++++++++++++ 9 files changed, 548 insertions(+), 490 deletions(-) create mode 100644 satpy/modifiers/_crefl.py create mode 100644 satpy/tests/modifier_tests/__init__.py create mode 100644 satpy/tests/modifier_tests/test_crefl.py diff --git a/satpy/composites/viirs.py b/satpy/composites/viirs.py index c467e1e454..576e27d16a 100644 --- a/satpy/composites/viirs.py +++ b/satpy/composites/viirs.py @@ -18,7 +18,6 @@ """Composite classes for the VIIRS instrument.""" import logging -import warnings import numpy as np import dask @@ -26,148 +25,11 @@ import xarray as xr from satpy.composites import CompositeBase, GenericCompositor -from satpy.modifiers import ModifierBase from satpy.dataset import combine_metadata -from satpy.utils import get_satpos -from satpy.aux_download import DataDownloadMixin, retrieve 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 - - class HistogramDNB(CompositeBase): """Histogram equalized DNB composite. diff --git a/satpy/etc/composites/abi.yaml b/satpy/etc/composites/abi.yaml index c792c4b906..a0cd4f18c2 100644 --- a/satpy/etc/composites/abi.yaml +++ b/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 diff --git a/satpy/etc/composites/modis.yaml b/satpy/etc/composites/modis.yaml index b4f775aa62..27524d81ee 100644 --- a/satpy/etc/composites/modis.yaml +++ b/satpy/etc/composites/modis.yaml @@ -1,7 +1,7 @@ sensor_name: visir/modis modifiers: rayleigh_corrected_crefl: - modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector + modifier: !!python/name:satpy.modifiers.atmosphere.ReflectanceCorrector url: "https://www.ssec.wisc.edu/~davidh/polar2grid/modis_crefl/tbase.hdf" known_hash: "sha256:ed5183cddce905361c1cac8ae6e3a447212875ea421a05747751efe76f8a068e" prerequisites: diff --git a/satpy/etc/composites/viirs.yaml b/satpy/etc/composites/viirs.yaml index 44fde79594..2862bc027f 100644 --- a/satpy/etc/composites/viirs.yaml +++ b/satpy/etc/composites/viirs.yaml @@ -2,7 +2,7 @@ sensor_name: visir/viirs modifiers: rayleigh_corrected_crefl: - modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector + modifier: !!python/name:satpy.modifiers.atmosphere.ReflectanceCorrector url: "https://www.ssec.wisc.edu/~davidh/polar2grid/viirs_crefl/CMGDEM.hdf" known_hash: "sha256:f33f1f867d79fff4fafe128f61c154236dd74fcc97bf418ea1437977a38d0604" prerequisites: @@ -16,7 +16,7 @@ modifiers: resolution: 742 rayleigh_corrected_crefl_iband: - modifier: !!python/name:satpy.composites.viirs.ReflectanceCorrector + modifier: !!python/name:satpy.modifiers.atmosphere.ReflectanceCorrector url: "https://www.ssec.wisc.edu/~davidh/polar2grid/viirs_crefl/CMGDEM.hdf" known_hash: "sha256:f33f1f867d79fff4fafe128f61c154236dd74fcc97bf418ea1437977a38d0604" prerequisites: diff --git a/satpy/modifiers/_crefl.py b/satpy/modifiers/_crefl.py new file mode 100644 index 0000000000..340f0ab0f2 --- /dev/null +++ b/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 . +"""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 diff --git a/satpy/modifiers/atmosphere.py b/satpy/modifiers/atmosphere.py index fb152f0b91..d49949508d 100644 --- a/satpy/modifiers/atmosphere.py +++ b/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__) diff --git a/satpy/tests/compositor_tests/test_viirs.py b/satpy/tests/compositor_tests/test_viirs.py index f2bb724aaf..a6604ce3b4 100644 --- a/satpy/tests/compositor_tests/test_viirs.py +++ b/satpy/tests/compositor_tests/test_viirs.py @@ -18,12 +18,6 @@ """Tests for VIIRS compositors.""" import unittest -from unittest import mock - -import pytest -import dask.array as da -import numpy as np -from pyresample.geometry import AreaDefinition class TestVIIRSComposites(unittest.TestCase): @@ -230,345 +224,3 @@ def test_hncc_dnb(self): unique, [3.48479712e-04, 6.96955799e-04, 1.04543189e-03, 4.75394738e-03, 9.50784532e-03, 1.42617433e-02, 1.50001560e+03, 3.00001560e+03, 4.50001560e+03]) - - -class TestViirsReflectanceCorrectorAnglesTest(unittest.TestCase): - """Tests for the VIIRS/MODIS Corrected Reflectance modifier handling angles.""" - - def setUp(self): - """Patch in-class imports.""" - self.astronomy = mock.MagicMock() - self.orbital = mock.MagicMock() - modules = { - 'pyorbital.astronomy': self.astronomy, - 'pyorbital.orbital': self.orbital, - } - self.module_patcher = mock.patch.dict('sys.modules', modules) - self.module_patcher.start() - - def tearDown(self): - """Unpatch in-class imports.""" - self.module_patcher.stop() - - @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() - 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'}) - - # Compute angles - psp = ReflectanceCorrector(name='dummy') - psp.get_angles(vis) - - # Check arguments of get_orbserver_look() call, especially the altitude - # unit conversion from meters to kilometers - 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) - - -class TestReflectanceCorrectorModifier: - """Test the CREFL modifier.""" - - @staticmethod - def data_area_ref_corrector(): - """Create test area definition and data.""" - rows = 5 - cols = 10 - area = AreaDefinition( - 'some_area_name', 'On-the-fly area', 'geosabii', - {'a': '6378137.0', 'b': '6356752.31414', 'h': '35786023.0', 'lon_0': '-89.5', 'proj': 'geos', 'sweep': 'x', - 'units': 'm'}, - cols, rows, - (-5434894.954752679, -5434894.964451744, 5434894.964451744, 5434894.954752679)) - - dnb = np.zeros((rows, cols)) + 25 - dnb[3, :] += 25 - dnb[4:, :] += 50 - dnb = da.from_array(dnb, chunks=100) - return area, dnb - - def test_reflectance_corrector_abi(self): - """Test ReflectanceCorrector modifier with ABI data.""" - import xarray as xr - import dask.array as da - import numpy as np - from satpy.composites.viirs import ReflectanceCorrector - from satpy.tests.utils import make_dsq - ref_cor = ReflectanceCorrector(optional_prerequisites=[ - make_dsq(name='satellite_azimuth_angle'), - make_dsq(name='satellite_zenith_angle'), - make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle')], name='C01', prerequisites=[], - wavelength=(0.45, 0.47, 0.49), resolution=1000, calibration='reflectance', - modifiers=('sunz_corrected', 'rayleigh_corrected_crefl',), sensor='abi') - - assert ref_cor.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl') - assert ref_cor.attrs['calibration'] == 'reflectance' - assert ref_cor.attrs['wavelength'] == (0.45, 0.47, 0.49) - assert ref_cor.attrs['name'] == 'C01' - assert ref_cor.attrs['resolution'] == 1000 - assert ref_cor.attrs['sensor'] == 'abi' - assert ref_cor.attrs['prerequisites'] == [] - assert ref_cor.attrs['optional_prerequisites'] == [ - make_dsq(name='satellite_azimuth_angle'), - make_dsq(name='satellite_zenith_angle'), - make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle')] - - area, dnb = self.data_area_ref_corrector() - c01 = xr.DataArray(dnb, - dims=('y', 'x'), - attrs={ - 'satellite_longitude': -89.5, 'satellite_latitude': 0.0, - 'satellite_altitude': 35786023.4375, 'platform_name': 'GOES-16', - 'calibration': 'reflectance', 'units': '%', 'wavelength': (0.45, 0.47, 0.49), - 'name': 'C01', 'resolution': 1000, 'sensor': 'abi', - 'start_time': '2017-09-20 17:30:40.800000', 'end_time': '2017-09-20 17:41:17.500000', - 'area': area, 'ancillary_variables': [] - }) - res = ref_cor([c01], []) - - assert isinstance(res, xr.DataArray) - assert isinstance(res.data, da.Array) - assert res.attrs['satellite_longitude'] == -89.5 - assert res.attrs['satellite_latitude'] == 0.0 - assert res.attrs['satellite_altitude'] == 35786023.4375 - assert res.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl') - assert res.attrs['platform_name'] == 'GOES-16' - assert res.attrs['calibration'] == 'reflectance' - assert res.attrs['units'] == '%' - assert res.attrs['wavelength'] == (0.45, 0.47, 0.49) - assert res.attrs['name'] == 'C01' - assert res.attrs['resolution'] == 1000 - assert res.attrs['sensor'] == 'abi' - assert res.attrs['start_time'] == '2017-09-20 17:30:40.800000' - assert res.attrs['end_time'] == '2017-09-20 17:41:17.500000' - assert res.attrs['area'] == area - assert res.attrs['ancillary_variables'] == [] - data = res.values - assert abs(np.nanmean(data) - 26.00760944144745) < 1e-10 - assert data.shape == (5, 10) - 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, - 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]) - - @pytest.mark.parametrize('url', [None, 'CMGDEM.hdf']) - def test_reflectance_corrector_viirs(self, tmpdir, url): - """Test ReflectanceCorrector modifier with VIIRS data.""" - import xarray as xr - import dask.array as da - import numpy as np - import datetime - from satpy.composites.viirs import ReflectanceCorrector - from satpy.tests.utils import make_dsq - - ref_cor = ReflectanceCorrector( - optional_prerequisites=[ - make_dsq(name='satellite_azimuth_angle'), - make_dsq(name='satellite_zenith_angle'), - make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle') - ], - name='I01', - prerequisites=[], - wavelength=(0.6, 0.64, 0.68), - resolution=371, - calibration='reflectance', - modifiers=('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband'), - sensor='viirs', - url=url) - - assert ref_cor.attrs['modifiers'] == ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband') - assert ref_cor.attrs['calibration'] == 'reflectance' - assert ref_cor.attrs['wavelength'] == (0.6, 0.64, 0.68) - assert ref_cor.attrs['name'] == 'I01' - assert ref_cor.attrs['resolution'] == 371 - assert ref_cor.attrs['sensor'] == 'viirs' - assert ref_cor.attrs['prerequisites'] == [] - assert ref_cor.attrs['optional_prerequisites'] == [ - make_dsq(name='satellite_azimuth_angle'), - make_dsq(name='satellite_zenith_angle'), - make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle')] - - area, dnb = self.data_area_ref_corrector() - - def make_xarray(name, standard_name, wavelength=None, units='degrees', calibration=None): - return xr.DataArray(dnb, dims=('y', 'x'), - attrs={ - 'start_orbit': 1708, 'end_orbit': 1708, 'wavelength': wavelength, 'level': None, - 'modifiers': None, 'calibration': calibration, - 'resolution': 371, 'name': name, - 'standard_name': standard_name, 'platform_name': 'Suomi-NPP', - 'polarization': None, 'sensor': 'viirs', 'units': units, - 'start_time': datetime.datetime(2012, 2, 25, 18, 1, 24, 570942), - 'end_time': datetime.datetime(2012, 2, 25, 18, 11, 21, 175760), 'area': area, - 'ancillary_variables': [] - }) - - c01 = make_xarray('I01', 'toa_bidirectional_reflectance', - wavelength=(0.6, 0.64, 0.68), units='%', - calibration='reflectance') - c02 = make_xarray('satellite_azimuth_angle', 'sensor_azimuth_angle') - c03 = make_xarray('satellite_zenith_angle', 'sensor_zenith_angle') - c04 = make_xarray('solar_azimuth_angle', 'solar_azimuth_angle') - c05 = make_xarray('solar_zenith_angle', 'solar_zenith_angle') - - rmock_obj = self._start_dem_mock(tmpdir, url) - res = ref_cor([c01], [c02, c03, c04, c05]) - self._stop_dem_mock(rmock_obj) - - assert isinstance(res, xr.DataArray) - assert isinstance(res.data, da.Array) - assert res.attrs['wavelength'] == (0.6, 0.64, 0.68) - assert res.attrs['modifiers'] == ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband') - assert res.attrs['calibration'] == 'reflectance' - assert res.attrs['resolution'] == 371 - assert res.attrs['name'] == 'I01' - assert res.attrs['standard_name'] == 'toa_bidirectional_reflectance' - assert res.attrs['platform_name'] == 'Suomi-NPP' - assert res.attrs['sensor'] == 'viirs' - assert res.attrs['units'] == '%' - assert res.attrs['start_time'] == datetime.datetime(2012, 2, 25, 18, 1, 24, 570942) - assert res.attrs['end_time'] == datetime.datetime(2012, 2, 25, 18, 11, 21, 175760) - assert res.attrs['area'] == area - assert res.attrs['ancillary_variables'] == [] - data = res.values - assert abs(np.mean(data) - 40.7578684169142) < 1e-10 - assert data.shape == (5, 10) - unique = np.unique(data) - np.testing.assert_allclose(unique, [25.20341702519979, 52.38819447051263, 75.79089653845898]) - - def test_reflectance_corrector_modis(self): - """Test ReflectanceCorrector modifier with MODIS data.""" - import xarray as xr - import dask.array as da - import numpy as np - import datetime - from satpy.composites.viirs import ReflectanceCorrector - from satpy.tests.utils import make_dsq - sataa_did = make_dsq(name='satellite_azimuth_angle') - satza_did = make_dsq(name='satellite_zenith_angle') - solaa_did = make_dsq(name='solar_azimuth_angle') - solza_did = make_dsq(name='solar_zenith_angle') - ref_cor = ReflectanceCorrector( - optional_prerequisites=[sataa_did, satza_did, solaa_did, solza_did], name='1', - prerequisites=[], wavelength=(0.62, 0.645, 0.67), resolution=250, calibration='reflectance', - modifiers=('sunz_corrected', 'rayleigh_corrected_crefl'), sensor='modis') - assert ref_cor.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl') - assert ref_cor.attrs['calibration'] == 'reflectance' - assert ref_cor.attrs['wavelength'] == (0.62, 0.645, 0.67) - assert ref_cor.attrs['name'] == '1' - assert ref_cor.attrs['resolution'] == 250 - assert ref_cor.attrs['sensor'] == 'modis' - assert ref_cor.attrs['prerequisites'] == [] - assert ref_cor.attrs['optional_prerequisites'] == [ - make_dsq(name='satellite_azimuth_angle'), - make_dsq(name='satellite_zenith_angle'), - make_dsq(name='solar_azimuth_angle'), - make_dsq(name='solar_zenith_angle')] - - area, dnb = self.data_area_ref_corrector() - - def make_xarray(name, calibration, wavelength=None, modifiers=None, resolution=1000, - file_type='hdf_eos_geo'): - return xr.DataArray(dnb, - dims=('y', 'x'), - attrs={ - 'wavelength': wavelength, 'level': None, 'modifiers': modifiers, - 'calibration': calibration, 'resolution': resolution, 'file_type': file_type, - 'name': name, 'coordinates': ['longitude', 'latitude'], - 'platform_name': 'EOS-Aqua', 'polarization': None, 'sensor': 'modis', - 'units': '%', 'start_time': datetime.datetime(2012, 8, 13, 18, 46, 1, 439838), - 'end_time': datetime.datetime(2012, 8, 13, 18, 57, 47, 746296), 'area': area, - 'ancillary_variables': [] - }) - - c01 = make_xarray('1', 'reflectance', wavelength=(0.62, 0.645, 0.67), modifiers='sunz_corrected', - resolution=500, file_type='hdf_eos_data_500m') - c02 = make_xarray('satellite_azimuth_angle', None) - c03 = make_xarray('satellite_zenith_angle', None) - c04 = make_xarray('solar_azimuth_angle', None) - c05 = make_xarray('solar_zenith_angle', None) - res = ref_cor([c01], [c02, c03, c04, c05]) - - assert isinstance(res, xr.DataArray) - assert isinstance(res.data, da.Array) - assert res.attrs['wavelength'] == (0.62, 0.645, 0.67) - assert res.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl',) - assert res.attrs['calibration'] == 'reflectance' - assert res.attrs['resolution'] == 500 - assert res.attrs['file_type'] == 'hdf_eos_data_500m' - assert res.attrs['name'] == '1' - assert res.attrs['platform_name'] == 'EOS-Aqua' - assert res.attrs['sensor'] == 'modis' - assert res.attrs['units'] == '%' - assert res.attrs['start_time'] == datetime.datetime(2012, 8, 13, 18, 46, 1, 439838) - assert res.attrs['end_time'] == datetime.datetime(2012, 8, 13, 18, 57, 47, 746296) - assert res.attrs['area'] == area - assert res.attrs['ancillary_variables'] == [] - data = res.values - if abs(np.mean(data) - 38.734365117099145) >= 1e-10: - raise AssertionError('{} is not within {} of {}'.format(np.mean(data), 1e-10, 38.734365117099145)) - assert data.shape == (5, 10) - unique = np.unique(data) - np.testing.assert_allclose(unique, [24.641586, 50.431692, 69.315375]) - - def test_reflectance_corrector_bad_prereqs(self): - """Test ReflectanceCorrector modifier with wrong number of inputs.""" - from satpy.composites.viirs import ReflectanceCorrector - ref_cor = ReflectanceCorrector("test") - pytest.raises(ValueError, ref_cor, [1], [2, 3, 4]) - pytest.raises(ValueError, ref_cor, [1, 2, 3, 4], []) - pytest.raises(ValueError, ref_cor, [], [1, 2, 3, 4]) - - def _start_dem_mock(self, tmpdir, url): - if not url: - return - rmock_obj = mock.patch('satpy.composites.viirs.retrieve') - rmock = rmock_obj.start() - dem_fn = str(tmpdir.join(url)) - rmock.return_value = dem_fn - from netCDF4 import Dataset - - nc = Dataset(dem_fn, 'w') - nc.createDimension('y', 10) - nc.createDimension('x', 10) - dem_var = nc.createVariable('averaged elevation', np.float32, - ('y', 'x')) - dem_var[:] = 0 - return rmock_obj - - def _stop_dem_mock(self, rmock_obj): - if rmock_obj: - rmock_obj.stop() diff --git a/satpy/tests/modifier_tests/__init__.py b/satpy/tests/modifier_tests/__init__.py new file mode 100644 index 0000000000..3094b02321 --- /dev/null +++ b/satpy/tests/modifier_tests/__init__.py @@ -0,0 +1,15 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2018 - 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. +"""Tests for modifiers.""" diff --git a/satpy/tests/modifier_tests/test_crefl.py b/satpy/tests/modifier_tests/test_crefl.py new file mode 100644 index 0000000000..a70d088417 --- /dev/null +++ b/satpy/tests/modifier_tests/test_crefl.py @@ -0,0 +1,364 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2018 - 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. +"""Tests for the CREFL ReflectanceCorrector modifier.""" +import unittest +from unittest import mock + +import numpy as np +import pytest +from dask import array as da +from pyresample.geometry import AreaDefinition + + +class TestViirsReflectanceCorrectorAnglesTest(unittest.TestCase): + """Tests for the VIIRS/MODIS Corrected Reflectance modifier handling angles.""" + + def setUp(self): + """Patch in-class imports.""" + self.astronomy = mock.MagicMock() + self.orbital = mock.MagicMock() + modules = { + 'pyorbital.astronomy': self.astronomy, + 'pyorbital.orbital': self.orbital, + } + self.module_patcher = mock.patch.dict('sys.modules', modules) + self.module_patcher.start() + + def tearDown(self): + """Unpatch in-class imports.""" + self.module_patcher.stop() + + @mock.patch('satpy.modifiers._crefl.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.modifiers._crefl 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() + 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'}) + + # Compute angles + psp = ReflectanceCorrector(name='dummy') + psp.get_angles(vis) + + # Check arguments of get_orbserver_look() call, especially the altitude + # unit conversion from meters to kilometers + 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) + + +class TestReflectanceCorrectorModifier: + """Test the CREFL modifier.""" + + @staticmethod + def data_area_ref_corrector(): + """Create test area definition and data.""" + rows = 5 + cols = 10 + area = AreaDefinition( + 'some_area_name', 'On-the-fly area', 'geosabii', + {'a': '6378137.0', 'b': '6356752.31414', 'h': '35786023.0', 'lon_0': '-89.5', 'proj': 'geos', 'sweep': 'x', + 'units': 'm'}, + cols, rows, + (-5434894.954752679, -5434894.964451744, 5434894.964451744, 5434894.954752679)) + + dnb = np.zeros((rows, cols)) + 25 + dnb[3, :] += 25 + dnb[4:, :] += 50 + dnb = da.from_array(dnb, chunks=100) + return area, dnb + + def test_reflectance_corrector_abi(self): + """Test ReflectanceCorrector modifier with ABI data.""" + import xarray as xr + import dask.array as da + import numpy as np + from satpy.modifiers._crefl import ReflectanceCorrector + from satpy.tests.utils import make_dsq + ref_cor = ReflectanceCorrector(optional_prerequisites=[ + make_dsq(name='satellite_azimuth_angle'), + make_dsq(name='satellite_zenith_angle'), + make_dsq(name='solar_azimuth_angle'), + make_dsq(name='solar_zenith_angle')], name='C01', prerequisites=[], + wavelength=(0.45, 0.47, 0.49), resolution=1000, calibration='reflectance', + modifiers=('sunz_corrected', 'rayleigh_corrected_crefl',), sensor='abi') + + assert ref_cor.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl') + assert ref_cor.attrs['calibration'] == 'reflectance' + assert ref_cor.attrs['wavelength'] == (0.45, 0.47, 0.49) + assert ref_cor.attrs['name'] == 'C01' + assert ref_cor.attrs['resolution'] == 1000 + assert ref_cor.attrs['sensor'] == 'abi' + assert ref_cor.attrs['prerequisites'] == [] + assert ref_cor.attrs['optional_prerequisites'] == [ + make_dsq(name='satellite_azimuth_angle'), + make_dsq(name='satellite_zenith_angle'), + make_dsq(name='solar_azimuth_angle'), + make_dsq(name='solar_zenith_angle')] + + area, dnb = self.data_area_ref_corrector() + c01 = xr.DataArray(dnb, + dims=('y', 'x'), + attrs={ + 'satellite_longitude': -89.5, 'satellite_latitude': 0.0, + 'satellite_altitude': 35786023.4375, 'platform_name': 'GOES-16', + 'calibration': 'reflectance', 'units': '%', 'wavelength': (0.45, 0.47, 0.49), + 'name': 'C01', 'resolution': 1000, 'sensor': 'abi', + 'start_time': '2017-09-20 17:30:40.800000', 'end_time': '2017-09-20 17:41:17.500000', + 'area': area, 'ancillary_variables': [] + }) + res = ref_cor([c01], []) + + assert isinstance(res, xr.DataArray) + assert isinstance(res.data, da.Array) + assert res.attrs['satellite_longitude'] == -89.5 + assert res.attrs['satellite_latitude'] == 0.0 + assert res.attrs['satellite_altitude'] == 35786023.4375 + assert res.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl') + assert res.attrs['platform_name'] == 'GOES-16' + assert res.attrs['calibration'] == 'reflectance' + assert res.attrs['units'] == '%' + assert res.attrs['wavelength'] == (0.45, 0.47, 0.49) + assert res.attrs['name'] == 'C01' + assert res.attrs['resolution'] == 1000 + assert res.attrs['sensor'] == 'abi' + assert res.attrs['start_time'] == '2017-09-20 17:30:40.800000' + assert res.attrs['end_time'] == '2017-09-20 17:41:17.500000' + assert res.attrs['area'] == area + assert res.attrs['ancillary_variables'] == [] + data = res.values + assert abs(np.nanmean(data) - 26.00760944144745) < 1e-10 + assert data.shape == (5, 10) + 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, + 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]) + + @pytest.mark.parametrize('url', [None, 'CMGDEM.hdf']) + def test_reflectance_corrector_viirs(self, tmpdir, url): + """Test ReflectanceCorrector modifier with VIIRS data.""" + import xarray as xr + import dask.array as da + import numpy as np + import datetime + from satpy.modifiers._crefl import ReflectanceCorrector + from satpy.tests.utils import make_dsq + + ref_cor = ReflectanceCorrector( + optional_prerequisites=[ + make_dsq(name='satellite_azimuth_angle'), + make_dsq(name='satellite_zenith_angle'), + make_dsq(name='solar_azimuth_angle'), + make_dsq(name='solar_zenith_angle') + ], + name='I01', + prerequisites=[], + wavelength=(0.6, 0.64, 0.68), + resolution=371, + calibration='reflectance', + modifiers=('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband'), + sensor='viirs', + url=url) + + assert ref_cor.attrs['modifiers'] == ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband') + assert ref_cor.attrs['calibration'] == 'reflectance' + assert ref_cor.attrs['wavelength'] == (0.6, 0.64, 0.68) + assert ref_cor.attrs['name'] == 'I01' + assert ref_cor.attrs['resolution'] == 371 + assert ref_cor.attrs['sensor'] == 'viirs' + assert ref_cor.attrs['prerequisites'] == [] + assert ref_cor.attrs['optional_prerequisites'] == [ + make_dsq(name='satellite_azimuth_angle'), + make_dsq(name='satellite_zenith_angle'), + make_dsq(name='solar_azimuth_angle'), + make_dsq(name='solar_zenith_angle')] + + area, dnb = self.data_area_ref_corrector() + + def make_xarray(name, standard_name, wavelength=None, units='degrees', calibration=None): + return xr.DataArray(dnb, dims=('y', 'x'), + attrs={ + 'start_orbit': 1708, 'end_orbit': 1708, 'wavelength': wavelength, 'level': None, + 'modifiers': None, 'calibration': calibration, + 'resolution': 371, 'name': name, + 'standard_name': standard_name, 'platform_name': 'Suomi-NPP', + 'polarization': None, 'sensor': 'viirs', 'units': units, + 'start_time': datetime.datetime(2012, 2, 25, 18, 1, 24, 570942), + 'end_time': datetime.datetime(2012, 2, 25, 18, 11, 21, 175760), 'area': area, + 'ancillary_variables': [] + }) + + c01 = make_xarray('I01', 'toa_bidirectional_reflectance', + wavelength=(0.6, 0.64, 0.68), units='%', + calibration='reflectance') + c02 = make_xarray('satellite_azimuth_angle', 'sensor_azimuth_angle') + c03 = make_xarray('satellite_zenith_angle', 'sensor_zenith_angle') + c04 = make_xarray('solar_azimuth_angle', 'solar_azimuth_angle') + c05 = make_xarray('solar_zenith_angle', 'solar_zenith_angle') + + rmock_obj = self._start_dem_mock(tmpdir, url) + res = ref_cor([c01], [c02, c03, c04, c05]) + self._stop_dem_mock(rmock_obj) + + assert isinstance(res, xr.DataArray) + assert isinstance(res.data, da.Array) + assert res.attrs['wavelength'] == (0.6, 0.64, 0.68) + assert res.attrs['modifiers'] == ('sunz_corrected_iband', 'rayleigh_corrected_crefl_iband') + assert res.attrs['calibration'] == 'reflectance' + assert res.attrs['resolution'] == 371 + assert res.attrs['name'] == 'I01' + assert res.attrs['standard_name'] == 'toa_bidirectional_reflectance' + assert res.attrs['platform_name'] == 'Suomi-NPP' + assert res.attrs['sensor'] == 'viirs' + assert res.attrs['units'] == '%' + assert res.attrs['start_time'] == datetime.datetime(2012, 2, 25, 18, 1, 24, 570942) + assert res.attrs['end_time'] == datetime.datetime(2012, 2, 25, 18, 11, 21, 175760) + assert res.attrs['area'] == area + assert res.attrs['ancillary_variables'] == [] + data = res.values + assert abs(np.mean(data) - 40.7578684169142) < 1e-10 + assert data.shape == (5, 10) + unique = np.unique(data) + np.testing.assert_allclose(unique, [25.20341702519979, 52.38819447051263, 75.79089653845898]) + + def test_reflectance_corrector_modis(self): + """Test ReflectanceCorrector modifier with MODIS data.""" + import xarray as xr + import dask.array as da + import numpy as np + import datetime + from satpy.modifiers._crefl import ReflectanceCorrector + from satpy.tests.utils import make_dsq + sataa_did = make_dsq(name='satellite_azimuth_angle') + satza_did = make_dsq(name='satellite_zenith_angle') + solaa_did = make_dsq(name='solar_azimuth_angle') + solza_did = make_dsq(name='solar_zenith_angle') + ref_cor = ReflectanceCorrector( + optional_prerequisites=[sataa_did, satza_did, solaa_did, solza_did], name='1', + prerequisites=[], wavelength=(0.62, 0.645, 0.67), resolution=250, calibration='reflectance', + modifiers=('sunz_corrected', 'rayleigh_corrected_crefl'), sensor='modis') + assert ref_cor.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl') + assert ref_cor.attrs['calibration'] == 'reflectance' + assert ref_cor.attrs['wavelength'] == (0.62, 0.645, 0.67) + assert ref_cor.attrs['name'] == '1' + assert ref_cor.attrs['resolution'] == 250 + assert ref_cor.attrs['sensor'] == 'modis' + assert ref_cor.attrs['prerequisites'] == [] + assert ref_cor.attrs['optional_prerequisites'] == [ + make_dsq(name='satellite_azimuth_angle'), + make_dsq(name='satellite_zenith_angle'), + make_dsq(name='solar_azimuth_angle'), + make_dsq(name='solar_zenith_angle')] + + area, dnb = self.data_area_ref_corrector() + + def make_xarray(name, calibration, wavelength=None, modifiers=None, resolution=1000, + file_type='hdf_eos_geo'): + return xr.DataArray(dnb, + dims=('y', 'x'), + attrs={ + 'wavelength': wavelength, 'level': None, 'modifiers': modifiers, + 'calibration': calibration, 'resolution': resolution, 'file_type': file_type, + 'name': name, 'coordinates': ['longitude', 'latitude'], + 'platform_name': 'EOS-Aqua', 'polarization': None, 'sensor': 'modis', + 'units': '%', 'start_time': datetime.datetime(2012, 8, 13, 18, 46, 1, 439838), + 'end_time': datetime.datetime(2012, 8, 13, 18, 57, 47, 746296), 'area': area, + 'ancillary_variables': [] + }) + + c01 = make_xarray('1', 'reflectance', wavelength=(0.62, 0.645, 0.67), modifiers='sunz_corrected', + resolution=500, file_type='hdf_eos_data_500m') + c02 = make_xarray('satellite_azimuth_angle', None) + c03 = make_xarray('satellite_zenith_angle', None) + c04 = make_xarray('solar_azimuth_angle', None) + c05 = make_xarray('solar_zenith_angle', None) + res = ref_cor([c01], [c02, c03, c04, c05]) + + assert isinstance(res, xr.DataArray) + assert isinstance(res.data, da.Array) + assert res.attrs['wavelength'] == (0.62, 0.645, 0.67) + assert res.attrs['modifiers'] == ('sunz_corrected', 'rayleigh_corrected_crefl',) + assert res.attrs['calibration'] == 'reflectance' + assert res.attrs['resolution'] == 500 + assert res.attrs['file_type'] == 'hdf_eos_data_500m' + assert res.attrs['name'] == '1' + assert res.attrs['platform_name'] == 'EOS-Aqua' + assert res.attrs['sensor'] == 'modis' + assert res.attrs['units'] == '%' + assert res.attrs['start_time'] == datetime.datetime(2012, 8, 13, 18, 46, 1, 439838) + assert res.attrs['end_time'] == datetime.datetime(2012, 8, 13, 18, 57, 47, 746296) + assert res.attrs['area'] == area + assert res.attrs['ancillary_variables'] == [] + data = res.values + if abs(np.mean(data) - 38.734365117099145) >= 1e-10: + raise AssertionError('{} is not within {} of {}'.format(np.mean(data), 1e-10, 38.734365117099145)) + assert data.shape == (5, 10) + unique = np.unique(data) + np.testing.assert_allclose(unique, [24.641586, 50.431692, 69.315375]) + + def test_reflectance_corrector_bad_prereqs(self): + """Test ReflectanceCorrector modifier with wrong number of inputs.""" + from satpy.modifiers._crefl import ReflectanceCorrector + ref_cor = ReflectanceCorrector("test") + pytest.raises(ValueError, ref_cor, [1], [2, 3, 4]) + pytest.raises(ValueError, ref_cor, [1, 2, 3, 4], []) + pytest.raises(ValueError, ref_cor, [], [1, 2, 3, 4]) + + def _start_dem_mock(self, tmpdir, url): + if not url: + return + rmock_obj = mock.patch('satpy.modifiers._crefl.retrieve') + rmock = rmock_obj.start() + dem_fn = str(tmpdir.join(url)) + rmock.return_value = dem_fn + from netCDF4 import Dataset + + nc = Dataset(dem_fn, 'w') + nc.createDimension('y', 10) + nc.createDimension('x', 10) + dem_var = nc.createVariable('averaged elevation', np.float32, + ('y', 'x')) + dem_var[:] = 0 + return rmock_obj + + def _stop_dem_mock(self, rmock_obj): + if rmock_obj: + rmock_obj.stop() From be132b876e0cbfc7b01023ed8b8dfe0b40b155c0 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Thu, 25 Mar 2021 16:39:05 -0500 Subject: [PATCH 7/7] Refactor ReflectanceCorrector modifier for cleaner code --- satpy/modifiers/_crefl.py | 105 +++++++++++++++++++++----------------- 1 file changed, 57 insertions(+), 48 deletions(-) diff --git a/satpy/modifiers/_crefl.py b/satpy/modifiers/_crefl.py index 340f0ab0f2..0dfe2f5b6f 100644 --- a/satpy/modifiers/_crefl.py +++ b/satpy/modifiers/_crefl.py @@ -36,7 +36,8 @@ class ReflectanceCorrector(ModifierBase, DataDownloadMixin): Uses a python rewrite of the C CREFL code written for VIIRS and MODIS. """ - def __init__(self, *args, **kwargs): + def __init__(self, *args, dem_filename=None, dem_sds="averaged elevation", + url=None, known_hash=None, **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 @@ -52,52 +53,54 @@ def __init__(self, *args, **kwargs): 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] + self.dem_sds = dem_sds + self.url = url + self.known_hash = known_hash + self.dem_cache_key = self._get_registered_dem_cache_key() + + def _get_registered_dem_cache_key(self): + if not self.url: + return + reg_files = self.register_data_files([{ + 'url': self.url, 'known_hash': self.known_hash} + ]) + return 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"] == "%" + from satpy.composites.crefl_utils import get_coefficients + refl_data, *angles = self._get_data_and_angles(datasets, optional_datasets) coefficients = get_coefficients(refl_data.attrs["sensor"], refl_data.attrs["wavelength"], refl_data.attrs["resolution"]) - use_abi = refl_data.attrs['sensor'] == 'abi' + results = self._call_crefl(refl_data, coefficients, angles) + info.update(refl_data.attrs) + info["rayleigh_corrected"] = True + results.attrs = info + self.apply_modifier_info(refl_data, results) + return results + + def _call_crefl(self, refl_data, coefficients, angles): + from satpy.composites.crefl_utils import run_crefl + avg_elevation = self._get_average_elevation() lons, lats = refl_data.attrs['area'].get_lonlats(chunks=refl_data.chunks) + is_percent = refl_data.attrs["units"] == "%" + use_abi = refl_data.attrs['sensor'] == 'abi' results = run_crefl(refl_data, coefficients, lons, lats, - sensor_aa, - sensor_za, - solar_aa, - solar_za, + *angles, 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): @@ -118,41 +121,47 @@ def _get_average_elevation(self): return avg_elevation def _get_data_and_angles(self, datasets, optional_datasets): + angles = self._extract_angle_data_arrays(datasets, optional_datasets) + angles = [xr.DataArray(dask_arr, dims=('y', 'x')) for dask_arr in angles] + return [datasets[0]] + angles + + def _extract_angle_data_arrays(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( + return self.get_angles(vis) + if len(all_datasets) == 5: + vis, *angles = 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 + return [data_arr.data for data_arr in angles] + raise ValueError("Not sure how to handle provided dependencies. " + "Either all 4 angles must be provided or none of " + "of them.") 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 = self._get_valid_lonlats(vis) + sun_angles = self._get_sun_angles(vis, lons, lats) + sat_angles = self._get_sensor_angles(vis, lons, lats) + # sata, satz, suna, sunz + return sat_angles + sun_angles + + def _get_valid_lonlats(self, vis): 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) + return lons, lats + + def _get_sun_angles(self, vis, lons, lats): + from pyorbital.astronomy import get_alt_az, sun_zenith_angle 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) + return suna, sunz + + def _get_sensor_angles(self, vis, lons, lats): + from pyorbital.orbital import get_observer_look sat_lon, sat_lat, sat_alt = get_satpos(vis) sata, satel = get_observer_look( sat_lon, @@ -161,4 +170,4 @@ def get_angles(self, vis): vis.attrs['start_time'], lons, lats, 0) satz = 90 - satel - return sata, satz, suna, sunz + return sata, satz