fix the flipping bug in day-night transmission flux averaging

docs/source/changelog.rst

@@ -7,8 +7,9 @@ Changelog
 0.4.18
 ------
 
-:Date: TBA
+:Date: 22 September 2023
 
+* Re-fix the bug in day-night transmission spectra averaging.
 * Add `lfric` functions in a submodule with the same name.
 * Add `esmf_regrid` as a dependency.
 * Move `calc_derived_cubes` from `pouch`, add `air_pressure()`.
@@ -41,7 +42,7 @@ Changelog
 0.4.15
 ------
 
-:Date: 22 March 2022
+:Date: 22 April 2022
 
 * Patch the calculation of the day-night average transmission flux.
 * Update citation.

src/aeolus/coord.py

@@ -933,7 +933,7 @@ def replace_z_coord(cube, model=um):
     return new_cube
 
 
-def roll_cube_0_360(cube_in, model=um):
+def roll_cube_0_360(cube_in, add_shift=0, model=um):
     """
     Take a cube spanning -180...180 degrees in longitude and roll it to 0...360 degrees.
 
@@ -943,6 +943,8 @@ def roll_cube_0_360(cube_in, model=um):
     ----------
     cube: iris.cube.Cube
         Cube with longitude and latitude coordinates.
+    add_shift: int
+        Additional shift of data (is not applied to the longitude coordinate).
     model: aeolus.model.Model, optional
         Model class with a relevant longitude coordinate name.
 
@@ -959,16 +961,17 @@ def roll_cube_0_360(cube_in, model=um):
     lon = cube.coord(coord_name)
     if (lon.points < 0.0).any():
         add = 180
-        cube.data = da.roll(cube.core_data(), len(lon.points) // 2, axis=-1)
+        cube.data = da.roll(cube.core_data(), len(lon.points) // 2 + add_shift, axis=-1)
         if lon.has_bounds():
             bounds = lon.bounds + add
         else:
             bounds = None
-        cube.replace_coord(lon.copy(points=lon.points + add, bounds=bounds))
+        new_points = lon.points + add
+        cube.replace_coord(lon.copy(points=new_points, bounds=bounds))
     return cube
 
 
-def roll_cube_pm180(cube_in, model=um):
+def roll_cube_pm180(cube_in, add_shift=0, model=um):
     """
     Take a cube spanning 0...360 degrees in longitude and roll it to -180...180 degrees.
 
@@ -978,6 +981,8 @@ def roll_cube_pm180(cube_in, model=um):
     ----------
     cube: iris.cube.Cube
         Cube with longitude and latitude coordinates.
+    add_shift: int
+        Additional shift of data (is not applied to the longitude coordinate).
     model: aeolus.model.Model, optional
         Model class with a relevant longitude coordinate name.
 
@@ -996,7 +1001,7 @@ def roll_cube_pm180(cube_in, model=um):
         assert is_regular(xcoord), "Operation is only valid for a regularly spaced coordinate."
         if _is_longitude_global(xcoord.points):
             # Shift data symmetrically only when dealing with global cubes
-            cube.data = da.roll(cube.core_data(), len(xcoord.points) // 2, axis=-1)
+            cube.data = da.roll(cube.core_data(), len(xcoord.points) // 2 + add_shift, axis=-1)
 
         if xcoord.has_bounds():
             bounds = wrap_lons(xcoord.bounds, -180, 360)  # + subtract

src/aeolus/synthobs.py

@@ -6,6 +6,7 @@
 from iris.coords import AuxCoord, DimCoord
 from iris.cube import Cube
 from iris.exceptions import CoordinateNotFoundError as CoNotFound
+from iris.util import reverse
 
 import numpy as np
 
@@ -23,9 +24,9 @@
 )
 
 
-def calc_geom_mean_mirrored(cube_a, cube_b, model=um):
+def calc_geom_mean_mirrored(cube_a, cube_b, add_shift=0, model=um):
     """
-    Calculate geometric mean of two cubes with one of them rolled along the x-axis.
+    Calculate geometric mean of two cubes with one of them flipped along the x-axis.
 
     This function can be used to get an average transmission flux
     calculated separately from the day- and night-side perspective.
@@ -34,19 +35,21 @@ def calc_geom_mean_mirrored(cube_a, cube_b, model=um):
         Cube with an x-coordinate.
     cube_b: iris.cube.Cube
         Another cube with an x-coordinate to be flipped before being averaged with cube A.
+    add_shift: int
+        Additional shift for the data along the x-coordinate.
     model: aeolus.model.Model, optional
         Model class with relevant variable names.
     """
     # Get x-coordinate from the 1st cube.
     x_coord = cube_a.coord(model.x)
 
-    # Roll the 2nd cube by 180 degrees
+    # Roll and reverse the 2nd cube by 180 degrees
     # This is specific to the UM output
-    cube_b_rolled = roll_cube_pm180(cube_b, model=model)
-    cube_b_rolled.replace_coord(x_coord)
+    cube_b_flipped = reverse(roll_cube_pm180(cube_b, add_shift=add_shift, model=model), model.x)
+    cube_b_flipped.replace_coord(x_coord)
 
     # Calculate geometric mean of the two cubes
-    geom_mean = (cube_a * cube_b_rolled) ** 0.5
+    geom_mean = (cube_a * cube_b_flipped) ** 0.5
     return geom_mean
 
 
@@ -211,6 +214,7 @@ def calc_transmission_spectrum(
 def calc_transmission_spectrum_day_night_average(
     trans_flux_day,
     trans_flux_night,
+    add_shift=0,
     spectral_file=None,
     stellar_constant_at_1_au=None,
     stellar_radius=None,
@@ -239,6 +243,8 @@ def calc_transmission_spectrum_day_night_average(
         Transmission flux on spectral bands and optionally latitudes and longitudes.
         Night-side perspective.
         In the Met Office Unified Model this is STASH items 555, 556, 755, 756 in section 1.
+    add_shift: int, optional
+        Additional shift of data along the x-coordinate.
 
     For other parameters, see the docstring of `aeolus.synthobs.calc_transmission_spectrum()`
 
@@ -247,14 +253,16 @@ def calc_transmission_spectrum_day_night_average(
     aeolus.synthobs.calc_transmission_spectrum
     """
     # Average the day and night flux
-    trans_flux = calc_geom_mean_mirrored(trans_flux_day, trans_flux_night, model=model)
+    trans_flux = calc_geom_mean_mirrored(
+        trans_flux_day, trans_flux_night, add_shift=add_shift, model=model
+    )
     return calc_transmission_spectrum(
         trans_flux,
         spectral_file=spectral_file,
         stellar_constant_at_1_au=stellar_constant_at_1_au,
         stellar_radius=stellar_radius,
         planet_top_of_atmosphere=planet_top_of_atmosphere,
-        model=um,
+        model=model,
     )
 
 

tests/test_synthobs.py

@@ -20,12 +20,12 @@
 
 @pytest.fixture(scope="module")
 def example_trans_day():
-    return iris.load_cube(str(TST_DATA / "netcdf" / "planet_transmission_day.nc"))
+    return iris.load_cube(TST_DATA / "netcdf" / "planet_transmission_day.nc")
 
 
 @pytest.fixture(scope="module")
 def example_trans_night():
-    return iris.load_cube(str(TST_DATA / "netcdf" / "planet_transmission_night.nc"))
+    return iris.load_cube(TST_DATA / "netcdf" / "planet_transmission_night.nc")
 
 
 def test_read_spectral_bands():
@@ -75,11 +75,36 @@ def test_calc_stellar_flux():
 
 
 def test_calc_geom_mean_mirrored(example_trans_day, example_trans_night):
+    # Default
     actual = synthobs.calc_geom_mean_mirrored(example_trans_day, example_trans_night)
-    npt.assert_allclose(actual.data.max(), 1.7307187876145394e-08)
+    npt.assert_allclose(actual.data.max(), 1.5901232591606386e-08)
     npt.assert_allclose(
-        actual.data[123, 1, 71:74],
-        [2.40508136e-11, 2.40748625e-11, 2.40532320e-11],
+        actual.data[123, 45, 102:108],
+        [
+            0.00000000e00,
+            5.28650299e-40,
+            1.09253009e-10,
+            7.16136321e-10,
+            5.28943014e-10,
+            0.00000000e00,
+        ],
+    )
+    assert actual.units == example_trans_day.units
+    assert actual.shape == example_trans_day.shape
+    # With an additional shift along the x-coordinate
+    actual = synthobs.calc_geom_mean_mirrored(example_trans_day, example_trans_night, add_shift=-1)
+    npt.assert_allclose(actual.data.max(), 1.76899099563551e-08)
+    npt.assert_allclose(
+        actual.data[123, 45, 102:109],
+        [
+            0.00000000e00,
+            1.54026439e-67,
+            2.12892023e-11,
+            6.61482043e-10,
+            7.37627923e-10,
+            3.77198057e-10,
+            0.00000000e00,
+        ],
     )
     assert actual.units == example_trans_day.units
     assert actual.shape == example_trans_day.shape
@@ -129,21 +154,40 @@ def test_calc_transmission_spectrum(example_trans_day):
 
 
 def test_calc_transmission_spectrum_day_night_average(example_trans_day, example_trans_night):
+    # Default
     expected_rp_eff_over_rs = [
-        0.9999097552081786,
-        0.9999249847710132,
-        0.9999267911653156,
-        0.9999376935302737,
-        0.9999389664305547,
+        0.1289366838346542,
+        0.1289351109234838,
+        0.1289354291294527,
+        0.1289365978678745,
+        0.1289363841112755,
     ]
     actual_rp_eff_over_rs = synthobs.calc_transmission_spectrum_day_night_average(
         example_trans_day,
         example_trans_night,
-        TST_DATA / "spectral" / "sp_sw_500ir_bd_hatp11",
-        123,
-        123,
-        123,
+        spectral_file=TST_DATA / "spectral" / "sp_sw_500ir_bd_hatp11",
+        stellar_constant_at_1_au=1272.86475403,
+        stellar_radius=7.302834e08,
+        planet_top_of_atmosphere=94193200,
     )
     assert isinstance(actual_rp_eff_over_rs, iris.cube.Cube)
     assert actual_rp_eff_over_rs.shape[0] == 500
-    npt.assert_allclose(expected_rp_eff_over_rs, actual_rp_eff_over_rs.data[:5])
+    npt.assert_allclose(expected_rp_eff_over_rs, actual_rp_eff_over_rs.data[-5:])
+    # With an additional shift
+    expected_rp_eff_over_rs = [
+        0.1289324245796913,
+        0.1289306558725381,
+        0.1289309758113692,
+        0.1289323585757668,
+        0.1289320864595353,
+    ]
+    actual_rp_eff_over_rs = synthobs.calc_transmission_spectrum_day_night_average(
+        example_trans_day,
+        example_trans_night,
+        add_shift=-1,
+        spectral_file=TST_DATA / "spectral" / "sp_sw_500ir_bd_hatp11",
+        stellar_constant_at_1_au=1272.86475403,
+        stellar_radius=7.302834e08,
+        planet_top_of_atmosphere=94193200,
+    )
+    npt.assert_allclose(expected_rp_eff_over_rs, actual_rp_eff_over_rs.data[-5:])