Fix Earth area weight computation for reduced Gaussian grids

CHANGES.txt

@@ -1,6 +1,10 @@
 Next Release
 =============
 
+Bugfixes
+---------
+- Fix Earth area computation for reduced Gaussian grid (PR #127)
+
 Internal changes
 -----------------
 - Move `get_template` into `helpers.files` (PR #125)

helpers/cubes.py

@@ -5,6 +5,7 @@
 import iris
 import iris.analysis.cartography
 import iris.cube
+import iris.warnings
 import numpy as np
 from iris.util import equalise_attributes
 from scriptengine.exceptions import (
@@ -185,7 +186,7 @@ def compute_reduced_grid_weights(cube):
     unique_lats, gridpoints_per_lat = unique_lats[0:-1], gridpoints_per_lat[0:-1]
     areas = []
     last_angle = 0
-    earth_radius = 6371
+    earth_radius = 6.371229e6  # m
     for latitude, amount in zip(unique_lats, gridpoints_per_lat):
         delta = latitude - last_angle
         current_angle = last_angle + 2 * delta
@@ -203,7 +204,14 @@ def compute_regular_grid_weights(cube):
     """compute area weights for a regular lat/lon grid"""
     cube.coord("latitude").guess_bounds()
     cube.coord("longitude").guess_bounds()
-    return iris.analysis.cartography.area_weights(cube)
+    with warnings.catch_warnings():
+        # Suppress default radius warning
+        warnings.filterwarnings(
+            action="ignore",
+            message="Using DEFAULT_SPHERICAL_EARTH_RADIUS",
+            category=iris.warnings.IrisDefaultingWarning,
+        )
+        return iris.analysis.cartography.area_weights(cube)
 
 
 def align_time_coords(new_cube, old_cube):

tests/test_oifs_timeseries.py

@@ -29,7 +29,7 @@ def test_oifs_global_mean_year_mean_timeseries_working(tmp_path):
     )
 
 
-def test_oifs_timeseries_compare_grids(tmp_path):
+def test_oifs_timeseries_compare_grids_mean(tmp_path):
     init = {
         "src": ["./tests/testdata/regular_grid_tas.nc"],
         "dst": str(tmp_path / "test_reg.nc"),
@@ -50,6 +50,27 @@ def test_oifs_timeseries_compare_grids(tmp_path):
     assert abs(cube_red.data - cube_reg.data) < 1e-3
 
 
+def test_oifs_timeseries_compare_grids_sum(tmp_path):
+    init = {
+        "src": ["./tests/testdata/regular_grid_tas.nc"],
+        "dst": str(tmp_path / "test_reg.nc"),
+        "varname": "tas",
+    }
+    atmo_ts = OifsGlobalSumYearMeanTimeseries(init)
+    atmo_ts.run(init)
+    cube_reg = iris.load_cube(init["dst"])
+
+    init = {
+        "src": ["./tests/testdata/reduced_grid_tas.nc"],
+        "dst": str(tmp_path / "test_red.nc"),
+        "varname": "tas",
+    }
+    atmo_ts = OifsGlobalSumYearMeanTimeseries(init)
+    atmo_ts.run(init)
+    cube_red = iris.load_cube(init["dst"])
+    assert abs(cube_red.data - cube_reg.data) < 1e-2 * abs(cube_red.data)
+
+
 def test_oifs_global_mean_year_mean_timeseries_wrong_varname(tmp_path):
     init = {
         "src": ["./tests/testdata/TES1_atm_1m_1990_2t.nc"],