Fix interpolating fieldlists with regular latlon grid (#32)

* Fix regular_ll fieldlist interpolation

docs/release_notes/version_0.3_updates.rst

@@ -2,12 +2,20 @@ Version 0.3 Updates
 /////////////////////////
 
 
+Version 0.3.3
+===============
+
+Fixes
+++++++++++++++++
+- fixed issue when failed to interpolate earthkit-data :xref:`fieldlist`\ s containing a regular latitude-longitude grid
+
+
 Version 0.3.2
 ===============
 
 Fixes
 ++++++++++++++++
-- fixed issue when  reading an interpolation matrix from the cache unnecessarily invoked checking of remote matrix files on download server
+- fixed issue when reading an interpolation matrix from the cache unnecessarily invoked checking of remote matrix files on download server
 
 
 Version 0.3.1

earthkit/regrid/db.py

@@ -11,7 +11,6 @@
 import logging
 import os
 from abc import ABCMeta, abstractmethod
-from functools import lru_cache
 
 from scipy.sparse import load_npz
 
@@ -64,15 +63,20 @@ def matrix_path(self, name):
 class UrlAccessor(MatrixAccessor):
     def __init__(self, url):
         self._url = url
+        self._index_path = None
 
     def path(self):
         return self._url
 
     def is_local(self):
         False
 
-    @lru_cache
     def index_path(self):
+        if self._index_path is None or not os.path.exists(self._index_path):
+            self._index_path = self.__index_path()
+        return self._index_path
+
+    def __index_path(self):
         from earthkit.regrid.utils.caching import cache_file
 
         url = os.path.join(self._url, _INDEX_FILENAME)
@@ -279,12 +283,11 @@ def find(self, gridspec_in, gridspec_out, method):
         gridspec_in = GridSpec.from_dict(gridspec_in)
         gridspec_out = GridSpec.from_dict(gridspec_out)
 
-        # print(f"{gridspec_in=} {gridspec_out=} {method=}")
         if gridspec_in is None or gridspec_out is None:
             return None
 
         for _, entry in self.items():
-            if self.match(entry, gridspec_in, gridspec_out, method):
+            if MatrixIndex.match(entry, gridspec_in, gridspec_out, method):
                 return entry
         return None
 

earthkit/regrid/gridspec.py

@@ -10,6 +10,8 @@
 import logging
 import re
 
+import numpy as np
+
 LOG = logging.getLogger(__name__)
 
 FULL_GLOBE = 360.0
@@ -77,14 +79,29 @@ def __eq__(self, o):
             return False
 
         for k in self.COMPARE_KEYS:
-            if str(self[k]) != str(o[k]):
+            if not self.compare_key(k, self[k], o[k]):
                 return False
 
         if "shape" in self and "shape" in o:
             if self["shape"] != o["shape"]:
                 return False
         return True
 
+    @staticmethod
+    def compare_key(key, v1, v2):
+        if isinstance(v1, str) and isinstance(v2, str):
+            return v1 == v2
+        elif key == "grid" and isinstance(v1, list) and isinstance(v2, list):
+            return np.allclose(np.array(v1), np.array(v2), atol=1e-6)
+        elif isinstance(v1, list) and isinstance(v2, list):
+            return v1 == v2
+        elif isinstance(v1, float) and isinstance(v2, float):
+            return np.isclose(v1, v2)
+        elif isinstance(v1, int) and isinstance(v2, int):
+            return v1 == v2
+        else:
+            return str(v1) == str(v2)
+
     @staticmethod
     def same_area(area1, area2, eps=DEGREE_EPS):
         if len(area1) == len(area2):

tests/test_fieldlist.py

@@ -0,0 +1,79 @@
+# (C) Copyright 2023 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+
+import os
+import sys
+
+import numpy as np
+import pytest
+
+from earthkit.regrid import interpolate
+
+here = os.path.dirname(__file__)
+sys.path.insert(0, here)
+from testing import get_test_data  # noqa: E402
+from testing import get_test_data_path  # noqa: E402
+
+try:
+    from earthkit.data import from_source  # noqa
+
+    NO_EKD = False
+except Exception:
+    NO_EKD = True
+
+
+@pytest.mark.download
+@pytest.mark.tmp_cache
+@pytest.mark.skipif(NO_EKD, reason="No access to earthkit-data")
+@pytest.mark.parametrize(
+    "_kwarg,method",
+    [
+        ({}, "linear"),
+        ({"method": "linear"}, "linear"),
+        ({"method": "nearest-neighbour"}, "nearest-neighbour"),
+        ({"method": "nn"}, "nearest-neighbour"),
+        ({"method": "nearest-neighbor"}, "nearest-neighbour"),
+    ],
+)
+def test_fieldlist_reg_ll(_kwarg, method):
+    ds = from_source("url", get_test_data_path("5x5.grib"))
+
+    f_ref = get_test_data(f"out_5x5_10x10_{method}.npz")
+    v_ref = np.load(f_ref)["arr_0"]
+
+    r = interpolate(ds, out_grid={"grid": [10, 10]}, **_kwarg)
+
+    assert len(r) == 1
+    assert r[0].shape == (19, 36)
+    assert np.allclose(r[0].values, v_ref)
+
+
+@pytest.mark.download
+@pytest.mark.tmp_cache
+@pytest.mark.skipif(NO_EKD, reason="No access to earthkit-data")
+@pytest.mark.parametrize(
+    "_kwarg,method",
+    [
+        ({}, "linear"),
+        ({"method": "linear"}, "linear"),
+        ({"method": "nearest-neighbour"}, "nearest-neighbour"),
+        ({"method": "nn"}, "nearest-neighbour"),
+        ({"method": "nearest-neighbor"}, "nearest-neighbour"),
+    ],
+)
+def test_fieldlist_gg(_kwarg, method):
+    ds = from_source("url", get_test_data_path("O32.grib"))
+
+    f_ref = get_test_data(f"out_O32_10x10_{method}.npz")
+    v_ref = np.load(f_ref)["arr_0"]
+
+    r = interpolate(ds, out_grid={"grid": [10, 10]}, **_kwarg)
+
+    assert len(r) == 1
+    assert r[0].shape == (19, 36)
+    assert np.allclose(r[0].values, v_ref)

tests/test_interpolate.py

@@ -7,42 +7,16 @@
 # nor does it submit to any jurisdiction.
 
 import os
+import sys
 
 import numpy as np
 import pytest
 
 from earthkit.regrid import interpolate
 
-PATH = os.path.dirname(__file__)
-
-URL_ROOT = "https://get.ecmwf.int/repository/test-data/earthkit-regrid/test-data"
-
-
-def simple_download(url, target):
-    import requests
-
-    r = requests.get(url, allow_redirects=True)
-    r.raise_for_status()
-    open(target, "wb").write(r.content)
-
-
-def get_test_data(filename, subfolder="global_0_360"):
-    if not isinstance(filename, list):
-        filename = [filename]
-
-    res = []
-    for fn in filename:
-        d_path = os.path.join(PATH, "data", subfolder)
-        os.makedirs(d_path, exist_ok=True)
-        f_path = os.path.join(d_path, fn)
-        if not os.path.exists(f_path):
-            simple_download(url=f"{URL_ROOT}/{subfolder}/{fn}", target=f_path)
-        res.append(f_path)
-
-    if len(res) == 1:
-        return res[0]
-    else:
-        return res
+here = os.path.dirname(__file__)
+sys.path.insert(0, here)
+from testing import get_test_data  # noqa: E402
 
 
 @pytest.mark.download

tests/testing.py

@@ -0,0 +1,44 @@
+# (C) Copyright 2023 ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+
+import os
+
+PATH = os.path.dirname(__file__)
+
+URL_ROOT = "https://get.ecmwf.int/repository/test-data/earthkit-regrid/test-data"
+
+
+def simple_download(url, target):
+    import requests
+
+    r = requests.get(url, allow_redirects=True)
+    r.raise_for_status()
+    open(target, "wb").write(r.content)
+
+
+def get_test_data_path(filename, subfolder="global_0_360"):
+    return os.path.join(URL_ROOT, subfolder, filename)
+
+
+def get_test_data(filename, subfolder="global_0_360"):
+    if not isinstance(filename, list):
+        filename = [filename]
+
+    res = []
+    for fn in filename:
+        d_path = os.path.join(PATH, "data", subfolder)
+        os.makedirs(d_path, exist_ok=True)
+        f_path = os.path.join(d_path, fn)
+        if not os.path.exists(f_path):
+            simple_download(url=f"{URL_ROOT}/{subfolder}/{fn}", target=f_path)
+        res.append(f_path)
+
+    if len(res) == 1:
+        return res[0]
+    else:
+        return res