fix cutout

anemoi/datasets/grids.py

@@ -121,7 +121,7 @@ def cutout_mask(
     global_lats,
     global_lons,
     cropping_distance=2.0,
-    min_distance_km=0.0,
+    min_distance_km=None,
     plot=None,
 ):
     """
@@ -131,8 +131,6 @@ def cutout_mask(
 
     # TODO: transform min_distance from lat/lon to xyz
 
-    min_distance = min_distance_km / 6371.0
-
     assert global_lats.ndim == 1
     assert global_lons.ndim == 1
     assert lats.ndim == 1
@@ -172,6 +170,30 @@ def cutout_mask(
     kdtree = KDTree(lam_points)
     distances, indices = kdtree.query(global_points, k=3)
 
+    if min_distance_km is not None:
+        min_distance = min_distance_km / 6371.0
+    else:
+        # Estimnation of the minimum distance between two grib points
+
+        glats = sorted(set(global_lats_masked))
+        glons = sorted(set(global_lons_masked))
+        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
+        centre = np.mean(lats), np.mean(lons)
+        centre_xyz = np.array(latlon_to_xyz(*centre))
+
+        pt1 = np.array(latlon_to_xyz(centre[0] + min_dlats, centre[1]))
+        pt2 = np.array(latlon_to_xyz(centre[0], centre[1] + min_dlons))
+        min_distance = (
+            min(
+                np.linalg.norm(pt1 - centre_xyz),
+                np.linalg.norm(pt2 - centre_xyz),
+            )
+            / 2.0
+        )
+
     zero = np.array([0.0, 0.0, 0.0])
     ok = []
     for i, (global_point, distance, index) in enumerate(zip(global_points, distances, indices)):
@@ -181,10 +203,10 @@ def cutout_mask(
         # from the point to the center of the Earth is not None
         # (the direction of the ray is not important)
 
-        intersect = t.intersect(zero, global_point) or t.intersect(global_point, zero)
+        intersect = t.intersect(zero, global_point)
         close = np.min(distance) <= min_distance
 
-        ok.append(intersect and not close)
+        ok.append(intersect or close)
 
     j = 0
     ok = np.array(ok)

tools/Makefile

@@ -0,0 +1,9 @@
+TARGETS=dataset1.zarr dataset2.zarr
+
+all: $(TARGETS)
+
+%.zarr: %.yaml
+	anemoi-datasets create $< $@ --overwrite
+
+
+.SUFFIXES: .zarr .yaml

tools/dataset1.yaml

@@ -0,0 +1,40 @@
+common:
+  mars_request: &mars_request
+    expver: "0001"
+    grid: 1/1
+
+dates:
+  start: 2024-01-01 00:00:00
+  end: 2024-01-01 18:00:00
+  frequency: 6h
+
+input:
+  join:
+  - mars:
+      <<: *mars_request
+      param: [2t, 10u, 10v, lsm]
+      levtype: sfc
+      stream: oper
+      type: an
+  - mars:
+      <<: *mars_request
+      param: [q, t, z]
+      levtype: pl
+      level: [50, 100]
+      stream: oper
+      type: an
+  - accumulations:
+      <<: *mars_request
+      levtype: sfc
+      param: [cp, tp]
+  - constants:
+      template: ${input.join.0.mars}
+      param:
+      - cos_latitude
+      - sin_latitude
+
+output:
+  order_by: [valid_datetime, param_level, number]
+  remapping:
+    param_level: "{param}_{levelist}"
+  statistics: param_level

tools/dataset2.yaml

@@ -0,0 +1,42 @@
+common:
+  mars_request: &mars_request
+    expver: "0001"
+    grid: 0.5/0.5
+    area: [28, 0, -14, 40]
+    rotation: [-20, -40]
+
+dates:
+  start: 2024-01-01 00:00:00
+  end: 2024-01-01 18:00:00
+  frequency: 6h
+
+input:
+  join:
+  - mars:
+      <<: *mars_request
+      param: [2t, 10u, 10v, lsm]
+      levtype: sfc
+      stream: oper
+      type: an
+  - mars:
+      <<: *mars_request
+      param: [q, t, z]
+      levtype: pl
+      level: [50, 100]
+      stream: oper
+      type: an
+  - accumulations:
+      <<: *mars_request
+      levtype: sfc
+      param: [cp, tp]
+  - constants:
+      template: ${input.join.0.mars}
+      param:
+      - cos_latitude
+      - sin_latitude
+
+output:
+  order_by: [valid_datetime, param_level, number]
+  remapping:
+    param_level: "{param}_{levelist}"
+  statistics: param_level