Repository: ecmwf/earthkit-data
Pull request #662 (2025-03-24 10:51:19)

Feature/merge 0.13.5 into develop


Patch:
diff --git a/docs/release_notes/version_0.13_updates.rst b/docs/release_notes/version_0.13_updates.rst
index 214b5e1..a874180 100644
--- a/docs/release_notes/version_0.13_updates.rst
+++ b/docs/release_notes/version_0.13_updates.rst
@@ -2,6 +2,16 @@ Version 0.13 Updates
 /////////////////////////
 
 
+
+Version 0.13.5
+===============
+
+Fixes
+++++++++++
+
+- Fixed issue when NetCDF data with 1D latitudes/longitudes appearing as variables was not parsed into a fieldlist when read with :func:`from_source` (:pr:`658`)
+
+
 Version 0.13.4
 ===============
 
diff --git a/src/earthkit/data/readers/netcdf/dataset.py b/src/earthkit/data/readers/netcdf/dataset.py
index bc2dc6b..75f83f7 100644
--- a/src/earthkit/data/readers/netcdf/dataset.py
+++ b/src/earthkit/data/readers/netcdf/dataset.py
@@ -41,8 +41,8 @@ class DataSet:
     def bbox(self, variable):
         data_array = self[variable]
 
-        keys, coords = self._get_xy_coords(data_array)
-        key = ("bbox", tuple(keys), tuple(coords))
+        keys, dims = self._get_xy_dims(data_array)
+        key = ("bbox", tuple(keys), tuple(dims))
         if key in self._cache:
             return self._cache[key]
 
@@ -59,7 +59,10 @@ class DataSet:
         self._cache[key] = (north, west, south, east)
         return self._cache[key]
 
-    def _get_xy_coords(self, data_array):
+    def _get_xy_dims(self, data_array):
+        # TODO: refactor it
+
+        # the last 2 dimensions are y,x
         if (
             len(data_array.dims) >= 2
             and data_array.dims[-1] in GEOGRAPHIC_COORDS["x"]
@@ -67,40 +70,68 @@ class DataSet:
         ):
             return ("y", "x"), (data_array.dims[-2], data_array.dims[-1])
 
+        # the last 2 dimensions are x, y
+        if (
+            len(data_array.dims) >= 2
+            and data_array.dims[-1] in GEOGRAPHIC_COORDS["y"]
+            and data_array.dims[-2] in GEOGRAPHIC_COORDS["x"]
+        ):
+            return ("x", "y"), (data_array.dims[-2], data_array.dims[-1])
+
+        # try to find the x and y dimensions
         keys = []
-        coords = []
+        dims = []
         axes = ("x", "y")
         for dim in data_array.dims:
             for ax in axes:
                 candidates = GEOGRAPHIC_COORDS.get(ax, [])
                 if dim in candidates:
                     keys.append(ax)
-                    coords.append(dim)
+                    dims.append(dim)
                 else:
                     ax = data_array.coords[dim].attrs.get("axis", "").lower()
                     if ax in axes:
                         keys.append(ax)
-                        coords.append(dim)
+                        dims.append(dim)
             if len(keys) == 2:
-                return tuple(keys), tuple(coords)
+                return tuple(keys), tuple(dims)
+
+        # 1D geographic coordinates using the dim='values' convention
+        if not keys or not dims:
+            dim = "values"
+            if data_array.dims and data_array.dims[-1] == dim:
+                for x, y in zip(GEOGRAPHIC_COORDS["x"], GEOGRAPHIC_COORDS["y"]):
+                    if (x in data_array.coords or x in self._ds.variables) and (
+                        y in data_array.coords or y in self._ds.variables
+                    ):
+                        if (
+                            len(self._ds[x].dims) == 1
+                            and len(self._ds[y].dims) == 1
+                            and self._ds[x].dims[-1] == dim
+                            and self._ds[y].dims[-1] == dim
+                        ):
+                            return tuple(["x"]), tuple([dim])
 
         for ax in axes:
             if ax not in keys:
-                raise ValueError(f"No coordinate found with axis '{ax}'")
+                raise ValueError(f"No dimension found with axis '{ax}'")
 
-        return keys, coords
+        return keys, dims
 
     def _get_xy(self, data_array, flatten=False, dtype=None):
-        keys, coords = self._get_xy_coords(data_array)
-        key = ("grid_points", tuple(keys), tuple(coords))
+        # TODO: refactor it
+
+        keys, dims = self._get_xy_dims(data_array)
+        key = ("grid_points", tuple(keys), tuple(dims))
 
         if key in self._cache:
             points = self._cache[key]
         else:
             points = dict()
-            v0, v1 = data_array.coords[coords[0]], data_array.coords[coords[1]]
-            points[keys[1]], points[keys[0]] = np.meshgrid(v1, v0)
-            self._cache[key] = points
+            if all(d in data_array.coords for d in dims):
+                v0, v1 = data_array.coords[dims[0]], data_array.coords[dims[1]]
+                points[keys[1]], points[keys[0]] = np.meshgrid(v1, v0)
+                self._cache[key] = points
 
         if flatten:
             points["x"] = points["x"].reshape(-1)
@@ -112,7 +143,8 @@ class DataSet:
             return points["x"], points["y"]
 
     def _get_latlon(self, data_array, flatten=False, dtype=None):
-        keys, coords = self._get_xy_coords(data_array)
+        # TODO: refactor it
+        keys, dims = self._get_xy_dims(data_array)
 
         points = dict()
 
@@ -126,22 +158,27 @@ class DataSet:
         latitude = _get_ll(lat_keys)
         longitude = _get_ll(lon_keys)
 
+        # lat-lon variable/coordinate available
         if latitude is not None and longitude is not None:
-            if latitude.dims == coords and longitude.dims == coords:
+            if latitude.dims == dims and longitude.dims == dims:
                 latitude = latitude.data
                 longitude = longitude.data
                 points["y"] = latitude
                 points["x"] = longitude
 
+        # lat-lon meshgrid
         if not points:
-            key = ("grid_points", tuple(keys), tuple(coords))
+            key = ("grid_points", tuple(keys), tuple(dims))
 
             if key in self._cache:
                 points = self._cache[key]
             else:
-                v0, v1 = data_array.coords[coords[0]], data_array.coords[coords[1]]
-                points[keys[1]], points[keys[0]] = np.meshgrid(v1, v0)
-                self._cache[key] = points
+                if all(d in data_array.coords for d in dims):
+                    v0, v1 = data_array.coords[dims[0]], data_array.coords[dims[1]]
+                    points[keys[1]], points[keys[0]] = np.meshgrid(v1, v0)
+                    self._cache[key] = points
+                else:
+                    raise ValueError(f"Could not find lat-lon coordinates for {data_array.name}")
 
         if flatten:
             points["x"] = points["x"].reshape(-1)
diff --git a/src/earthkit/data/readers/netcdf/field.py b/src/earthkit/data/readers/netcdf/field.py
index 5114890..589f247 100644
--- a/src/earthkit/data/readers/netcdf/field.py
+++ b/src/earthkit/data/readers/netcdf/field.py
@@ -47,8 +47,8 @@ class XArrayFieldGeography(Geography):
         return self.ds._get_xy(self.data_array, flatten=True, dtype=dtype)[1]
 
     def shape(self):
-        _, coords = self.ds._get_xy_coords(self.data_array)
-        return tuple([self.data_array.coords[v].size for v in coords])
+        _, dims = self.ds._get_xy_dims(self.data_array)
+        return tuple([self.ds.sizes[v] for v in dims])
 
     def _unique_grid_id(self):
         return self.shape()
diff --git a/src/earthkit/data/readers/netcdf/fieldlist.py b/src/earthkit/data/readers/netcdf/fieldlist.py
index ad33113..e291f5d 100644
--- a/src/earthkit/data/readers/netcdf/fieldlist.py
+++ b/src/earthkit/data/readers/netcdf/fieldlist.py
@@ -28,6 +28,19 @@ from .field import XArrayField
 LOG = logging.getLogger(__name__)
 
 
+def get_geo_vars(ds):
+    "" "Return the names and dims if the geographic coordinates are variables in the dataset" ""
+    for x, y in zip(GEOGRAPHIC_COORDS["x"], GEOGRAPHIC_COORDS["y"]):
+        if x in ds.variables and y in ds.variables:
+            dims_x = set(ds[x].dims)
+            dims_y = set(ds[y].dims)
+            dims = set()
+            if dims_x == dims_y:
+                dims = dims_x
+            return (x, y), dims
+    return None, None
+
+
 def get_fields_from_ds(
     ds,
     field_type=None,
@@ -43,6 +56,11 @@ def get_fields_from_ds(
         if isinstance(attr_val, str):
             skip.update(attr_val.split(" "))
 
+    skip_vars = set()
+    geo_vars, geo_vars_dims = get_geo_vars(ds)
+    if geo_vars:
+        skip_vars.update(geo_vars)
+
     for name in ds.data_vars:
         v = ds[name]
         _skip_attr(v, "coordinates")
@@ -50,9 +68,12 @@ def get_fields_from_ds(
         _skip_attr(v, "grid_mapping")
 
     for name in ds.data_vars:
+        if name in skip_vars:
+            continue
+
         # Select only geographical variables
-        has_lat = False
-        has_lon = False
+        has_lat = None
+        has_lon = None
 
         if name in skip:
             continue
@@ -66,6 +87,7 @@ def get_fields_from_ds(
         info = [value for value in v.coords if value not in v.dims]
         non_dim_coords = {}
         for coord in v.coords:
+            # print(f"{coord=}")
             if coord not in v.dims:
                 non_dim_coords[coord] = ds[coord].values
                 continue
@@ -81,7 +103,7 @@ def get_fields_from_ds(
 
             # LOG.debug(f"{standard_name=} {long_name=} {axis=} {coord_name}")
             use = False
-
+            # print(f"{standard_name=} {long_name=} {axis=} {coord_name}")
             if (
                 standard_name.lower() in GEOGRAPHIC_COORDS["x"]
                 or (long_name == "longitude")
@@ -134,6 +156,12 @@ def get_fields_from_ds(
             if not use:
                 coordinates.append(OtherCoordinate(c, coord in info))
 
+        if not (has_lat and has_lon):
+            # handle the case where the lat/lon are variables
+            if geo_vars and geo_vars_dims and all([d in v.dims for d in geo_vars_dims]):
+                has_lat = True
+                has_lon = True
+
         if not (has_lat and has_lon):
             # self.log.info("NetCDFReader: skip %s (Not a 2 field)", name)
             continue


Test patch:
diff --git a/tests/netcdf/test_netcdf_geography.py b/tests/netcdf/test_netcdf_geography.py
index fa9aae4..e342ec8 100644
--- a/tests/netcdf/test_netcdf_geography.py
+++ b/tests/netcdf/test_netcdf_geography.py
@@ -12,6 +12,7 @@
 import numpy as np
 import pytest
 
+from earthkit.data import from_object
 from earthkit.data import from_source
 from earthkit.data.testing import earthkit_examples_file
 from earthkit.data.testing import earthkit_remote_test_data_file
@@ -251,12 +252,217 @@ def test_netcdf_to_latlon_laea():
             assert np.isclose(v["lat"][x], ref[i]), f"{i=}, {x=}"
 
 
-def test_netcdf_forecast_reference_time():
-    ds = from_source("url", earthkit_remote_test_data_file("test-data", "fa_ta850.nc"))
+@pytest.mark.parametrize("lat_name,lon_name", [("lat", "lon"), ("latitude", "longitude")])
+def test_netcdf_geography_2d_1(lat_name, lon_name):
+    # Dimensions:  (level: 2, lat: 3, lon: 3)
+    # Coordinates:
+    #   * level    (level) int64 16B 700 500
+    #   * lat      (lat) int64 24B 50 40 30
+    #   * lon      (lon) int64 24B 0 10 20
+    # Data variables:
+    #     a        (level, lat, lon) int64 144B 11 12 13 21 22 23 ... 25 26 34 35 36
 
-    assert len(ds) == 37
-    assert ds[0].metadata("valid_datetime") == "2020-01-23T00:00:00"
-    assert ds[5].metadata("valid_datetime") == "2020-01-23T05:00:00"
+    import xarray as xr
+
+    dims = {"level": 2, lat_name: 3, lon_name: 3}
+    coords = {
+        "level": np.array([700, 500]),
+        lat_name: np.array([50, 40, 30]),
+        lon_name: np.array([0, 10, 20]),
+    }
+
+    lats = [[50, 50, 50], [40, 40, 40], [30, 30, 30]]
+    lons = [[0, 10, 20], [0, 10, 20], [0, 10, 20]]
+
+    data = np.array(
+        [
+            [[11, 12, 13], [21, 22, 23], [31, 32, 33]],
+            [[14, 15, 16], [24, 25, 26], [34, 35, 36]],
+        ]
+    )
+
+    a = xr.Variable(dims, data)
+    v = {"a": a}
+    ds_in = xr.Dataset(v, coords=coords)
+
+    ds = from_object(ds_in)
+    assert len(ds) == 2
+    assert np.allclose(ds.metadata("level"), coords["level"])
+
+    for ll in [ds[0].to_latlon(), ds.to_latlon()]:
+        assert ll["lat"].shape == (3, 3)
+        assert ll["lon"].shape == (3, 3)
+        assert np.allclose(ll["lat"], lats)
+        assert np.allclose(ll["lon"], lons)
+
+
+@pytest.mark.parametrize("lat_name,lon_name", [("lat", "lon"), ("latitude", "longitude")])
+def test_netcdf_geography_2d_2(lat_name, lon_name):
+    # Dimensions:  (level: 2, y: 3, x: 2)
+    # Coordinates:
+    #   * level    (level) int64 16B 700 500
+    #     lat      (y, x) int64 48B 50 50 40 40 30 30
+    #     lon      (y, x) int64 48B 0 10 0 10 0 10
+    # Dimensions without coordinates: y, x
+    # Data variables:
+    #     a        (level, y, x) int64 96B 11 12 21 22 31 32 14 15 24 25 34 35
+
+    import xarray as xr
+
+    dims = {"level": 2, "y": 3, "x": 2}
+    coords = {
+        "level": np.array([700, 500]),
+        lat_name: (["y", "x"], np.array([[50, 50], [40, 40], [30, 30]])),
+        lon_name: (["y", "x"], np.array([[0, 10], [0, 10], [0, 10]])),
+    }
+
+    data = np.array(
+        [
+            [[11, 12], [21, 22], [31, 32]],
+            [[14, 15], [24, 25], [34, 35]],
+        ]
+    )
+
+    a = xr.Variable(dims, data)
+    v = {"a": a}
+    ds_in = xr.Dataset(v, coords=coords)
+
+    ds = from_object(ds_in)
+    assert len(ds) == 2
+    assert np.allclose(ds.metadata("level"), coords["level"])
+
+    for ll in [ds[0].to_latlon(), ds.to_latlon()]:
+        assert ll["lat"].shape == (3, 2)
+        assert ll["lon"].shape == (3, 2)
+        assert np.allclose(ll["lat"], coords[lat_name][1])
+        assert np.allclose(ll["lon"], coords[lon_name][1])
+
+
+@pytest.mark.parametrize("lat_name,lon_name", [("lat", "lon"), ("latitude", "longitude")])
+def test_netcdf_geography_2d_3(lat_name, lon_name):
+    # Dimensions:  (level: 2, y: 3, x: 2)
+    # Coordinates:
+    #   * level    (level) int64 16B 700 500
+    # Dimensions without coordinates: y, x
+    # Data variables:
+    #     a        (level, y, x) int64 96B 11 12 21 22 31 32 14 15 24 25 34 35
+    #     lat      (y, x) int64 48B 50 50 40 40 30 30
+    #     lon      (y, x) int64 48B 0 10 0 10 0 10
+
+    import xarray as xr
+
+    dims = {"level": 2, "y": 3, "x": 2}
+    coords = {
+        "level": np.array([700, 500]),
+    }
+
+    data = np.array(
+        [
+            [[11, 12], [21, 22], [31, 32]],
+            [[14, 15], [24, 25], [34, 35]],
+        ]
+    )
+
+    a = xr.Variable(dims, data)
+    lat = xr.Variable({"y": 3, "x": 3}, np.array([[50, 50], [40, 40], [30, 30]]))
+    lon = xr.Variable({"y": 3, "x": 3}, np.array([[0, 10], [0, 10], [0, 10]]))
+    v = {"a": a, lat_name: lat, lon_name: lon}
+    ds_in = xr.Dataset(v, coords=coords)
+
+    ds = from_object(ds_in)
+    assert len(ds) == 2
+    assert np.allclose(ds.metadata("level"), coords["level"])
+
+    for ll in [ds[0].to_latlon(), ds.to_latlon()]:
+        assert ll["lat"].shape == (3, 2)
+        assert ll["lon"].shape == (3, 2)
+        assert np.allclose(ll["lat"], lat.data)
+        assert np.allclose(ll["lon"], lon.data)
+
+
+@pytest.mark.parametrize("lat_name,lon_name", [("lat", "lon"), ("latitude", "longitude")])
+def test_netcdf_geography_1d_1(lat_name, lon_name):
+    # Dimensions:  (level: 2, values: 9)
+    # Coordinates:
+    #   * level    (level) int64 16B 700 500
+    #     lat      (values) int64 72B 50 50 50 40 40 40 30 30 30
+    #     lon      (values) int64 72B 0 10 20 0 10 20 0 10 20
+    # Dimensions without coordinates: values
+    # Data variables:
+    #     a        (level, values) int64 144B 11 12 13 21 22 23 ... 24 25 26 34 35 36
+
+    import xarray as xr
+
+    dims = {"level": 2, "values": 9}
+    coords = {
+        "level": np.array([700, 500]),
+        lat_name: ("values", np.array([50, 50, 50, 40, 40, 40, 30, 30, 30])),
+        lon_name: ("values", np.array([0, 10, 20, 0, 10, 20, 0, 10, 20])),
+    }
+
+    data = np.array(
+        [
+            [11, 12, 13, 21, 22, 23, 31, 32, 33],
+            [14, 15, 16, 24, 25, 26, 34, 35, 36],
+        ]
+    )
+
+    a = xr.Variable(dims, data)
+    v = {"a": a}
+    ds_in = xr.Dataset(v, coords=coords)
+
+    ds = from_object(ds_in)
+    assert len(ds) == 2
+    assert np.allclose(ds.metadata("level"), coords["level"])
+
+    for ll in [ds[0].to_latlon(), ds.to_latlon()]:
+        assert ll["lat"].shape == (9,)
+        assert ll["lon"].shape == (9,)
+        assert np.allclose(ll["lat"], coords[lat_name][1])
+        assert np.allclose(ll["lon"], coords[lon_name][1])
+
+
+@pytest.mark.parametrize("lat_name,lon_name", [("lat", "lon"), ("latitude", "longitude")])
+def test_netcdf_geography_1d_2(lat_name, lon_name):
+    # Dimensions:  (level: 2, values: 9)
+    # Coordinates:
+    # * level    (level) int64 16B 700 500
+    # Dimensions without coordinates: values
+    # Data variables:
+    #     a        (level, values) int64 144B 11 12 13 21 22 23 ... 24 25 26 34 35 36
+    #     lat      (values) int64 72B 50 50 50 40 40 40 30 30 30
+    #     lon      (values) int64 72B 0 10 20 0 10 20 0 10 20
+
+    import xarray as xr
+
+    dims = {"level": 2, "values": 9}
+    coords = {
+        "level": np.array([700, 500]),
+    }
+
+    data = np.array(
+        [
+            [11, 12, 13, 21, 22, 23, 31, 32, 33],
+            [14, 15, 16, 24, 25, 26, 34, 35, 36],
+        ]
+    )
+
+    a = xr.Variable(dims, data)
+    lat = xr.Variable({"values": 9}, np.array([50, 50, 50, 40, 40, 40, 30, 30, 30]))
+    lon = xr.Variable({"values": 9}, np.array([0, 10, 20, 0, 10, 20, 0, 10, 20]))
+
+    v = {"a": a, lat_name: lat, lon_name: lon}
+    ds_in = xr.Dataset(v, coords=coords)
+
+    ds = from_object(ds_in)
+    assert len(ds) == 2
+    assert np.allclose(ds.metadata("level"), coords["level"])
+
+    for ll in [ds[0].to_latlon(), ds.to_latlon()]:
+        assert ll["lat"].shape == (9,)
+        assert ll["lon"].shape == (9,)
+        assert np.allclose(ll["lat"], lat.data)
+        assert np.allclose(ll["lon"], lon.data)
 
 
 if __name__ == "__main__":
diff --git a/tests/netcdf/test_netcdf_metadata.py b/tests/netcdf/test_netcdf_metadata.py
index 81cd3aa..2225719 100644
--- a/tests/netcdf/test_netcdf_metadata.py
+++ b/tests/netcdf/test_netcdf_metadata.py
@@ -15,6 +15,7 @@ import pytest
 
 from earthkit.data import from_source
 from earthkit.data.testing import earthkit_examples_file
+from earthkit.data.testing import earthkit_remote_test_data_file
 from earthkit.data.testing import load_nc_or_xr_source
 
 
@@ -81,6 +82,14 @@ def test_netcdf_valid_datetime(mode):
     assert ds[0].metadata("valid_datetime") == "2020-05-13T12:00:00"
 
 
+def test_netcdf_forecast_reference_time():
+    ds = from_source("url", earthkit_remote_test_data_file("test-data", "fa_ta850.nc"))
+
+    assert len(ds) == 37
+    assert ds[0].metadata("valid_datetime") == "2020-01-23T00:00:00"
+    assert ds[5].metadata("valid_datetime") == "2020-01-23T05:00:00"
+
+
 if __name__ == "__main__":
     from earthkit.data.testing import main