From e5ea6c51da2658b5d20c4ea3c4287794e192a498 Mon Sep 17 00:00:00 2001
From: Nabil Freij <nabil.freij@gmail.com>
Date: Thu, 11 Jun 2026 14:32:30 -0700
Subject: [PATCH 1/2] Support single N-D Time and Quantity lookup tables in
 extra_coords

A single N-D table now represents one coordinate varying over N pixel
dimensions, e.g. a 2-D Time table indexed by (raster scan, raster step).
N-D table models expose their inputs in pixel order to match the APE-14
convention assumed by NDCube world-coordinate generation, and slicing
drops integer-sliced axes from a table's axes when the table loses pixel
dimensions. Also fixes a latent Length1Tabular dispatch bug for tables of
shape (1, n).
---
 changelog/941.feature.rst                     |   1 +
 ndcube/extra_coords/extra_coords.py           |  20 ++-
 ndcube/extra_coords/table_coord.py            | 126 ++++++++++++++----
 .../extra_coords/tests/test_extra_coords.py   |  32 +++++
 .../tests/test_lookup_table_coord.py          |  93 +++++++++++++
 5 files changed, 245 insertions(+), 27 deletions(-)
 create mode 100644 changelog/941.feature.rst

diff --git a/changelog/941.feature.rst b/changelog/941.feature.rst
new file mode 100644
index 000000000..b93850f0c
--- /dev/null
+++ b/changelog/941.feature.rst
@@ -0,0 +1 @@
+`~ndcube.extra_coords.TimeTableCoordinate` and `~ndcube.extra_coords.QuantityTableCoordinate` now accept a single N-D table, representing one coordinate that varies over several pixel dimensions. This allows, for example, a 2-D ``Time`` table indexed by (raster scan, raster step) to be attached to an `~ndcube.NDCube` via ``extra_coords.add`` and sliced along either axis.
diff --git a/ndcube/extra_coords/extra_coords.py b/ndcube/extra_coords/extra_coords.py
index 85243597f..5672f7fd3 100644
--- a/ndcube/extra_coords/extra_coords.py
+++ b/ndcube/extra_coords/extra_coords.py
@@ -260,9 +260,17 @@ def mapping(self):
 
         # The mapping is from the array index (position in the list) to the
         # pixel dimensions (numbers in the list)
-        lts = [list([lt[0]] if isinstance(lt[0], Integral) else lt[0]) for lt in self._lookup_tables]
         converter = partial(convert_between_array_and_pixel_axes, naxes=len(self._ndcube.shape))
-        pixel_indicies = [list(converter(np.array(ids))) for ids in lts]
+        pixel_indicies = []
+        for lut_axis, lut in self._lookup_tables:
+            ids = [lut_axis] if isinstance(lut_axis, Integral) else list(lut_axis)
+            pixel_ids = list(converter(np.array(ids)))
+            if getattr(lut, "_model_inputs_are_pixel_ordered", False):
+                # Single N-D tables expose their model inputs in pixel order,
+                # i.e. reversed with respect to the array-ordered axes given
+                # to `add`.
+                pixel_ids = pixel_ids[::-1]
+            pixel_indicies.append(pixel_ids)
         return tuple(reduce(list.__add__, pixel_indicies))
 
     @mapping.setter
@@ -360,7 +368,6 @@ def _getitem_lookup_tables(self, item):
             n_dropped_dims = np.cumsum([isinstance(i, Integral) for i in item])
         for lut_axis, lut in self._lookup_tables:
             lut_axes = (lut_axis,) if not isinstance(lut_axis, tuple) else lut_axis
-            new_lut_axes = tuple(ax - n_dropped_dims[ax] for ax in lut_axes)
             lut_slice = tuple(item[i] for i in lut_axes)
             if isinstance(lut_slice, tuple) and len(lut_slice) == 1:
                 lut_slice = lut_slice[0]
@@ -370,6 +377,13 @@ def _getitem_lookup_tables(self, item):
             if sliced_lut.is_scalar():
                 dropped_tables.add(sliced_lut)
             else:
+                kept_axes = lut_axes
+                if sliced_lut.n_inputs < len(lut_axes):
+                    # The sliced table lost pixel dimensions (e.g. an N-D
+                    # table sliced with an integer), so drop the
+                    # integer-sliced axes from the table's axes.
+                    kept_axes = tuple(ax for ax in lut_axes if not isinstance(item[ax], Integral))
+                new_lut_axes = tuple(ax - n_dropped_dims[ax] for ax in kept_axes)
                 new_lookup_tables.add((new_lut_axes, sliced_lut))
         new_extra_coords = type(self)()
         new_extra_coords._lookup_tables = list(new_lookup_tables)
diff --git a/ndcube/extra_coords/table_coord.py b/ndcube/extra_coords/table_coord.py
index e975393d4..a9b0e6f74 100644
--- a/ndcube/extra_coords/table_coord.py
+++ b/ndcube/extra_coords/table_coord.py
@@ -170,7 +170,7 @@ def _generate_tabular(lookup_table, interpolation='linear', points_unit=u.pix, *
               'method': interpolation,
               **kwargs}
 
-    if len(lookup_table) == 1:
+    if lookup_table.shape == (1,):
         t = Length1Tabular(points, lookup_table, **kwargs)
     else:
         t = TabularND(points, lookup_table, **kwargs)
@@ -282,6 +282,18 @@ def model(self):
         Generate the Astropy Model for this LookupTable.
         """
 
+    @property
+    def _model_inputs_are_pixel_ordered(self):
+        """
+        True when this coordinate's model inputs are in pixel order.
+
+        Single N-D tables span several pixel dimensions with one model. Their
+        model inputs are exposed in pixel order (reversed array order) so the
+        resulting WCS follows the APE-14 convention expected by
+        `~ndcube.NDCube`.
+        """
+        return False
+
     @property
     def wcs(self):
         """
@@ -301,15 +313,19 @@ class QuantityTableCoordinate(BaseTableCoordinate):
     """
     A lookup table up built on `~astropy.units.Quantity`.
 
-    Quantities must be 1-D but more than one can be provided to represent
-    different dimensions of an N-D coordinate.
+    Either a single N-D Quantity, or multiple 1-D Quantities can be provided.
+    A single N-D Quantity represents one physical coordinate which varies over
+    N pixel dimensions. Multiple 1-D Quantities represent the different
+    dimensions of an N-D coordinate, with each table corresponding to one
+    pixel dimension.
 
     Parameters
     ----------
     tables: one or more `~astropy.units.Quantity`
-        The coordinates. Must be 1 dimensionsal. If coordinate system is >1D,
-        multiple 1-D Quantities can be provided representing the different
-        dimensions
+        The coordinates. Either a single Quantity of any dimensionality
+        representing one coordinate varying over that many pixel dimensions,
+        or multiple 1-D Quantities representing the different dimensions of
+        an N-D coordinate system.
 
     names: `str` or `list` of `str`
         Custom names for the components of the QuantityTableCoord. If provided,
@@ -329,10 +345,10 @@ def __init__(self, *tables, names=None, physical_types=None):
             raise u.UnitsError("All tables must have equivalent units.")
         ndim = len(tables)
         dims = np.array([t.ndim for t in tables])
-        if any(dims > 1):
+        if len(tables) > 1 and any(dims > 1):
             raise ValueError(
-                "Currently all tables must be 1-D. If you need >1D support, please "
-                "raise an issue at https://github.con/sunpy/ndcube/issues")
+                "Multiple tables can only be provided if they are all 1-D. "
+                "A single N-D table representing one coordinate is also supported.")
 
         if isinstance(names, str):
             names = [names]
@@ -379,12 +395,26 @@ def _slice_table(self, i, table, item, new_components, whole_slice):
         if self.physical_types:
             new_components["physical_types"].append(self.physical_types[i])
 
+    @property
+    def _single_nd_table(self):
+        return len(self.table) == 1 and self.table[0].ndim > 1
+
     def __getitem__(self, item):
         if isinstance(item, (slice, Integral)):
             item = (item,)
         if not (len(item) == len(self.table) or len(item) == self.table[0].ndim):
             raise ValueError("Can not slice with incorrect length")
 
+        if self._single_nd_table:
+            # A single N-D table represents one world coordinate, so slicing
+            # reduces the table but never splits or drops individual world
+            # components.
+            ret_table = type(self)(self.table[0][item],
+                                   names=self.names,
+                                   physical_types=self.physical_types)
+            ret_table._dropped_world_dimensions = copy.deepcopy(self._dropped_world_dimensions)
+            return ret_table
+
         new_components = defaultdict(list)
         new_components["dropped_world_dimensions"] = copy.deepcopy(self._dropped_world_dimensions)
 
@@ -400,6 +430,8 @@ def __getitem__(self, item):
 
     @property
     def n_inputs(self):
+        if self._single_nd_table:
+            return self.table[0].ndim
         return len(self.table)
 
     def is_scalar(self):
@@ -412,12 +444,22 @@ def frame(self):
         """
         return _generate_generic_frame(len(self.table), self.unit, self.names, self.physical_types)
 
+    @property
+    def _model_inputs_are_pixel_ordered(self):
+        # Docstring inherited.
+        return self._single_nd_table
+
     @property
     def model(self):
         """
         Generate the Astropy Model for this LookupTable.
         """
-        return _model_from_quantity(self.table, True)
+        model = _model_from_quantity(self.table, True)
+        if self._single_nd_table:
+            # Expose the inputs of the N-D table in pixel order, i.e. reversed
+            # with respect to the table's (array-ordered) dimensions.
+            model = models.Mapping(tuple(range(model.n_inputs))[::-1]) | model
+        return model
 
     @property
     def ndim(self):
@@ -427,6 +469,8 @@ def ndim(self):
         Note this may be different from the number of the dimensions in the
         underlying table(s) if different tables represent different dimensions.
         """
+        if self._single_nd_table:
+            return self.table[0].ndim
         return len(self.table)
 
     @property
@@ -437,6 +481,8 @@ def shape(self):
         Note this may be different from the shape of the underlying table(s)
         if different tables represent a different dimensions.
         """
+        if self._single_nd_table:
+            return self.table[0].shape
         return tuple(len(t) for t in self.table)
 
     def interpolate(self, *new_array_grids, **kwargs):
@@ -472,10 +518,16 @@ def interpolate(self, *new_array_grids, **kwargs):
             raise ValueError("New array grids must all be same shape.")
         # Build array grids for non-interpolated table.
         old_array_grids = tuple(np.arange(d) for d in self.shape)
-        # Iterate through tables and interpolate each.
-        new_tables = [
-            np.interp(new_grid, old_grid, t.value, **kwargs) * t.unit
-            for new_grid, old_grid, t in zip(new_array_grids, old_array_grids, self.table)]
+        if self._single_nd_table:
+            table = self.table[0]
+            new_values = scipy.interpolate.interpn(
+                old_array_grids, table.value, np.stack(new_array_grids, axis=-1), **kwargs)
+            new_tables = [new_values * table.unit]
+        else:
+            # Iterate through tables and interpolate each.
+            new_tables = [
+                np.interp(new_grid, old_grid, t.value, **kwargs) * t.unit
+                for new_grid, old_grid, t in zip(new_array_grids, old_array_grids, self.table)]
         # Rebuild return interpolated coord.
         new_coord = type(self)(*new_tables, names=self.names, physical_types=self.physical_types)
         new_coord._dropped_world_dimensions = self._dropped_world_dimensions
@@ -699,12 +751,16 @@ def interpolate(self, *new_array_grids, mesh_output=None, **kwargs):
 
 class TimeTableCoordinate(BaseTableCoordinate):
     """
-    A lookup table based on a `~astropy.time.Time`, will always be one dimensional.
+    A lookup table based on a `~astropy.time.Time`.
+
+    The table represents a single time coordinate which can vary over one or
+    more pixel dimensions, i.e. the input `~astropy.time.Time` can be N-D.
 
     Parameters
     ----------
     table: `~astropy.time.Time`
-        Time coordinates. Only one can be provided and must be 1D.
+        Time coordinates. Only one can be provided. An N-D table corresponds
+        to N pixel dimensions.
 
     names: `str` or `list` of `str`
         Custom names for the components of the SkyCoord. If provided, a name must
@@ -735,11 +791,15 @@ def __init__(self, *tables, names=None, physical_types=None, reference_time=None
 
         super().__init__(*tables, mesh=False, names=names, physical_types=physical_types)
         self.table = self.table[0]
-        self.reference_time = reference_time or self.table[0]
+        self.reference_time = reference_time or self.table.ravel()[0]
 
     def __getitem__(self, item):
-        if not (isinstance(item, (slice, Integral)) or len(item) == 1):
+        if isinstance(item, (slice, Integral)):
+            item = (item,)
+        if len(item) != max(self.table.ndim, 1):
             raise ValueError("Can not slice with incorrect length")
+        if len(item) == 1:
+            item = item[0]
 
         return type(self)(self.table[item],
                           names=self.names,
@@ -748,7 +808,7 @@ def __getitem__(self, item):
 
     @property
     def n_inputs(self):
-        return 1  # The time table has to be one dimensional
+        return max(self.table.ndim, 1)
 
     def is_scalar(self):
         return self.table.shape == ()
@@ -763,6 +823,11 @@ def frame(self):
                                 axes_names=self.names,
                                 name="TemporalFrame")
 
+    @property
+    def _model_inputs_are_pixel_ordered(self):
+        # Docstring inherited.
+        return self.table.ndim > 1
+
     @property
     def model(self):
         """
@@ -771,9 +836,14 @@ def model(self):
         time = self.table
         deltas = (time - self.reference_time).to(u.s)
 
-        return _model_from_quantity((deltas,), mesh=False)
+        model = _model_from_quantity((deltas,), mesh=False)
+        if deltas.ndim > 1:
+            # Expose the inputs of the N-D table in pixel order, i.e. reversed
+            # with respect to the table's (array-ordered) dimensions.
+            model = models.Mapping(tuple(range(model.n_inputs))[::-1]) | model
+        return model
 
-    def interpolate(self, new_array_grids, **kwargs):
+    def interpolate(self, *new_array_grids, **kwargs):
         """
         Interpolate TimeTableCoordinate to new array index grids.
 
@@ -794,10 +864,18 @@ def interpolate(self, new_array_grids, **kwargs):
         """
         if self.is_scalar():
             raise ValueError("Cannot interpolate a scalar TimeTableCoordinate.")
-        # Build pixel grids for current TimeTableCoord.
-        old_array_grids = np.arange(len(self.table))
+        if len(new_array_grids) != max(self.table.ndim, 1):
+            raise ValueError(
+                f"A new array grid must be given for each array axis, i.e. {self.table.ndim}")
         # Interpolate using MJD format and convert back to a Time object.
-        new_table = np.interp(new_array_grids, old_array_grids, self.table.mjd, **kwargs)
+        if self.table.ndim == 1:
+            # Build pixel grids for current TimeTableCoord.
+            old_array_grids = np.arange(len(self.table))
+            new_table = np.interp(new_array_grids[0], old_array_grids, self.table.mjd, **kwargs)
+        else:
+            old_array_grids = tuple(np.arange(d) for d in self.table.shape)
+            new_table = scipy.interpolate.interpn(
+                old_array_grids, self.table.mjd, np.stack(new_array_grids, axis=-1), **kwargs)
         new_table = Time(new_table, scale=self.table.scale, format="mjd")
         new_table.format = self.table.format
         # Rebuild new TimeTableCoord and return.
diff --git a/ndcube/extra_coords/tests/test_extra_coords.py b/ndcube/extra_coords/tests/test_extra_coords.py
index e48190ace..266d56415 100644
--- a/ndcube/extra_coords/tests/test_extra_coords.py
+++ b/ndcube/extra_coords/tests/test_extra_coords.py
@@ -560,3 +560,35 @@ def test_length1_extra_coord(wave_lut):
     sec = ec[item]
     assert (sec.wcs.pixel_to_world(0) == wave_lut[item]).all()
     assert (sec.wcs.world_to_pixel(wave_lut[item])[0] == [0]).all()
+
+
+def test_2d_time_extra_coord_through_cube(wcs_3d_lt_ln_l):
+    cube = NDCube(np.zeros((3, 4, 5)), wcs=wcs_3d_lt_ln_l)
+    times = Time("2020-01-01T00:00:00") + np.arange(12).reshape(3, 4) * u.s
+    cube.extra_coords.add("time", (0, 1), times, physical_types="time")
+
+    (world_times,) = cube.axis_world_coords("time", wcs=cube.extra_coords)
+    assert world_times.shape == (3, 4)
+    assert (world_times == times).all()
+
+    # Slicing with ranges keeps the table 2-D.
+    sub = cube[1:3, 0:2]
+    (sub_times,) = sub.axis_world_coords("time", wcs=sub.extra_coords)
+    assert sub_times.shape == (2, 2)
+    assert (sub_times == times[1:3, 0:2]).all()
+
+    # Integer slicing drops the corresponding table dimension.
+    row = cube[1]
+    (row_times,) = row.axis_world_coords("time", wcs=row.extra_coords)
+    assert row_times.shape == (4,)
+    assert (row_times == times[1]).all()
+
+    column = cube[:, 2]
+    (column_times,) = column.axis_world_coords("time", wcs=column.extra_coords)
+    assert column_times.shape == (3,)
+    assert (column_times == times[:, 2]).all()
+
+    # Slicing away both table dimensions drops the coordinate.
+    point = cube[1, 2]
+    assert point.extra_coords.is_empty
+    assert len(point.extra_coords._dropped_tables) == 1
diff --git a/ndcube/extra_coords/tests/test_lookup_table_coord.py b/ndcube/extra_coords/tests/test_lookup_table_coord.py
index 35ae13886..2a29ec243 100644
--- a/ndcube/extra_coords/tests/test_lookup_table_coord.py
+++ b/ndcube/extra_coords/tests/test_lookup_table_coord.py
@@ -722,3 +722,96 @@ def assert_lutc_ancilliary_data_same(lutc1, lutc2):
     assert lutc1.names == lutc2.names
     assert lutc1.physical_types == lutc2.physical_types
     assert lutc1._dropped_world_dimensions == lutc2._dropped_world_dimensions
+
+
+@pytest.fixture
+def timetable_2d():
+    times = Time("2020-01-01T00:00:00") + np.arange(12).reshape(3, 4) * u.min
+    return TimeTableCoordinate(times, names="time", physical_types="time")
+
+
+@pytest.fixture
+def quantitytable_2d():
+    return QuantityTableCoordinate(np.arange(12).reshape(3, 4) * u.km,
+                                   names="distance", physical_types="pos.distance")
+
+
+def test_2d_time_table(timetable_2d):
+    assert timetable_2d.n_inputs == 2
+    assert not timetable_2d.is_scalar()
+
+    twcs = timetable_2d.wcs
+    assert twcs.pixel_n_dim == 2
+    assert twcs.world_n_dim == 1
+    # Model inputs are in pixel order, i.e. reversed array order.
+    assert twcs.pixel_to_world(1, 2) == timetable_2d.table[2, 1]
+
+
+def test_2d_time_table_slicing(timetable_2d):
+    sub = timetable_2d[1:3, 0:2]
+    assert isinstance(sub, TimeTableCoordinate)
+    assert sub.table.shape == (2, 2)
+    assert sub.n_inputs == 2
+    assert (sub.table == timetable_2d.table[1:3, 0:2]).all()
+    # Reference time must be preserved through slicing.
+    assert sub.reference_time == timetable_2d.reference_time
+
+    row = timetable_2d[1, :]
+    assert row.table.shape == (4,)
+    assert row.n_inputs == 1
+
+    scalar = timetable_2d[1, 2]
+    assert scalar.is_scalar()
+
+    with pytest.raises(ValueError, match="Can not slice with incorrect length"):
+        timetable_2d[1]
+
+
+def test_2d_time_table_interpolate(timetable_2d):
+    new_grid0 = np.array([0.5, 1.5])
+    new_grid1 = np.array([1.0, 2.0])
+    new_coord = timetable_2d.interpolate(new_grid0, new_grid1)
+    assert isinstance(new_coord, TimeTableCoordinate)
+    expected = timetable_2d.table.mjd[0:2, 1:3].diagonal() + 0.5 * (1 * u.min).to_value(u.day) * 4
+    np.testing.assert_allclose(new_coord.table.mjd, expected)
+
+
+def test_2d_quantity_table(quantitytable_2d):
+    assert quantitytable_2d.n_inputs == 2
+    assert quantitytable_2d.ndim == 2
+    assert quantitytable_2d.shape == (3, 4)
+    assert not quantitytable_2d.is_scalar()
+
+    qwcs = quantitytable_2d.wcs
+    assert qwcs.pixel_n_dim == 2
+    assert qwcs.world_n_dim == 1
+    # Model inputs are in pixel order, i.e. reversed array order.
+    assert qwcs.pixel_to_world(1, 2) == quantitytable_2d.table[0][2, 1]
+
+
+def test_2d_quantity_table_slicing(quantitytable_2d):
+    sub = quantitytable_2d[0:2, 1:3]
+    assert isinstance(sub, QuantityTableCoordinate)
+    assert sub.table[0].shape == (2, 2)
+    assert sub.n_inputs == 2
+    assert (sub.table[0] == quantitytable_2d.table[0][0:2, 1:3]).all()
+
+    row = quantitytable_2d[0, :]
+    assert row.table[0].shape == (4,)
+    assert row.n_inputs == 1
+
+    scalar = quantitytable_2d[0, 1]
+    assert scalar.is_scalar()
+
+
+def test_2d_quantity_table_interpolate(quantitytable_2d):
+    new_grid0 = np.array([0.0, 1.0])
+    new_grid1 = np.array([1.0, 2.0])
+    new_coord = quantitytable_2d.interpolate(new_grid0, new_grid1)
+    assert isinstance(new_coord, QuantityTableCoordinate)
+    np.testing.assert_allclose(new_coord.table[0].to_value(u.km), [1.0, 6.0])
+
+
+def test_multiple_nd_tables_rejected():
+    with pytest.raises(ValueError, match="Multiple tables can only be provided if they are all 1-D"):
+        QuantityTableCoordinate(np.ones((2, 2)) * u.km, np.ones((2, 2)) * u.km)

From cae3f77c8b825f2f3340f2a5a0f68c423008faaa Mon Sep 17 00:00:00 2001
From: Nabil Freij <nabil.freij@gmail.com>
Date: Fri, 12 Jun 2026 13:39:53 -0700
Subject: [PATCH 2/2] Clean up

---
 ndcube/extra_coords/extra_coords.py           |  2 +-
 ndcube/extra_coords/table_coord.py            | 21 +++---
 .../extra_coords/tests/test_extra_coords.py   |  3 +-
 .../tests/test_lookup_table_coord.py          | 70 -------------------
 ndcube/tests/test_ndcube_axis_world_coords.py |  4 +-
 5 files changed, 14 insertions(+), 86 deletions(-)

diff --git a/ndcube/extra_coords/extra_coords.py b/ndcube/extra_coords/extra_coords.py
index 5672f7fd3..eee15980e 100644
--- a/ndcube/extra_coords/extra_coords.py
+++ b/ndcube/extra_coords/extra_coords.py
@@ -265,7 +265,7 @@ def mapping(self):
         for lut_axis, lut in self._lookup_tables:
             ids = [lut_axis] if isinstance(lut_axis, Integral) else list(lut_axis)
             pixel_ids = list(converter(np.array(ids)))
-            if getattr(lut, "_model_inputs_are_pixel_ordered", False):
+            if lut._model_inputs_are_pixel_ordered:
                 # Single N-D tables expose their model inputs in pixel order,
                 # i.e. reversed with respect to the array-ordered axes given
                 # to `add`.
diff --git a/ndcube/extra_coords/table_coord.py b/ndcube/extra_coords/table_coord.py
index a9b0e6f74..8179e329b 100644
--- a/ndcube/extra_coords/table_coord.py
+++ b/ndcube/extra_coords/table_coord.py
@@ -294,6 +294,13 @@ def _model_inputs_are_pixel_ordered(self):
         """
         return False
 
+    @staticmethod
+    def _reorder_inputs_to_pixel(model):
+        """
+        Reverse a model's inputs from array order to pixel order.
+        """
+        return models.Mapping(tuple(range(model.n_inputs))[::-1]) | model
+
     @property
     def wcs(self):
         """
@@ -430,9 +437,7 @@ def __getitem__(self, item):
 
     @property
     def n_inputs(self):
-        if self._single_nd_table:
-            return self.table[0].ndim
-        return len(self.table)
+        return self.ndim
 
     def is_scalar(self):
         return all(t.shape == () for t in self.table)
@@ -456,9 +461,7 @@ def model(self):
         """
         model = _model_from_quantity(self.table, True)
         if self._single_nd_table:
-            # Expose the inputs of the N-D table in pixel order, i.e. reversed
-            # with respect to the table's (array-ordered) dimensions.
-            model = models.Mapping(tuple(range(model.n_inputs))[::-1]) | model
+            model = self._reorder_inputs_to_pixel(model)
         return model
 
     @property
@@ -838,9 +841,7 @@ def model(self):
 
         model = _model_from_quantity((deltas,), mesh=False)
         if deltas.ndim > 1:
-            # Expose the inputs of the N-D table in pixel order, i.e. reversed
-            # with respect to the table's (array-ordered) dimensions.
-            model = models.Mapping(tuple(range(model.n_inputs))[::-1]) | model
+            model = self._reorder_inputs_to_pixel(model)
         return model
 
     def interpolate(self, *new_array_grids, **kwargs):
@@ -864,7 +865,7 @@ def interpolate(self, *new_array_grids, **kwargs):
         """
         if self.is_scalar():
             raise ValueError("Cannot interpolate a scalar TimeTableCoordinate.")
-        if len(new_array_grids) != max(self.table.ndim, 1):
+        if len(new_array_grids) != self.table.ndim:
             raise ValueError(
                 f"A new array grid must be given for each array axis, i.e. {self.table.ndim}")
         # Interpolate using MJD format and convert back to a Time object.
diff --git a/ndcube/extra_coords/tests/test_extra_coords.py b/ndcube/extra_coords/tests/test_extra_coords.py
index 266d56415..47c18f223 100644
--- a/ndcube/extra_coords/tests/test_extra_coords.py
+++ b/ndcube/extra_coords/tests/test_extra_coords.py
@@ -283,11 +283,10 @@ def test_extra_coords_index(skycoord_2d_lut, time_lut):
     assert sub_ec.wcs.world_axis_names == ("exposure_time",)
 
 
-@pytest.mark.xfail(reason=">1D Tables not supported")
 def test_extra_coords_2d_quantity(quantity_2d_lut):
     ec = ExtraCoords()
     ec.add("velocity", (0, 1), quantity_2d_lut)
-    assert ec.wcs.pixel_to_world(0, 0)
+    assert u.allclose(ec.wcs.pixel_to_world(1, 2), quantity_2d_lut[2, 1])
 
 
 # Extra Coords with NDCube
diff --git a/ndcube/extra_coords/tests/test_lookup_table_coord.py b/ndcube/extra_coords/tests/test_lookup_table_coord.py
index 2a29ec243..8aaf80a0b 100644
--- a/ndcube/extra_coords/tests/test_lookup_table_coord.py
+++ b/ndcube/extra_coords/tests/test_lookup_table_coord.py
@@ -92,22 +92,6 @@ def test_3d_distance(lut_3d_distance_mesh):
     assert u.allclose(ltc.wcs.world_to_pixel(0*u.km, 10*u.km, 20*u.km), (0, 0, 0))
 
 
-@pytest.mark.xfail(reason=">1D Tables not supported")
-def test_2d_nout_1_no_mesh(lut_2d_distance_no_mesh):
-    ltc = lut_2d_distance_no_mesh
-    assert ltc.wcs.world_n_dim == 2
-    assert ltc.wcs.pixel_n_dim == 2
-
-    assert ltc.model.n_inputs == 2
-    assert ltc.model.n_outputs == 2
-
-    assert u.allclose(ltc.wcs.pixel_to_world(0*u.pix, 0*u.pix),
-                      (0, 9)*u.km)
-
-    # TODO: this model is not invertable
-    assert u.allclose(ltc.wcs.world_to_pixel(0*u.km, 9*u.km), (0, 0))
-
-
 def test_1d_skycoord_no_mesh(lut_1d_skycoord_no_mesh):
     ltc = lut_1d_skycoord_no_mesh
 
@@ -207,7 +191,6 @@ def test_join_3d(lut_2d_skycoord_mesh, lut_1d_wave):
     assert u.allclose(ltc.wcs.world_to_pixel(*world), (0, 0, 0))
 
 
-@pytest.mark.xfail(reason=">1D Tables not supported")
 def test_2d_quantity():
     shape = (3, 3)
     data = np.arange(np.prod(shape)).reshape(shape) * u.m / u.s
@@ -336,14 +319,6 @@ def test_3d_distance_slice(lut_3d_distance_mesh):
     assert len(sub_ltc.table[2]) == 7
 
 
-@pytest.mark.xfail(reason=">1D Tables not supported")
-def test_2d_nout_1_no_mesh_slice(lut_2d_distance_no_mesh):
-    ltc = lut_2d_distance_no_mesh
-    sub_ltc = ltc[0:2, 0:2]
-    assert sub_ltc.table[0].shape == (2, 2)
-    assert sub_ltc.table[1].shape == (2, 2)
-
-
 def test_1d_skycoord_no_mesh_slice(lut_1d_skycoord_no_mesh):
     sub_ltc = lut_1d_skycoord_no_mesh[0:4]
     assert sub_ltc.table.shape == (4, )
@@ -466,33 +441,6 @@ def test_mtc_dropped_table_skycoord_join(lut_1d_time, lut_2d_skycoord_mesh):
     assert dwd["value"] == [0, 0]
 
 
-@pytest.mark.xfail(reason=">1D Tables not supported")
-def test_mtc_dropped_quantity_table(lut_1d_time, lut_2d_distance_no_mesh):
-    mtc = MultipleTableCoordinate(lut_1d_time, lut_2d_distance_no_mesh)
-    sub = mtc[:, 0, 0]
-
-    assert len(sub._table_coords) == 1
-    assert len(sub._dropped_coords) == 1
-
-    pytest.importorskip("gwcs", minversion="0.17")
-
-    dwd = sub.dropped_world_dimensions
-    assert isinstance(dwd, dict)
-    wao_classes = dwd.pop("world_axis_object_classes")
-    assert all(isinstance(value, list) for value in dwd.values())
-    assert dwd
-    assert all(len(value) == 2 for value in dwd.values())
-
-    assert dwd["world_axis_names"] == ["", ""]
-    assert all(isinstance(u, str) for u in dwd["world_axis_units"])
-    assert dwd["world_axis_units"] == ["km", "km"]
-    assert dwd["world_axis_physical_types"] == ["custom:SPATIAL", "custom:SPATIAL"]
-    assert dwd["world_axis_object_components"] == [("SPATIAL", 0, "value"), ("SPATIAL1", 0, "value")]
-    assert wao_classes["SPATIAL"][0] is u.Quantity
-    assert wao_classes["SPATIAL1"][0] is u.Quantity
-    assert dwd["value"] == [0*u.km, 9*u.km]
-
-
 def test_mtc_dropped_quantity_inside_table(lut_3d_distance_mesh):
     sub = lut_3d_distance_mesh[:, 0, :]
 
@@ -519,24 +467,6 @@ def test_mtc_dropped_quantity_inside_table(lut_3d_distance_mesh):
     assert all(len(value) == 2 for value in dwd.values())
 
 
-@pytest.mark.xfail(reason=">1D Tables not supported")
-def test_mtc_dropped_quantity_inside_table_no_mesh(lut_2d_distance_no_mesh):
-    """
-    When not meshing, we don't drop a coord, as the coordinate for the sliced
-    out axis can still vary along the remaining coordinate.
-    """
-    sub = lut_2d_distance_no_mesh[:, 0]
-
-    assert len(sub.table) == 2
-
-    pytest.importorskip("gwcs", minversion="0.17")
-
-    dwd = sub.dropped_world_dimensions
-    assert isinstance(dwd, dict)
-    dwd.pop("world_axis_object_classes")
-    assert not dwd
-
-
 def test_mtc_dropped_quantity_join_drop_table(lut_1d_time, lut_3d_distance_mesh):
     mtc = MultipleTableCoordinate(lut_1d_time, lut_3d_distance_mesh)
     sub = mtc[:, 0, :, :]
diff --git a/ndcube/tests/test_ndcube_axis_world_coords.py b/ndcube/tests/test_ndcube_axis_world_coords.py
index 2d20a5cd2..17758658f 100644
--- a/ndcube/tests/test_ndcube_axis_world_coords.py
+++ b/ndcube/tests/test_ndcube_axis_world_coords.py
@@ -64,14 +64,12 @@ def test_axis_world_coords_empty_ec(ndcube_3d_l_ln_lt_ectime):
     assert awc == ()
 
 
-@pytest.mark.xfail(reason=">1D Tables not supported")
 def test_axis_world_coords_complex_ec(ndcube_4d_ln_lt_l_t):
     cube = ndcube_4d_ln_lt_l_t
     ec_shape = cube.data.shape[1:3]
     data = np.arange(np.prod(ec_shape)).reshape(ec_shape) * u.m / u.s
 
-    # The lookup table has to be in world order so transpose it.
-    cube.extra_coords.add('velocity', (2, 1), data.T)
+    cube.extra_coords.add('velocity', (1, 2), data)
 
     coords = cube.axis_world_coords(wcs=cube.extra_coords)
     assert len(coords) == 1