diff --git a/changelog/9999.feature.1.rst b/changelog/9999.feature.1.rst
new file mode 100644
index 000000000..e76a20117
--- /dev/null
+++ b/changelog/9999.feature.1.rst
@@ -0,0 +1 @@
+Added the ``NDCube._extra_attrs_to_copy`` extension point. Subclasses can declare instance attributes that are automatically propagated to cubes derived through arithmetic operations (``_new_instance``) and through ``to_nddata`` when the target type is a subclass carrying the same attributes.
diff --git a/changelog/9999.feature.2.rst b/changelog/9999.feature.2.rst
new file mode 100644
index 000000000..5c15dffef
--- /dev/null
+++ b/changelog/9999.feature.2.rst
@@ -0,0 +1 @@
+The error raised by `~ndcube.NDCube.crop` when a point has the wrong number of components now lists the expected world objects and their order.
diff --git a/changelog/9999.feature.rst b/changelog/9999.feature.rst
new file mode 100644
index 000000000..b93850f0c
--- /dev/null
+++ b/changelog/9999.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)
diff --git a/ndcube/ndcube.py b/ndcube/ndcube.py
index a9bcbb0ed..d2d4c8b42 100644
--- a/ndcube/ndcube.py
+++ b/ndcube/ndcube.py
@@ -378,6 +378,15 @@ class NDCubeBase(NDCubeABC, astropy.nddata.NDData, NDCubeSlicingMixin):
     _extra_coords = NDCubeLinkedDescriptor(ExtraCoords)
     _global_coords = NDCubeLinkedDescriptor(GlobalCoords)
 
+    # Names of additional instance attributes which subclasses want propagated
+    # (by reference) to instances derived through `_new_instance` (e.g.
+    # arithmetic operations) and `to_nddata` when the target type carries the
+    # same attributes. Subclasses should override this with a tuple of
+    # attribute names. Note these attributes are not automatically modified
+    # when a cube is sliced; subclasses with shape-dependent attributes must
+    # handle slicing themselves.
+    _extra_attrs_to_copy = ()
+
     def __init__(self, data, wcs=None, uncertainty=None, mask=None, meta=None,
                  unit=None, copy=False, psf=None, *, extra_coords=None, global_coords=None, **kwargs):
 
@@ -642,8 +651,13 @@ def _get_crop_item(self, *points, wcs=None, keepdims=False):
         classes = [wcs.world_axis_object_classes[c][0] for c in comp]
         for i, point in enumerate(points):
             if len(point) != len(comp):
+                component_names = ", ".join(
+                    f"{name} ({cls.__name__})" for name, cls in zip(comp, classes))
                 raise ValueError(f"{len(point)} components in point {i} do not match "
-                                 f"WCS with {len(comp)} components.")
+                                 f"WCS with {len(comp)} components. Each point must "
+                                 "have one entry per world object (use None for a "
+                                 "component that should not be cropped), in order: "
+                                 f"{component_names}.")
             for j, value in enumerate(point):
                 if not (value is None or isinstance(value, classes[j])):
                     raise TypeError(f"{type(value)} of component {j} in point {i} is "
@@ -963,6 +977,9 @@ def _new_instance(self, **kwargs):
             new_cube._extra_coords = deepcopy(self.extra_coords)
         if self.global_coords is not None:
             new_cube._global_coords = deepcopy(self.global_coords)
+        for attr in self._extra_attrs_to_copy:
+            if hasattr(self, attr):
+                setattr(new_cube, attr, getattr(self, attr))
         return new_cube
 
     def __neg__(self):
@@ -1635,6 +1652,12 @@ def to_nddata(self,
           array([[1., 1., 1.],
                  [1., 1., 1.]])
         """
+        # If the target type carries the same subclass-specific attributes as
+        # this cube, copy them by default unless explicitly overridden.
+        if isinstance(nddata_type, type) and issubclass(nddata_type, type(self)):
+            for attr in self._extra_attrs_to_copy:
+                if hasattr(self, attr):
+                    kwargs.setdefault(attr, "copy")
         # Put all NDData kwargs in a dict
         user_kwargs = {"data": data,
                        "wcs": wcs,
diff --git a/ndcube/tests/test_ndcube.py b/ndcube/tests/test_ndcube.py
index 7719b842c..1592d216f 100644
--- a/ndcube/tests/test_ndcube.py
+++ b/ndcube/tests/test_ndcube.py
@@ -293,3 +293,47 @@ def __init__(self, data, *, spam=None, **kwargs):
     assert new_ndd.spam == "Eggs"
     assert new_ndd.data is ndc.data
     assert new_ndd.wcs is ndc.wcs
+
+
+class SidecarCube(NDCube):
+    """NDCube subclass declaring extra attributes to propagate to derived cubes."""
+
+    _extra_attrs_to_copy = ("observer", "calibration_level")
+
+    def __init__(self, *args, **kwargs):
+        self.observer = kwargs.pop("observer", None)
+        self.calibration_level = kwargs.pop("calibration_level", 0)
+        super().__init__(*args, **kwargs)
+
+
+@pytest.fixture
+def sidecar_cube(wcs_3d_lt_ln_l):
+    cube = SidecarCube(np.ones((2, 3, 4)), wcs=wcs_3d_lt_ln_l)
+    cube.observer = "earth"
+    cube.calibration_level = 2
+    return cube
+
+
+def test_extra_attrs_to_copy_propagate_through_arithmetic(sidecar_cube):
+    doubled = sidecar_cube * 2
+    assert type(doubled) is SidecarCube
+    assert doubled.observer == "earth"
+    assert doubled.calibration_level == 2
+
+    negated = -sidecar_cube
+    assert negated.observer == "earth"
+    assert negated.calibration_level == 2
+
+
+def test_extra_attrs_to_copy_propagate_through_to_nddata(sidecar_cube):
+    copied = sidecar_cube.to_nddata(nddata_type=SidecarCube)
+    assert copied.observer == "earth"
+    assert copied.calibration_level == 2
+
+    # Explicit kwargs override the automatic copy.
+    overridden = sidecar_cube.to_nddata(nddata_type=SidecarCube, observer="sdo")
+    assert overridden.observer == "sdo"
+
+    # Types which do not carry the attributes are unaffected.
+    plain = sidecar_cube.to_nddata(nddata_type=astropy.nddata.NDData)
+    assert not hasattr(plain, "observer")
diff --git a/ndcube/tests/test_ndcube_slice_and_crop.py b/ndcube/tests/test_ndcube_slice_and_crop.py
index 56b7cd4a4..81ca3d745 100644
--- a/ndcube/tests/test_ndcube_slice_and_crop.py
+++ b/ndcube/tests/test_ndcube_slice_and_crop.py
@@ -617,6 +617,13 @@ def test_crop_all_points_beyond_cube_extent_error(points):
         cube.crop(*points, keepdims=True)
 
 
+def test_crop_length_error_names_world_components(ndcube_4d_ln_lt_l_t):
+    cube = ndcube_4d_ln_lt_l_t
+    interval0 = cube.wcs.array_index_to_world([1, 2], [0, 1], [0, 1], [0, 2])[0]
+    with pytest.raises(ValueError, match=r"one entry per world object"):
+        cube.crop([interval0[0], None], [interval0[-1], None])
+
+
 def test_crop_by_values_quantity_table_coordinate():
     # Regression: QuantityTableCoordinate-based WCS raised
     # "High Level objects are not supported with the native API" because