Skip to content

Commit

Permalink
Fix rayleigh correction not handling angles as required inputs
Browse files Browse the repository at this point in the history
  • Loading branch information
@djhoese
djhoese committed May 10, 2023
1 parent 2aac9e4 commit 5e2cb88
Show file tree
Hide file tree
Showing 2 changed files with 82 additions and 37 deletions.
6 changes: 3 additions & 3 deletions satpy/modifiers/atmosphere.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,12 +74,12 @@ def __call__(self, projectables, optional_datasets=None, **info):
Uses pyspectral.
"""
from pyspectral.rayleigh import Rayleigh
if not optional_datasets or len(optional_datasets) != 4:
projectables = projectables + (optional_datasets or [])
if len(projectables) != 6:
vis, red = self.match_data_arrays(projectables)
sata, satz, suna, sunz = get_angles(vis)
else:
vis, red, sata, satz, suna, sunz = self.match_data_arrays(
projectables + optional_datasets)
vis, red, sata, satz, suna, sunz = self.match_data_arrays(projectables)
# First make sure the two azimuth angles are in the range 0-360:
sata = sata % 360.
suna = suna % 360.
Expand Down
113 changes: 79 additions & 34 deletions satpy/tests/test_modifiers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -395,7 +395,7 @@ def test_compositor(self, calculator, apply_modifier_info, sza):
class TestPSPRayleighReflectance:
"""Test the pyspectral-based Rayleigh correction modifier."""

def make_data_area(self):
def _make_data_area(self):
"""Create test area definition and data."""
rows = 3
cols = 5
Expand All @@ -412,39 +412,8 @@ def make_data_area(self):
data = da.from_array(data, chunks=2)
return area, data

@pytest.mark.parametrize(
("name", "wavelength", "resolution", "aerosol_type", "reduce_lim_low", "reduce_lim_high", "reduce_strength",
"exp_mean", "exp_unique"),
[
("B01", (0.45, 0.47, 0.49), 1000, "rayleigh_only", 70, 95, 1, 41.540239,
np.array([9.22630464, 10.67844368, 13.58057226, 37.92186549, 40.13822472, 44.66259518,
44.92748445, 45.03917091, 69.5821722, 70.11226943, 71.07352559])),
("B02", (0.49, 0.51, 0.53), 1000, "rayleigh_only", 70, 95, 1, 43.663805,
np.array([13.15770104, 14.26526104, 16.49084485, 40.88633902, 42.60682921, 46.04288,
46.2356062, 46.28276282, 70.92799823, 71.33561614, 72.07001693])),
("B03", (0.62, 0.64, 0.66), 500, "rayleigh_only", 70, 95, 1, 46.916187,
np.array([19.22922328, 19.76884762, 20.91027446, 45.51075967, 46.39925968, 48.10221156,
48.15715058, 48.18698356, 73.01115816, 73.21552816, 73.58666477])),
("B01", (0.45, 0.47, 0.49), 1000, "rayleigh_only", -95, -70, -1, 41.540239,
np.array([9.22630464, 10.67844368, 13.58057226, 37.92186549, 40.13822472, 44.66259518,
44.92748445, 45.03917091, 69.5821722, 70.11226943, 71.07352559])),
]
)
def test_rayleigh_corrector(self, name, wavelength, resolution, aerosol_type, reduce_lim_low, reduce_lim_high,
reduce_strength, exp_mean, exp_unique):
"""Test PSPRayleighReflectance with fake data."""
from satpy.modifiers.atmosphere import PSPRayleighReflectance
ray_cor = PSPRayleighReflectance(name=name, atmosphere='us-standard', aerosol_types=aerosol_type,
reduce_lim_low=reduce_lim_low, reduce_lim_high=reduce_lim_high,
reduce_strength=reduce_strength)
assert ray_cor.attrs['name'] == name
assert ray_cor.attrs['atmosphere'] == 'us-standard'
assert ray_cor.attrs['aerosol_types'] == aerosol_type
assert ray_cor.attrs['reduce_lim_low'] == reduce_lim_low
assert ray_cor.attrs['reduce_lim_high'] == reduce_lim_high
assert ray_cor.attrs['reduce_strength'] == reduce_strength

area, dnb = self.make_data_area()
def _create_test_data(self, name, wavelength, resolution):
area, dnb = self._make_data_area()
input_band = xr.DataArray(dnb,
dims=('y', 'x'),
attrs={
Expand Down Expand Up @@ -476,7 +445,52 @@ def test_rayleigh_corrector(self, name, wavelength, resolution, aerosol_type, re
'satellite_nominal_altitude': 35786023.4375,
},
})
fake_angle_data = da.ones_like(dnb, dtype=np.float32) * 90.0
angle1 = xr.DataArray(fake_angle_data,
dims=('y', 'x'),
attrs={
'platform_name': 'Himawari-8',
'calibration': 'reflectance', 'units': '%', 'wavelength': wavelength,
'name': "satellite_azimuth_angle", 'resolution': resolution, 'sensor': 'ahi',
'start_time': '2017-09-20 17:30:40.800000',
'end_time': '2017-09-20 17:41:17.500000',
'area': area, 'ancillary_variables': [],
})
return input_band, red_band, angle1, angle1, angle1, angle1

@pytest.mark.parametrize(
("name", "wavelength", "resolution", "aerosol_type", "reduce_lim_low", "reduce_lim_high", "reduce_strength",
"exp_mean", "exp_unique"),
[
("B01", (0.45, 0.47, 0.49), 1000, "rayleigh_only", 70, 95, 1, 41.540239,
np.array([9.22630464, 10.67844368, 13.58057226, 37.92186549, 40.13822472, 44.66259518,
44.92748445, 45.03917091, 69.5821722, 70.11226943, 71.07352559])),
("B02", (0.49, 0.51, 0.53), 1000, "rayleigh_only", 70, 95, 1, 43.663805,
np.array([13.15770104, 14.26526104, 16.49084485, 40.88633902, 42.60682921, 46.04288,
46.2356062, 46.28276282, 70.92799823, 71.33561614, 72.07001693])),
("B03", (0.62, 0.64, 0.66), 500, "rayleigh_only", 70, 95, 1, 46.916187,
np.array([19.22922328, 19.76884762, 20.91027446, 45.51075967, 46.39925968, 48.10221156,
48.15715058, 48.18698356, 73.01115816, 73.21552816, 73.58666477])),
("B01", (0.45, 0.47, 0.49), 1000, "rayleigh_only", -95, -70, -1, 41.540239,
np.array([9.22630464, 10.67844368, 13.58057226, 37.92186549, 40.13822472, 44.66259518,
44.92748445, 45.03917091, 69.5821722, 70.11226943, 71.07352559])),
]
)
def test_rayleigh_corrector(self, name, wavelength, resolution, aerosol_type, reduce_lim_low, reduce_lim_high,
reduce_strength, exp_mean, exp_unique):
"""Test PSPRayleighReflectance with fake data."""
from satpy.modifiers.atmosphere import PSPRayleighReflectance
ray_cor = PSPRayleighReflectance(name=name, atmosphere='us-standard', aerosol_types=aerosol_type,
reduce_lim_low=reduce_lim_low, reduce_lim_high=reduce_lim_high,
reduce_strength=reduce_strength)
assert ray_cor.attrs['name'] == name
assert ray_cor.attrs['atmosphere'] == 'us-standard'
assert ray_cor.attrs['aerosol_types'] == aerosol_type
assert ray_cor.attrs['reduce_lim_low'] == reduce_lim_low
assert ray_cor.attrs['reduce_lim_high'] == reduce_lim_high
assert ray_cor.attrs['reduce_strength'] == reduce_strength

input_band, red_band, *_ = self._create_test_data(name, wavelength, resolution)
res = ray_cor([input_band, red_band])

assert isinstance(res, xr.DataArray)
Expand All @@ -488,6 +502,37 @@ def test_rayleigh_corrector(self, name, wavelength, resolution, aerosol_type, re
assert data.shape == (3, 5)
np.testing.assert_allclose(unique, exp_unique, rtol=1e-5)

@pytest.mark.parametrize("as_optionals", [False, True])
def test_rayleigh_with_angles(self, as_optionals):
"""Test PSPRayleighReflectance with angles provided."""
from satpy.modifiers.atmosphere import PSPRayleighReflectance
aerosol_type = "rayleigh_only"
ray_cor = PSPRayleighReflectance(name="B01", atmosphere='us-standard', aerosol_types=aerosol_type)
prereqs, opt_prereqs = self._get_angles_prereqs_and_opts(as_optionals)
with mock.patch("satpy.modifiers.atmosphere.get_angles") as get_angles:
res = ray_cor(prereqs, opt_prereqs)
get_angles.assert_not_called()

assert isinstance(res, xr.DataArray)
assert isinstance(res.data, da.Array)

data = res.values
unique = np.unique(data[~np.isnan(data)])
np.testing.assert_allclose(unique, np.array([-75.0, -37.71298492, 31.14350754]), rtol=1e-5)
assert data.shape == (3, 5)

def _get_angles_prereqs_and_opts(self, as_optionals):
wavelength = (0.45, 0.47, 0.49)
resolution = 1000
input_band, red_band, *angles = self._create_test_data("B01", wavelength, resolution)
prereqs = [input_band, red_band]
opt_prereqs = []
if as_optionals:
opt_prereqs = angles
else:
prereqs += angles
return prereqs, opt_prereqs


class TestPSPAtmosphericalCorrection(unittest.TestCase):
"""Test the pyspectral-based atmospheric correction modifier."""
Expand Down

0 comments on commit 5e2cb88

Please sign in to comment.