Fix/grid consistency check

src/earthkit/data/core/field.py

@@ -525,6 +525,9 @@ def from_components(
             )
         elif geography is not None:
             components[_GEOGRAPHY] = _COMPONENT_MAKER.default_cls(_GEOGRAPHY).from_any(geography)
+        elif geography is None and shape_hint is not None:
+            # create EmptyGeography with shape_hint
+            components[_GEOGRAPHY] = _COMPONENT_MAKER.default_cls(_GEOGRAPHY).from_dict({}, shape_hint=shape_hint)
         return cls(**components)
 
     @property

src/earthkit/data/indexing/indexed.py

@@ -431,7 +431,16 @@ def geography(self):
         elif len(self) == 0:
             raise ValueError("Cannot determine geography of an empty FieldList")
         else:
-            raise ValueError("Fields do not have the same grid geometry")
+            field_iter = iter(self)
+            field0_grid_id = next(field_iter).geography.unique_grid_id()
+            for i, f in enumerate(field_iter, start=1):
+                other_field_grid_id = f.geography.unique_grid_id()
+                if other_field_grid_id != field0_grid_id:
+                    break
+            raise ValueError(
+                f"Fields do not have the same grid geometry: field 0 has grid id {field0_grid_id}, "
+                f"but field {i} has grid id {other_field_grid_id}; field list length is {len(self)}"
+            )
 
     @thread_safe_cached_property
     def _has_shared_geography(self):

src/earthkit/data/xr_engine/builder.py

@@ -643,11 +643,12 @@ def parse(self, ds, profile=None, full=False):
             ds,
             keys=profile.index_keys,
             remapping=remapping,
-            # component=profile.add_earthkit_attrs
-            # PW: check the role of the `component` kwarg
         )
         # LOG.debug(f"{ds.db=}")
 
+        # check grid consistency across fields
+        ds_xr.geography.unique_grid_id()
+
         # LOG.debug(f"before update: {profile.dim_keys=}")
         profile.update(ds_xr)
         # LOG.debug(f"after update: {profile.dim_keys=}")
@@ -691,10 +692,10 @@ def __init__(self, *args, from_xr=False, **kwargs):
         super().__init__(*args, **kwargs)
 
         if self.split_dims:
-            raise ValueError("SingleDatasetMaker does not support splitting")
+            raise ValueError("SingleDatasetBuilder does not support splitting")
 
         if from_xr and self.direct_backend:
-            raise ValueError("SingleDatasetMaker does not support direct_backend=True when invoked from xarray")
+            raise ValueError("SingleDatasetBuilder does not support direct_backend=True when invoked from xarray")
 
     def build(self):
         ds_sorted, _ = self.parse(self.ds, self.profile)
@@ -726,7 +727,7 @@ def __init__(self, *args, backend_kwargs=None, other_kwargs=None):
             self.xr_open_dataset_kwargs["backend_kwargs"] = backend_kwargs
 
         if not self.split_dims:
-            raise ValueError("SplitDatasetMaker requires split_dims")
+            raise ValueError("SplitDatasetBuilder requires split_dims")
 
     def prepare(self, keys):
         from .fieldlist import XArrayInputFieldList
@@ -741,13 +742,23 @@ def prepare(self, keys):
 
         return ds_xr, vals
 
-    def build(self):
-        from .splitter import Splitter
+    def split(self):
+        ds_xr, dims = self.prepare(self.split_dims)
 
-        splitter = Splitter.make(self.split_dims)
+        for x in itertools.product(*dims.values()):
+            y = dict(zip(dims.keys(), x))
+            ds_sel = ds_xr.sel(y)
+            if len(ds_sel) == 0:
+                continue
+            ds_sort, profile = self.parse(ds_sel, None)
+            if len(ds_sort) == 0:
+                raise ValueError(f"No field found for selection={y}")
+            yield ds_sort, profile, y
+
+    def build(self):
         datasets = []
         split_coords_list = []
-        for ds, profile, split_coords in splitter.split(self):
+        for ds, profile, split_coords in self.split():
             dims = profile.dims.to_list()
             split_coords_list.append(dict(split_coords))
             LOG.debug(f"splitting {dims=} type of s_ds={type(ds)} {split_coords=}")

tests/netcdf/test_netcdf_output.py

@@ -80,7 +80,7 @@ def test_netcdf_fieldlist_multi_subset_save_bad():
     assert len(ds) == 6
 
     with temp_file() as tmp:
-        with pytest.raises((AttributeError, IndexError)):
+        with pytest.raises((AttributeError, IndexError, ValueError)):
             ds.to_target("file", tmp)
 
 

tests/xr_engine/test_xr_engine_grid.py

@@ -114,3 +114,13 @@ def test_xr_engine_gridspec(allow_holes, lazy_load):
     gs = r.earthkit.grid_spec
     #  assert gs["type"] == "reduced_gg"
     assert gs["grid"] == "O32"
+
+
+@pytest.mark.cache
+@pytest.mark.parametrize("allow_holes", [False, True])
+@pytest.mark.parametrize("lazy_load", [True, False])
+def test_xr_engine_different_grids(allow_holes, lazy_load):
+    url = earthkit_remote_test_data_file("xr_engine", "grid", "different_grid_resolutions.grib")
+    ds = from_source("url", url).to_fieldlist()
+    with pytest.raises(ValueError, match="Fields do not have the same grid geometry"):
+        ds.to_xarray(allow_holes=allow_holes, lazy_load=lazy_load)

tests/xr_engine/test_xr_engine_lod.py

@@ -52,6 +52,28 @@ def xr_lod_nongeo():
     return ds
 
 
+@pytest.fixture
+def xr_lod_nongeo_many_shapes():
+    prototype1 = {
+        "values": [1, 2, 3, 4, 5, 6],
+        "time": {"valid_datetime": "2018-08-01T09:00:00Z"},
+    }
+
+    prototype2 = {
+        "values": [1, 2, 3, 4, 5, 6, 7],
+        "time": {"valid_datetime": "2018-08-01T09:00:00Z"},
+    }
+
+    d = [
+        {"parameter": {"variable": "t"}, "vertical": {"level": 500}, **prototype1},
+        {"parameter": {"variable": "t"}, "vertical": {"level": 850}, **prototype1},
+        {"parameter": {"variable": "u"}, "vertical": {"level": 500}, **prototype2},
+        {"parameter": {"variable": "u"}, "vertical": {"level": 850}, **prototype2},
+    ]
+    ds = from_source("list-of-dicts", d).to_fieldlist()
+    return ds
+
+
 @pytest.fixture
 def xr_lod_forecast():
     prototype = {
@@ -187,6 +209,16 @@ def test_xr_engine_lod_nongeo(allow_holes, lazy_load, xr_lod_nongeo):
     assert np.allclose(ds["u"].values, ref)
 
 
+@pytest.mark.parametrize("allow_holes", [False, True])
+@pytest.mark.parametrize("lazy_load", [True, False])
+def test_xr_engine_lod_nongeo_many_shapes(allow_holes, lazy_load, xr_lod_nongeo_many_shapes):
+    ds_in = xr_lod_nongeo_many_shapes
+    with pytest.raises(ValueError, match="Fields do not have the same grid geometry"):
+        ds_in.to_xarray(
+            profile="earthkit", time_dims=["date", "time", "step"], allow_holes=allow_holes, lazy_load=lazy_load
+        )
+
+
 @pytest.mark.parametrize("allow_holes", [False, True])
 @pytest.mark.parametrize("lazy_load", [True, False])
 def test_xr_engine_lod_forecast(allow_holes, lazy_load, xr_lod_forecast):