From 70c62531dbd6ff1618cfc8768ddb10ea129629a7 Mon Sep 17 00:00:00 2001
From: tomvothecoder <tomvothecoder@gmail.com>
Date: Wed, 3 Dec 2025 12:02:02 -0800
Subject: [PATCH 01/12] Initial code to support multidimensional coordinates

---
 xcdat/axis.py               | 12 ++++++++++++
 xcdat/regridder/accessor.py | 28 ++++++++++++++++++++--------
 2 files changed, 32 insertions(+), 8 deletions(-)

diff --git a/xcdat/axis.py b/xcdat/axis.py
index 8cab0657..44a9391f 100644
--- a/xcdat/axis.py
+++ b/xcdat/axis.py
@@ -5,6 +5,7 @@
 
 from typing import Literal
 
+import cf_xarray as cfxr  # noqa: F401
 import numpy as np
 import xarray as xr
 
@@ -122,6 +123,17 @@ def get_dim_coords(
     # Example: ["lat", "lon", "time"]
     index_keys = obj.indexes.keys()
 
+    # FIXME: Attempt -- If index_keys is None, attempt to retrieve keys using
+    # cf_xarray for the axis. This is a fallback for objects with
+    # multidimensional coordiantes.
+    if not index_keys:
+        try:
+            obj = obj.cf[axis]
+        except KeyError:
+            raise KeyError(
+                f"Could not find coordinate for dimension '{axis}'"
+            ) from None
+
     # Attempt to map the axis it all of its coordinate variable(s) using the
     # axis and coordinate names in the object attributes (if they are set).
     # Example: Returns ["time", "time_centered"] with `axis="T"`
diff --git a/xcdat/regridder/accessor.py b/xcdat/regridder/accessor.py
index 40062aed..82da0fd4 100644
--- a/xcdat/regridder/accessor.py
+++ b/xcdat/regridder/accessor.py
@@ -166,6 +166,7 @@ def horizontal(
                 f"Tool {e!s} does not exist, valid choices {list(HORIZONTAL_REGRID_TOOLS)}"
             ) from e
 
+        # FIXME: Line 170 -- Can't add bounds for multidimensional coordinates
         input_grid = _get_input_grid(self._ds, data_var, ["X", "Y"])
         regridder = regrid_tool(input_grid, output_grid, **options)
         output_ds = regridder.horizontal(data_var, self._ds)
@@ -236,13 +237,8 @@ def vertical(
                 f"Tool {e!s} does not exist, valid choices "
                 f"{list(VERTICAL_REGRID_TOOLS)}"
             ) from e
-        input_grid = _get_input_grid(
-            self._ds,
-            data_var,
-            [
-                "Z",
-            ],
-        )
+
+        input_grid = _get_input_grid(self._ds, data_var, ["Z"])
         regridder = regrid_tool(input_grid, output_grid, **options)
         output_ds = regridder.vertical(data_var, self._ds)
 
@@ -310,8 +306,10 @@ def _obj_to_grid_ds(obj: xr.Dataset | xr.DataArray) -> xr.Dataset:
     # same axis (e.g., (nlat, lat) for X and (nlon, lon) for Y). We only
     # need lat_bnds and lon_bnds for the X and Y axes, respectively, and not
     # nlat_bnds and nlon_bnds.
+
     for axis, has_bounds in axis_has_bounds.items():
         if not has_bounds:
+            # FIXME: Line 313 --Can't add bounds for multidimensional coordinates
             output_ds = output_ds.bounds.add_bounds(axis=axis)
 
     return output_ds
@@ -374,7 +372,20 @@ def _get_input_grid(ds: xr.Dataset, data_var: str, dup_check_dims: list[CFAxisKe
     all_coords = set(ds.coords.keys())
 
     for dimension in dup_check_dims:
-        coords = get_dim_coords(ds, dimension)
+        try:
+            coords = get_dim_coords(ds, dimension)
+        except KeyError:
+            coords = None
+
+        # If index_keys is None, attempt to retrieve keys using cf_xarray for the axis.
+        # This is a fallback for objects with multidimensional coordiantes
+        if not coords:
+            try:
+                coords = ds.cf.axes[dimension]
+            except KeyError:
+                raise KeyError(
+                    f"Could not find coordinate for dimension '{dimension}'"
+                ) from None
 
         if isinstance(coords, xr.Dataset):
             coord = set([get_dim_coords(ds[data_var], dimension).name])
@@ -387,6 +398,7 @@ def _get_input_grid(ds: xr.Dataset, data_var: str, dup_check_dims: list[CFAxisKe
     input_grid = ds.drop_dims(to_drop)
 
     # drops extra dimensions on input grid
+    # FIXME: Line 402 --Can't add bounds for multidimensional coordinates
     grid = input_grid.regridder.grid
 
     # preserve mask on grid

From 717df67064d45c33fa855ea6e40b408c115fbaed Mon Sep 17 00:00:00 2001
From: tomvothecoder <tomvothecoder@gmail.com>
Date: Wed, 3 Dec 2025 12:03:57 -0800
Subject: [PATCH 02/12] Remove unused code

---
 xcdat/regridder/accessor.py | 10 ----------
 1 file changed, 10 deletions(-)

diff --git a/xcdat/regridder/accessor.py b/xcdat/regridder/accessor.py
index 82da0fd4..1203994b 100644
--- a/xcdat/regridder/accessor.py
+++ b/xcdat/regridder/accessor.py
@@ -377,16 +377,6 @@ def _get_input_grid(ds: xr.Dataset, data_var: str, dup_check_dims: list[CFAxisKe
         except KeyError:
             coords = None
 
-        # If index_keys is None, attempt to retrieve keys using cf_xarray for the axis.
-        # This is a fallback for objects with multidimensional coordiantes
-        if not coords:
-            try:
-                coords = ds.cf.axes[dimension]
-            except KeyError:
-                raise KeyError(
-                    f"Could not find coordinate for dimension '{dimension}'"
-                ) from None
-
         if isinstance(coords, xr.Dataset):
             coord = set([get_dim_coords(ds[data_var], dimension).name])
 

From 00635c202133eb6c0404b54adcf70465391078c2 Mon Sep 17 00:00:00 2001
From: tomvothecoder <tomvothecoder@gmail.com>
Date: Wed, 3 Dec 2025 12:06:35 -0800
Subject: [PATCH 03/12] Remove unused code

---
 xcdat/regridder/accessor.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/xcdat/regridder/accessor.py b/xcdat/regridder/accessor.py
index 1203994b..5c092198 100644
--- a/xcdat/regridder/accessor.py
+++ b/xcdat/regridder/accessor.py
@@ -372,10 +372,7 @@ def _get_input_grid(ds: xr.Dataset, data_var: str, dup_check_dims: list[CFAxisKe
     all_coords = set(ds.coords.keys())
 
     for dimension in dup_check_dims:
-        try:
-            coords = get_dim_coords(ds, dimension)
-        except KeyError:
-            coords = None
+        coords = get_dim_coords(ds, dimension)
 
         if isinstance(coords, xr.Dataset):
             coord = set([get_dim_coords(ds[data_var], dimension).name])

From 048a07ae228d8921b754e06e1f33048374f1ac8d Mon Sep 17 00:00:00 2001
From: Tom Vo <tomvothecoder@gmail.com>
Date: Wed, 3 Dec 2025 12:12:50 -0800
Subject: [PATCH 04/12] Apply suggestion from @tomvothecoder

---
 xcdat/regridder/accessor.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/xcdat/regridder/accessor.py b/xcdat/regridder/accessor.py
index 5c092198..84855154 100644
--- a/xcdat/regridder/accessor.py
+++ b/xcdat/regridder/accessor.py
@@ -385,7 +385,7 @@ def _get_input_grid(ds: xr.Dataset, data_var: str, dup_check_dims: list[CFAxisKe
     input_grid = ds.drop_dims(to_drop)
 
     # drops extra dimensions on input grid
-    # FIXME: Line 402 --Can't add bounds for multidimensional coordinates
+    # FIXME: Line 389 --Can't add bounds for multidimensional coordinates
     grid = input_grid.regridder.grid
 
     # preserve mask on grid

From bc35ab9e2bab223ea3445f944aad2e272bf68502 Mon Sep 17 00:00:00 2001
From: Jason Boutte <boutte3@llnl.gov>
Date: Fri, 20 Mar 2026 16:26:25 -0700
Subject: [PATCH 05/12] Adds support for multidimensional coordinates in
 regridders

---
 xcdat/axis.py               | 25 +++++++++----------------
 xcdat/regridder/accessor.py | 29 +++++++++++++++++++++--------
 xcdat/regridder/base.py     |  6 ++++++
 xcdat/regridder/regrid2.py  |  2 ++
 xcdat/regridder/xesmf.py    |  2 ++
 xcdat/regridder/xgcm.py     |  2 ++
 6 files changed, 42 insertions(+), 24 deletions(-)

diff --git a/xcdat/axis.py b/xcdat/axis.py
index 44a9391f..3e399097 100644
--- a/xcdat/axis.py
+++ b/xcdat/axis.py
@@ -74,7 +74,7 @@ def get_dim_keys(obj: xr.Dataset | xr.DataArray, axis: CFAxisKey) -> str | list[
 
 
 def get_dim_coords(
-    obj: xr.Dataset | xr.DataArray, axis: CFAxisKey
+    obj: xr.Dataset | xr.DataArray, axis: CFAxisKey, multidim: bool = False
 ) -> xr.Dataset | xr.DataArray:
     """Gets the dimension coordinates for an axis.
 
@@ -118,21 +118,14 @@ def get_dim_coords(
     ----------
     .. [1] https://cf-xarray.readthedocs.io/en/latest/coord_axes.html#axes-and-coordinates
     """
-    # Get the object's index keys, with each being a dimension.
-    # NOTE: xarray does not include multidimensional coordinates as index keys.
-    # Example: ["lat", "lon", "time"]
-    index_keys = obj.indexes.keys()
-
-    # FIXME: Attempt -- If index_keys is None, attempt to retrieve keys using
-    # cf_xarray for the axis. This is a fallback for objects with
-    # multidimensional coordiantes.
-    if not index_keys:
-        try:
-            obj = obj.cf[axis]
-        except KeyError:
-            raise KeyError(
-                f"Could not find coordinate for dimension '{axis}'"
-            ) from None
+    if multidim:
+        # multidimensional coordinates cannot be indexes, use all coords
+        index_keys = list(obj.coords)
+    else:
+        # Get the object's index keys, with each being a dimension.
+        # NOTE: xarray does not include multidimensional coordinates as index keys.
+        # Example: ["lat", "lon", "time"]
+        index_keys = obj.indexes.keys()
 
     # Attempt to map the axis it all of its coordinate variable(s) using the
     # axis and coordinate names in the object attributes (if they are set).
diff --git a/xcdat/regridder/accessor.py b/xcdat/regridder/accessor.py
index 84855154..78c639f5 100644
--- a/xcdat/regridder/accessor.py
+++ b/xcdat/regridder/accessor.py
@@ -166,8 +166,9 @@ def horizontal(
                 f"Tool {e!s} does not exist, valid choices {list(HORIZONTAL_REGRID_TOOLS)}"
             ) from e
 
-        # FIXME: Line 170 -- Can't add bounds for multidimensional coordinates
-        input_grid = _get_input_grid(self._ds, data_var, ["X", "Y"])
+        input_grid = _get_input_grid(
+            self._ds, data_var, ["X", "Y"], multidim=regrid_tool.can_handle_multidim()
+        )
         regridder = regrid_tool(input_grid, output_grid, **options)
         output_ds = regridder.horizontal(data_var, self._ds)
 
@@ -245,7 +246,9 @@ def vertical(
         return output_ds
 
 
-def _obj_to_grid_ds(obj: xr.Dataset | xr.DataArray) -> xr.Dataset:
+def _obj_to_grid_ds(
+    obj: xr.Dataset | xr.DataArray, multidim: bool = False
+) -> xr.Dataset:
     """
     Convert an xarray object to a new Dataset containing axis coordinates and
     bounds.
@@ -300,6 +303,10 @@ def _obj_to_grid_ds(obj: xr.Dataset | xr.DataArray) -> xr.Dataset:
         attrs=obj.attrs,
     )
 
+    # Multidimensional coordinates bounds generation is not supported
+    if multidim:
+        return output_ds
+
     # Add bounds only for axes that do not already have them. This
     # prevents multiple sets of bounds being added for the same axis.
     # For example, curvilinear grids can have multiple coordinates for the
@@ -345,7 +352,12 @@ def _get_axis_coord_and_bounds(
     return coord_var, bounds_var
 
 
-def _get_input_grid(ds: xr.Dataset, data_var: str, dup_check_dims: list[CFAxisKey]):
+def _get_input_grid(
+    ds: xr.Dataset,
+    data_var: str,
+    dup_check_dims: list[CFAxisKey],
+    multidim: bool = False,
+):
     """
     Extract the grid from ``ds``.
 
@@ -372,10 +384,12 @@ def _get_input_grid(ds: xr.Dataset, data_var: str, dup_check_dims: list[CFAxisKe
     all_coords = set(ds.coords.keys())
 
     for dimension in dup_check_dims:
-        coords = get_dim_coords(ds, dimension)
+        coords = get_dim_coords(ds, dimension, multidim=multidim)
 
         if isinstance(coords, xr.Dataset):
-            coord = set([get_dim_coords(ds[data_var], dimension).name])
+            coord = set(
+                [get_dim_coords(ds[data_var], dimension, multidim=multidim).name]
+            )
 
             dimension_coords = set(ds.cf[[dimension]].coords.keys())
 
@@ -385,8 +399,7 @@ def _get_input_grid(ds: xr.Dataset, data_var: str, dup_check_dims: list[CFAxisKe
     input_grid = ds.drop_dims(to_drop)
 
     # drops extra dimensions on input grid
-    # FIXME: Line 389 --Can't add bounds for multidimensional coordinates
-    grid = input_grid.regridder.grid
+    grid = _obj_to_grid_ds(input_grid, multidim=multidim)
 
     # preserve mask on grid
     if "mask" in ds:
diff --git a/xcdat/regridder/base.py b/xcdat/regridder/base.py
index 0458abad..36871f64 100644
--- a/xcdat/regridder/base.py
+++ b/xcdat/regridder/base.py
@@ -95,6 +95,12 @@ def _drop_axis(ds: xr.Dataset, axis: list[CFAxisKey]) -> xr.Dataset:
 class BaseRegridder(abc.ABC):
     """BaseRegridder."""
 
+    supports_multidim: bool
+
+    @classmethod
+    def can_handle_multidim(cls) -> bool:
+        return cls.supports_multidim
+
     def __init__(self, input_grid: xr.Dataset, output_grid: xr.Dataset, **options: Any):
         self._input_grid = input_grid
         self._output_grid = output_grid
diff --git a/xcdat/regridder/regrid2.py b/xcdat/regridder/regrid2.py
index bb702c4b..dfa90eab 100644
--- a/xcdat/regridder/regrid2.py
+++ b/xcdat/regridder/regrid2.py
@@ -11,6 +11,8 @@
 
 
 class Regrid2Regridder(BaseRegridder):
+    supports_multidim = False
+
     def __init__(
         self,
         input_grid: xr.Dataset,
diff --git a/xcdat/regridder/xesmf.py b/xcdat/regridder/xesmf.py
index f4fa8961..325da558 100644
--- a/xcdat/regridder/xesmf.py
+++ b/xcdat/regridder/xesmf.py
@@ -19,6 +19,8 @@
 
 
 class XESMFRegridder(BaseRegridder):
+    supports_multidim = True
+
     def __init__(
         self,
         input_grid: xr.Dataset,
diff --git a/xcdat/regridder/xgcm.py b/xcdat/regridder/xgcm.py
index c0619755..337c4788 100644
--- a/xcdat/regridder/xgcm.py
+++ b/xcdat/regridder/xgcm.py
@@ -13,6 +13,8 @@
 
 
 class XGCMRegridder(BaseRegridder):
+    supports_multidim = False
+
     def __init__(
         self,
         input_grid: xr.Dataset,

From 6399d55162be94f554e6f524f675c605b71c71c2 Mon Sep 17 00:00:00 2001
From: Jason Boutte <boutte3@llnl.gov>
Date: Fri, 3 Apr 2026 14:48:13 -0700
Subject: [PATCH 06/12] Fix mypy error

---
 xcdat/regridder/accessor.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/xcdat/regridder/accessor.py b/xcdat/regridder/accessor.py
index 78c639f5..25802146 100644
--- a/xcdat/regridder/accessor.py
+++ b/xcdat/regridder/accessor.py
@@ -5,16 +5,17 @@
 from xcdat.axis import CFAxisKey, get_coords_by_name, get_dim_coords
 from xcdat.bounds import create_bounds
 from xcdat.regridder import regrid2, xesmf, xgcm
+from xcdat.regridder.base import BaseRegridder
 from xcdat.regridder.grid import _validate_grid_has_single_axis_dim
 
 HorizontalRegridTools = Literal["xesmf", "regrid2"]
-HORIZONTAL_REGRID_TOOLS = {
+HORIZONTAL_REGRID_TOOLS: dict[str, type[BaseRegridder]] = {
     "regrid2": regrid2.Regrid2Regridder,
     "xesmf": xesmf.XESMFRegridder,
 }
 
 VerticalRegridTools = Literal["xgcm"]
-VERTICAL_REGRID_TOOLS = {"xgcm": xgcm.XGCMRegridder}
+VERTICAL_REGRID_TOOLS: dict[str, type[BaseRegridder]] = {"xgcm": xgcm.XGCMRegridder}
 
 
 @xr.register_dataset_accessor(name="regridder")

From d13a72ab2c09f00c2b704163d018922aa41a8f7f Mon Sep 17 00:00:00 2001
From: Jason Boutte <boutte3@llnl.gov>
Date: Fri, 3 Apr 2026 14:49:21 -0700
Subject: [PATCH 07/12] Fixes where coordinates are sourced when allowing
 multi-dim coords

---
 xcdat/axis.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/xcdat/axis.py b/xcdat/axis.py
index 3e399097..cc4e12d2 100644
--- a/xcdat/axis.py
+++ b/xcdat/axis.py
@@ -120,12 +120,12 @@ def get_dim_coords(
     """
     if multidim:
         # multidimensional coordinates cannot be indexes, use all coords
-        index_keys = list(obj.coords)
+        index_keys = list([y for x in obj.cf.coordinates.values() for y in x])
     else:
         # Get the object's index keys, with each being a dimension.
         # NOTE: xarray does not include multidimensional coordinates as index keys.
         # Example: ["lat", "lon", "time"]
-        index_keys = obj.indexes.keys()
+        index_keys = list(obj.indexes.keys())
 
     # Attempt to map the axis it all of its coordinate variable(s) using the
     # axis and coordinate names in the object attributes (if they are set).

From f9272262d4547b7b3adc9ef1ce63bc4c0ba6b5d5 Mon Sep 17 00:00:00 2001
From: Jason Boutte <boutte3@llnl.gov>
Date: Fri, 1 May 2026 11:22:06 -0700
Subject: [PATCH 08/12] Removes old comment

---
 xcdat/regridder/accessor.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/xcdat/regridder/accessor.py b/xcdat/regridder/accessor.py
index 25802146..060744ba 100644
--- a/xcdat/regridder/accessor.py
+++ b/xcdat/regridder/accessor.py
@@ -317,7 +317,6 @@ def _obj_to_grid_ds(
 
     for axis, has_bounds in axis_has_bounds.items():
         if not has_bounds:
-            # FIXME: Line 313 --Can't add bounds for multidimensional coordinates
             output_ds = output_ds.bounds.add_bounds(axis=axis)
 
     return output_ds

From 24d2c25ecc0ecc9252c80eb43a8f6218d940d8f1 Mon Sep 17 00:00:00 2001
From: Jason Boutte <boutte3@llnl.gov>
Date: Fri, 1 May 2026 11:28:30 -0700
Subject: [PATCH 09/12] fix: default supports_multidim to False on
 BaseRegridder

Set a default value for supports_multidim so subclasses are not required
to explicitly declare it. Added tests to verify the default and override
behavior.
---
 tests/test_regrid.py    | 22 ++++++++++++++++++++++
 xcdat/regridder/base.py |  2 +-
 2 files changed, 23 insertions(+), 1 deletion(-)

diff --git a/tests/test_regrid.py b/tests/test_regrid.py
index ab9b619e..f7b5e8d4 100644
--- a/tests/test_regrid.py
+++ b/tests/test_regrid.py
@@ -1677,3 +1677,25 @@ def horizontal(self, data_var, ds):
         ds_out = regridder.horizontal("ts", ds_in)
 
         assert ds_in == ds_out
+
+    def test_supports_multidim_defaults_to_false(self):
+        """Test that BaseRegridder subclasses default supports_multidim to False."""
+
+        class MinimalRegridder(base.BaseRegridder):
+            def horizontal(self, data_var, ds):
+                return ds
+
+        assert MinimalRegridder.supports_multidim is False
+        assert MinimalRegridder.can_handle_multidim() is False
+
+    def test_supports_multidim_can_be_overridden(self):
+        """Test that subclasses can override supports_multidim to True."""
+
+        class MultidimRegridder(base.BaseRegridder):
+            supports_multidim = True
+
+            def horizontal(self, data_var, ds):
+                return ds
+
+        assert MultidimRegridder.supports_multidim is True
+        assert MultidimRegridder.can_handle_multidim() is True
diff --git a/xcdat/regridder/base.py b/xcdat/regridder/base.py
index 36871f64..571e1255 100644
--- a/xcdat/regridder/base.py
+++ b/xcdat/regridder/base.py
@@ -95,7 +95,7 @@ def _drop_axis(ds: xr.Dataset, axis: list[CFAxisKey]) -> xr.Dataset:
 class BaseRegridder(abc.ABC):
     """BaseRegridder."""
 
-    supports_multidim: bool
+    supports_multidim: bool = False
 
     @classmethod
     def can_handle_multidim(cls) -> bool:

From eb9cc42f32a30e20b2173705d9355c122f05d302 Mon Sep 17 00:00:00 2001
From: Jason Boutte <boutte3@llnl.gov>
Date: Fri, 1 May 2026 11:30:24 -0700
Subject: [PATCH 10/12] fix: include cf axes in multidim coord discovery in
 get_dim_coords

The multidim=True branch only used obj.cf.coordinates to build
index_keys, missing coordinates discoverable only via obj.cf.axes.
Combine both sources to avoid incorrect KeyError when an axis
coordinate exists but is not in cf.coordinates.
---
 tests/test_axis.py | 26 ++++++++++++++++++++++++++
 xcdat/axis.py      | 11 +++++++++--
 2 files changed, 35 insertions(+), 2 deletions(-)

diff --git a/tests/test_axis.py b/tests/test_axis.py
index b85fb344..df604f2b 100644
--- a/tests/test_axis.py
+++ b/tests/test_axis.py
@@ -247,6 +247,32 @@ def test_returns_dataarray_dimension_coordinate_var_using_standard_name_attr(sel
 
         assert result.identical(expected)
 
+    def test_multidim_finds_coords_only_in_cf_axes(self):
+        """Test that multidim=True finds coords discoverable via cf axes."""
+        # Create a dataset where 'lat' has axis="Y" attribute but no
+        # standard_name or recognized coordinate name that cf_xarray would
+        # put in obj.cf.coordinates.
+        lat = xr.DataArray(
+            data=np.array([[0, 1], [2, 3]]),
+            dims=["y", "x"],
+            attrs={"axis": "Y"},
+        )
+        lon = xr.DataArray(
+            data=np.array([[10, 11], [12, 13]]),
+            dims=["y", "x"],
+            attrs={"axis": "X"},
+        )
+        ds = xr.Dataset(coords={"lat": lat, "lon": lon})
+
+        # Without multidim, this would raise KeyError because lat/lon are
+        # multidimensional and not in indexes.
+        with pytest.raises(KeyError):
+            get_dim_coords(ds, "Y", multidim=False)
+
+        # With multidim=True, it should find lat via cf axes.
+        result = get_dim_coords(ds, "Y", multidim=True)
+        xr.testing.assert_identical(result, ds["lat"])
+
 
 class TestGetCoordsByName:
     def test_raises_error_if_coordinate_not_found(self):
diff --git a/xcdat/axis.py b/xcdat/axis.py
index cc4e12d2..ff88ebf4 100644
--- a/xcdat/axis.py
+++ b/xcdat/axis.py
@@ -119,8 +119,15 @@ def get_dim_coords(
     .. [1] https://cf-xarray.readthedocs.io/en/latest/coord_axes.html#axes-and-coordinates
     """
     if multidim:
-        # multidimensional coordinates cannot be indexes, use all coords
-        index_keys = list([y for x in obj.cf.coordinates.values() for y in x])
+        # multidimensional coordinates cannot be indexes, use all coords.
+        # Combine both obj.cf.coordinates and obj.cf.axes to avoid missing
+        # coordinates that are only discoverable via CF axis metadata.
+        cf_keys: set[str] = set()
+        for keys in obj.cf.coordinates.values():
+            cf_keys.update(keys)
+        for keys in obj.cf.axes.values():
+            cf_keys.update(keys)
+        index_keys = list(cf_keys)
     else:
         # Get the object's index keys, with each being a dimension.
         # NOTE: xarray does not include multidimensional coordinates as index keys.

From 37cc65786c718cfc5123a0fdac99fca26c41733a Mon Sep 17 00:00:00 2001
From: Jason Boutte <boutte3@llnl.gov>
Date: Fri, 1 May 2026 11:33:55 -0700
Subject: [PATCH 11/12] test: add regression test for #816 multidim coord
 regridding

Exercises ds.regridder.horizontal(..., tool='xesmf') on a dataset with
2D lat/lon coordinates and no 1D index coord vars, ensuring the
multidim code path correctly discovers coordinates via CF axis metadata.
---
 tests/test_regrid.py | 45 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)

diff --git a/tests/test_regrid.py b/tests/test_regrid.py
index f7b5e8d4..2654678e 100644
--- a/tests/test_regrid.py
+++ b/tests/test_regrid.py
@@ -1617,6 +1617,51 @@ def test_vertical_tool_check(self, _get_input_grid):
         ):
             self.ac.vertical("ts", mock_data, tool="dummy", target_data=None)  # type: ignore
 
+    def test_horizontal_with_multidim_coords_issue_816(self):
+        """Regression test for #816: dataset with 2D lat/lon and no 1D coord vars.
+
+        Ensures that ds.regridder.horizontal(..., tool="xesmf") works on
+        datasets shaped like curvilinear/unstructured grids where lat and lon
+        are multidimensional coordinates (not index dimensions).
+        """
+        ny, nx = 10, 20
+
+        lat_2d = np.linspace(-90, 90, ny * nx).reshape(ny, nx)
+        lon_2d = np.linspace(-180, 180, ny * nx).reshape(ny, nx)
+
+        lat = xr.DataArray(
+            data=lat_2d,
+            dims=["y", "x"],
+            attrs={"units": "degrees_north", "axis": "Y"},
+        )
+        lon = xr.DataArray(
+            data=lon_2d,
+            dims=["y", "x"],
+            attrs={"units": "degrees_east", "axis": "X"},
+        )
+
+        ts = xr.DataArray(
+            data=np.random.default_rng(42).random((ny, nx)),
+            dims=["y", "x"],
+            attrs={"units": "K"},
+        )
+
+        ds = xr.Dataset(
+            data_vars={"ts": ts},
+            coords={"lat": lat, "lon": lon},
+        )
+
+        output_grid = grid.create_uniform_grid(-90, 90, 30.0, -180, 180, 60.0)
+
+        # This should not raise KeyError; the multidim path must discover
+        # lat/lon via CF axis metadata.
+        output = ds.regridder.horizontal(
+            "ts", output_grid, tool="xesmf", method="bilinear"
+        )
+
+        assert "ts" in output
+        assert output.ts.dims == ("lat", "lon")
+
 
 class TestBase:
     def test_preserve_bounds(self):

From e0d344cfcb1faf9347a02282efc689ee5f34204e Mon Sep 17 00:00:00 2001
From: Jason Boutte <boutte3@llnl.gov>
Date: Fri, 1 May 2026 11:46:55 -0700
Subject: [PATCH 12/12] fix: plumb multidim flag through
 _get_axis_coord_and_bounds

Datasets with non-standard 2D coordinate names (e.g., nav_lat/nav_lon)
that are only CF-discoverable via axis metadata were failing because
_get_axis_coord_and_bounds called get_dim_coords without multidim=True.
Pass the flag through so the full coord discovery path works.
---
 tests/test_regrid.py        | 45 +++++++++++++++++++++++++++++++++++++
 xcdat/regridder/accessor.py |  6 ++---
 2 files changed, 48 insertions(+), 3 deletions(-)

diff --git a/tests/test_regrid.py b/tests/test_regrid.py
index 2654678e..b0c251fc 100644
--- a/tests/test_regrid.py
+++ b/tests/test_regrid.py
@@ -1662,6 +1662,51 @@ def test_horizontal_with_multidim_coords_issue_816(self):
         assert "ts" in output
         assert output.ts.dims == ("lat", "lon")
 
+    def test_horizontal_with_nonstandard_multidim_coord_names_issue_816(self):
+        """Regression test for #816: multidim coords with non-standard names.
+
+        Ensures that ds.regridder.horizontal(..., tool="xesmf") works when
+        2D coordinates use non-standard names (e.g., nav_lat/nav_lon) that
+        are not in VAR_NAME_MAP but are discoverable via CF axis attributes.
+        """
+        ny, nx = 10, 20
+
+        lat_2d = np.linspace(-90, 90, ny * nx).reshape(ny, nx)
+        lon_2d = np.linspace(-180, 180, ny * nx).reshape(ny, nx)
+
+        nav_lat = xr.DataArray(
+            data=lat_2d,
+            dims=["y", "x"],
+            attrs={"units": "degrees_north", "axis": "Y"},
+        )
+        nav_lon = xr.DataArray(
+            data=lon_2d,
+            dims=["y", "x"],
+            attrs={"units": "degrees_east", "axis": "X"},
+        )
+
+        ts = xr.DataArray(
+            data=np.random.default_rng(42).random((ny, nx)),
+            dims=["y", "x"],
+            attrs={"units": "K"},
+        )
+
+        ds = xr.Dataset(
+            data_vars={"ts": ts},
+            coords={"nav_lat": nav_lat, "nav_lon": nav_lon},
+        )
+
+        output_grid = grid.create_uniform_grid(-90, 90, 30.0, -180, 180, 60.0)
+
+        # This should not raise KeyError; the multidim flag must be plumbed
+        # through _get_axis_coord_and_bounds to get_dim_coords.
+        output = ds.regridder.horizontal(
+            "ts", output_grid, tool="xesmf", method="bilinear"
+        )
+
+        assert "ts" in output
+        assert output.ts.dims == ("lat", "lon")
+
 
 class TestBase:
     def test_preserve_bounds(self):
diff --git a/xcdat/regridder/accessor.py b/xcdat/regridder/accessor.py
index 060744ba..132d8d9d 100644
--- a/xcdat/regridder/accessor.py
+++ b/xcdat/regridder/accessor.py
@@ -286,7 +286,7 @@ def _obj_to_grid_ds(
 
     with xr.set_options(keep_attrs=True):
         for axis in axis_names:
-            coord, bounds = _get_axis_coord_and_bounds(obj, axis)
+            coord, bounds = _get_axis_coord_and_bounds(obj, axis, multidim=multidim)
 
             if coord is not None:
                 axis_coords[str(coord.name)] = coord
@@ -323,7 +323,7 @@ def _obj_to_grid_ds(
 
 
 def _get_axis_coord_and_bounds(
-    obj: xr.Dataset | xr.DataArray, axis: CFAxisKey
+    obj: xr.Dataset | xr.DataArray, axis: CFAxisKey, multidim: bool = False
 ) -> tuple[xr.DataArray | None, xr.DataArray | None]:
     try:
         coord_var = get_coords_by_name(obj, axis)
@@ -333,7 +333,7 @@ def _get_axis_coord_and_bounds(
             )
     except (ValueError, KeyError):
         try:
-            coord_var = get_dim_coords(obj, axis)  # type: ignore
+            coord_var = get_dim_coords(obj, axis, multidim=multidim)  # type: ignore
             _validate_grid_has_single_axis_dim(axis, coord_var)
         except KeyError:
             coord_var = None