Use earthkit-data instead of climetlab

docs/building/sources.rst

@@ -15,4 +15,4 @@ The following `sources` are currently available:
    sources/opendap
    sources/forcings
    sources/accumulations
-   sources/perturbations
+   sources/recentre

docs/building/sources/perturbations.rst → docs/building/sources/recentre.rst

@@ -1,29 +1,29 @@
-.. _perturbations:
+.. _recentre:
 
-###############
- perturbations
-###############
+##########
+ recentre
+##########
 
 Perturbations refers to the small variations centered around a nominal
 value of a parameter. When dealing with `ensemble forecasting`_, the
 perturbations are related to the difference between `ensemble members`
-and their given `center`.
+and their given `centre`.
 
-The `perturbations` function computes a set of new ensemble members
-centered on a different center from previous ensemble members using the
-following formula:
+The `recentre` function computes a set of new ensemble members centered
+on a different centre from previous ensemble members using the following
+formula:
 
 .. math::
 
-   members_{new} = center + ( members - \overline{members} )
+   members_{new} = centre + ( members - \overline{members} )
 
 Additionally, some variables must be non-negative to have a physical
 meaning (e.g. accumulated variables or `specific humidity`). To ensure
 this, positive clipping is performed using the alternative fomula :
 
 .. math::
 
-   members_{new} = max(0, center + ( members - \overline{members} ))
+   members_{new} = max(0, centre + ( members - \overline{members} ))
 
 The current implementation enforces that following variables are
 positive when using the `perturbations` function :
@@ -45,12 +45,12 @@ It uses the following arguments:
 members
    A :ref:`reference <yaml-reference>` to the ensemble members.
 
-center
+centre
    A :ref:`reference <yaml-reference>` to the new center requested.
 
 Examples
 
-.. literalinclude:: yaml/perturbations.yaml
+.. literalinclude:: yaml/recentre.yaml
    :language: yaml
 
 .. _convective precipitation: https://codes.ecmwf.int/grib/param-db/?id=143

src/anemoi/datasets/compute/perturbations.py → src/anemoi/datasets/compute/recentre.py

@@ -32,7 +32,7 @@
 SKIP = ("class", "stream", "type", "number", "expver", "_leg_number", "anoffset")
 
 
-def check_compatible(f1, f2, center_field_as_mars, ensemble_field_as_mars):
+def check_compatible(f1, f2, centre_field_as_mars, ensemble_field_as_mars):
     assert f1.mars_grid == f2.mars_grid, (f1.mars_grid, f2.mars_grid)
     assert f1.mars_area == f2.mars_area, (f1.mars_area, f2.mars_area)
     assert f1.shape == f2.shape, (f1.shape, f2.shape)
@@ -43,21 +43,22 @@ def check_compatible(f1, f2, center_field_as_mars, ensemble_field_as_mars):
         f2.metadata("valid_datetime"),
     )
 
-    for k in set(center_field_as_mars.keys()) | set(ensemble_field_as_mars.keys()):
+    for k in set(centre_field_as_mars.keys()) | set(ensemble_field_as_mars.keys()):
         if k in SKIP:
             continue
-        assert center_field_as_mars[k] == ensemble_field_as_mars[k], (
+        assert centre_field_as_mars[k] == ensemble_field_as_mars[k], (
             k,
-            center_field_as_mars[k],
+            centre_field_as_mars[k],
             ensemble_field_as_mars[k],
         )
 
 
-def perturbations(
+def recentre(
     *,
     members,
-    center,
+    centre,
     clip_variables=CLIP_VARIABLES,
+    alpha=1.0,
     output=None,
 ):
 
@@ -70,16 +71,16 @@ def perturbations(
 
     LOG.info("Ordering fields")
     members = members.order_by(*keys)
-    center = center.order_by(*keys)
+    centre = centre.order_by(*keys)
     LOG.info("Done")
 
-    if len(center) * n_numbers != len(members):
-        LOG.error("%s %s %s", len(center), n_numbers, len(members))
+    if len(centre) * n_numbers != len(members):
+        LOG.error("%s %s %s", len(centre), n_numbers, len(members))
         for f in members:
             LOG.error("Member: %r", f)
-        for f in center:
-            LOG.error("Center: %r", f)
-        raise ValueError(f"Inconsistent number of fields: {len(center)} * {n_numbers} != {len(members)}")
+        for f in centre:
+            LOG.error("centre: %r", f)
+        raise ValueError(f"Inconsistent number of fields: {len(centre)} * {n_numbers} != {len(members)}")
 
     if output is None:
         # prepare output tmp file so we can read it back
@@ -93,49 +94,44 @@ def perturbations(
 
     seen = set()
 
-    for i, center_field in enumerate(center):
-        param = center_field.metadata("param")
-        center_field_as_mars = center_field.metadata(namespace="mars")
+    for i, centre_field in enumerate(centre):
+        param = centre_field.metadata("param")
+        centre_field_as_mars = centre_field.metadata(namespace="mars")
 
-        # load the center field
-        center_np = center_field.to_numpy()
+        # load the centre field
+        centre_np = centre_field.to_numpy()
 
         # load the ensemble fields and compute the mean
-        members_np = np.zeros((n_numbers, *center_np.shape))
+        members_np = np.zeros((n_numbers, *centre_np.shape))
 
         for j in range(n_numbers):
             ensemble_field = members[i * n_numbers + j]
             ensemble_field_as_mars = ensemble_field.metadata(namespace="mars")
-            check_compatible(
-                center_field,
-                ensemble_field,
-                center_field_as_mars,
-                ensemble_field_as_mars,
-            )
+            check_compatible(centre_field, ensemble_field, centre_field_as_mars, ensemble_field_as_mars)
             members_np[j] = ensemble_field.to_numpy()
 
             ensemble_field_as_mars = tuple(sorted(ensemble_field_as_mars.items()))
             assert ensemble_field_as_mars not in seen, ensemble_field_as_mars
             seen.add(ensemble_field_as_mars)
 
-        # cmin=np.amin(center_np)
+        # cmin=np.amin(centre_np)
         # emin=np.amin(members_np)
 
         # if cmin < 0 and emin >= 0:
         #     LOG.warning(f"Negative values in {param} cmin={cmin} emin={emin}")
-        #     LOG.warning(f"Center: {center_field_as_mars}")
+        #     LOG.warning(f"centre: {centre_field_as_mars}")
 
         mean_np = members_np.mean(axis=0)
 
         for j in range(n_numbers):
             template = members[i * n_numbers + j]
             e = members_np[j]
             m = mean_np
-            c = center_np
+            c = centre_np
 
             assert e.shape == c.shape == m.shape, (e.shape, c.shape, m.shape)
 
-            x = c - m + e
+            x = c + (e - m) * alpha
 
             if param in clip_variables:
                 # LOG.warning(f"Clipping {param} to be positive")

src/anemoi/datasets/create/functions/filters/rotate_winds.py

@@ -70,11 +70,7 @@ def execute(
             lons,
             x.to_numpy(flatten=True),
             y.to_numpy(flatten=True),
-            (
-                source_projection
-                if source_projection is not None
-                else CRS.from_cf(x.grid_mapping)
-            ),
+            (source_projection if source_projection is not None else CRS.from_cf(x.grid_mapping)),
             target_projection,
         )
 

src/anemoi/datasets/create/functions/sources/perturbations.py → src/anemoi/datasets/create/functions/sources/recentre.py

@@ -8,7 +8,7 @@
 #
 from copy import deepcopy
 
-from anemoi.datasets.compute.perturbations import perturbations as compute_perturbations
+from anemoi.datasets.compute.recentre import recentre as compute_perturbations
 
 from .mars import mars
 
@@ -50,10 +50,10 @@ def load_if_needed(context, dates, dict_or_dataset):
     return dict_or_dataset
 
 
-def perturbations(context, dates, members, center, remapping={}, patches={}):
+def recentre(context, dates, members, centre, alpha=1.0, remapping={}, patches={}):
     members = load_if_needed(context, dates, members)
-    center = load_if_needed(context, dates, center)
-    return compute_perturbations(members, center)
+    centre = load_if_needed(context, dates, centre)
+    return compute_perturbations(members, centre, alpha=alpha)
 
 
-execute = perturbations
+execute = recentre

src/anemoi/datasets/grids.py

@@ -178,7 +178,7 @@ def cutout_mask(
         min_dlats = np.min(np.diff(glats))
         min_dlons = np.min(np.diff(glons))
 
-        # Use the center of the LAM grid as the reference point
+        # Use the centre of the LAM grid as the reference point
         centre = np.mean(lats), np.mean(lons)
         centre_xyz = np.array(latlon_to_xyz(*centre))
 
@@ -198,7 +198,7 @@ def cutout_mask(
         t = Triangle3D(lam_points[index[0]], lam_points[index[1]], lam_points[index[2]])
         # distance = np.min(distance)
         # The point is inside the triangle if the intersection with the ray
-        # from the point to the center of the Earth is not None
+        # from the point to the centre of the Earth is not None
         # (the direction of the ray is not important)
 
         intersect = t.intersect(zero, global_point)

tests/create/perturbations.yaml → tests/create/recentre.yaml

@@ -36,7 +36,7 @@ common:
     type: an
     number: [1, 2, 4]
     # number: [1, 2, 3, 4, 5, 6, 7, 8, 9]
-  center: &center
+  centre: &centre
     stream: oper
     type: an
 
@@ -53,27 +53,27 @@ data_sources:
       - accumulations:
           <<: *ensembles
           <<: *acc
-  center:
+  centre:
     join:
       - mars:
-          <<: *center
+          <<: *centre
           <<: *sfc
       - mars:
-          <<: *center
+          <<: *centre
           <<: *pl
       - accumulations:
-          <<: *center
+          <<: *centre
           <<: *acc
 
 input:
   join:
-    - perturbations:
+    - recentre:
         # the ensemble data which has one additional dimension
         members: ${data_sources.ensembles}
-        # the new center of the data
-        center: ${data_sources.center}
+        # the new centre of the data
+        centre: ${data_sources.centre}
     - forcings:
-        template: ${input.join.0.perturbations}
+        template: ${input.join.0.recentre}
         param:
           - cos_latitude
           - cos_longitude

tests/create/test_create.py

@@ -29,8 +29,7 @@
 HERE = os.path.dirname(__file__)
 # find_yamls
 NAMES = sorted([os.path.basename(path).split(".")[0] for path in glob.glob(os.path.join(HERE, "*.yaml"))])
-
-SKIP = ["perturbations"]
+SKIP = ["recentre"]
 NAMES = [name for name in NAMES if name not in SKIP]
 assert NAMES, "No yaml files found in " + HERE