Merge pull request #2175 from simonrp84/raa_calc

satpy/modifiers/angles.py

@@ -317,6 +317,25 @@ def _geo_dask_to_data_array(arr: da.Array) -> xr.DataArray:
     return xr.DataArray(arr, dims=('y', 'x'))
 
 
+def compute_relative_azimuth(sat_azi: xr.DataArray, sun_azi: xr.DataArray) -> xr.DataArray:
+    """Compute the relative azimuth angle.
+
+    Args:
+        sat_azi: DataArray for the satellite azimuth angles, typically in 0-360 degree range.
+        sun_azi: DataArray for the solar azimuth angles, should be in same range as sat_azi.
+    Returns:
+        A DataArray containing the relative azimuth angle in the 0-180 degree range.
+
+    NOTE: Relative azimuth is defined such that:
+    Relative azimuth is 0 when sun and satellite are aligned on one side of a pixel (back scatter).
+    Relative azimuth is 180 when sun and satellite are directly opposite each other (forward scatter).
+    """
+    ssadiff = np.absolute(sun_azi - sat_azi)
+    ssadiff = np.minimum(ssadiff, 360 - ssadiff)
+
+    return ssadiff
+
+
 def get_angles(data_arr: xr.DataArray) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray]:
     """Get sun and satellite azimuth and zenith angles.
 

satpy/modifiers/atmosphere.py

@@ -25,7 +25,7 @@
 
 from satpy.modifiers import ModifierBase
 from satpy.modifiers._crefl import ReflectanceCorrector  # noqa
-from satpy.modifiers.angles import get_angles, get_satellite_zenith_angle
+from satpy.modifiers.angles import compute_relative_azimuth, get_angles, get_satellite_zenith_angle
 
 logger = logging.getLogger(__name__)
 
@@ -45,18 +45,17 @@ def __call__(self, projectables, optional_datasets=None, **info):
         else:
             vis, red, sata, satz, suna, sunz = self.match_data_arrays(
                 projectables + optional_datasets)
+            # First make sure the two azimuth angles are in the range 0-360:
+            sata = sata % 360.
+            suna = suna % 360.
 
         # get the dask array underneath
         sata = sata.data
         satz = satz.data
         suna = suna.data
         sunz = sunz.data
 
-        # First make sure the two azimuth angles are in the range 0-360:
-        sata = sata % 360.
-        suna = suna % 360.
-        ssadiff = da.absolute(suna - sata)
-        ssadiff = da.minimum(ssadiff, 360 - ssadiff)
+        ssadiff = compute_relative_azimuth(sata, suna)
         del sata, suna
 
         atmosphere = self.attrs.get('atmosphere', 'us-standard')

satpy/tests/modifier_tests/test_angles.py

@@ -324,6 +324,19 @@ def test_no_cache_dir_fails(self, tmp_path):
                 satpy.config.set(cache_lonlats=True, cache_sensor_angles=True, cache_dir=None):
             _get_sensor_angles_from_sat_pos.cache_clear()
 
+    def test_relative_azimuth_calculation(self):
+        """Test relative azimuth calculation."""
+        from satpy.modifiers.angles import compute_relative_azimuth
+
+        saa = xr.DataArray(np.array([-120, 40., 0.04, 179.4, 94.2, 12.1]))
+        vaa = xr.DataArray(np.array([60., 57.7, 175.1, 234.18, 355.4, 12.1]))
+
+        expected_raa = xr.DataArray(np.array([180., 17.7, 175.06, 54.78, 98.8, 0.]))
+
+        raa = compute_relative_azimuth(vaa, saa)
+        assert isinstance(raa, xr.DataArray)
+        np.testing.assert_allclose(expected_raa, raa)
+
     def test_solazi_correction(self):
         """Test that solar azimuth angles are corrected into the right range."""
         from datetime import datetime