diff --git a/conftest.py b/conftest.py
index 2206968d..eefc2a68 100644
--- a/conftest.py
+++ b/conftest.py
@@ -1,11 +1,33 @@
 import pytest
 
+from sasdata.metadata import Instrument, Metadata, Source
+from sasdata.quantities.quantity import Quantity
+from sasdata.quantities.units import angstroms
+
 
 def pytest_addoption(parser):
-    parser.addoption(
-        "--show_plots", action="store_true", default=False, help="Display diagnostic plots during tests"
-    )
+    parser.addoption("--show_plots", action="store_true", default=False, help="Display diagnostic plots during tests")
+
 
 @pytest.fixture
 def show_plots(request):
     return request.config.getoption("--show_plots")
+
+
+@pytest.fixture
+def basic_metadata():
+    """A minimal metadata object."""
+    wavelength = Quantity(1.0, angstroms)
+    source = Source(
+        radiation=None,
+        beam_shape=None,
+        beam_size=None,
+        wavelength=wavelength,
+        wavelength_max=None,
+        wavelength_min=None,
+        wavelength_spread=None,
+    )
+    instrument = Instrument(collimations=[], source=source, detector=[])
+    metadata = Metadata(title=None, run=[], definition=None, process=[], sample=None, instrument=instrument, raw=None)
+
+    return metadata
diff --git a/sasdata/abscissa.py b/sasdata/abscissa.py
new file mode 100644
index 00000000..c0bb8519
--- /dev/null
+++ b/sasdata/abscissa.py
@@ -0,0 +1,121 @@
+from abc import ABC, abstractmethod
+
+import numpy as np
+from numpy._typing import ArrayLike
+
+from sasdata.exceptions import InterpretationError
+from sasdata.quantities.quantity import Quantity
+from sasdata.util import is_increasing
+
+
+class Abscissa(ABC):
+    def __init__(self, axes: list[Quantity]):
+        self._axes = axes
+        self._dimensionality = len(axes)
+
+    @property
+    def dimensionality(self) -> int:
+        """Dimensionality of this data"""
+        return self._dimensionality
+
+    @property
+    @abstractmethod
+    def is_grid(self) -> bool:
+        """Are these coordinates using a grid representation
+        (as opposed to a general list representation)
+
+        is_grid = True: implies that the corresponding ordinate is n-dimensional tensor
+        is_grid = False: implies that the corresponding ordinate is a 1D list
+
+        If the data is one dimensional, is_grid=True
+        """
+
+    @property
+    def axes(self) -> list[Quantity]:
+        """Axes of the data:
+
+        If it's an (n1-by-n2-by-n3...) grid (is_grid=True): give the values for each axis, returning a list like
+          [Quantity(length n1), Quantity(length n2), Quantity(length n3) ... ]
+
+        If it is not grid data (is_grid=False), but n points on a general mesh, give one array for each dimension
+          [Quantity(length n), Quantity(length n), Quantity(length n) ... ]
+        """
+        return self._axes
+
+    @staticmethod
+    def _determine_error_message(axis_arrays: list[np.ndarray], ordinate_shape: tuple):
+        """Error message for the `.determine` function"""
+
+        shape_string = ", ".join([str(axis.shape) for axis in axis_arrays])
+
+        return f"Cannot interpret array shapes axis: [{shape_string}], ordinate: {ordinate_shape}"
+
+    @staticmethod
+    def determine(axis_data: list[Quantity[ArrayLike]], ordinate_data: Quantity[ArrayLike]) -> "Abscissa":
+        """Get an Abscissa object that fits the combination of axes and data"""
+
+        # Different posibilites:
+        #   1: axes_data[i].shape == axes_data[j].shape == ordinate_data.shape
+        #    1a: axis_data[i] is 1D =>
+        #      1a-i:  len(axes_data) == 1 => Grid type or Scatter type depending on sortedness
+        #      1a-ii: len(axes_data) != 1 => Scatter type
+        #    1b: axis_data[i] dimensionality matches len(axis_data) => Meshgrid type
+        #    1c: other => Error
+        #   2: (len(axes_data[0]), len(axes_data[1])... ) == ordinate_data.shape => Grid type
+        #   3: None of the above => Error
+
+        ordinate_shape = np.array(ordinate_data.value).shape
+        axis_arrays = [np.array(axis.value) for axis in axis_data]
+
+        # 1:
+        if all([axis.shape == ordinate_shape for axis in axis_arrays]):
+            # 1a:
+            if all([len(axis.shape) == 1 for axis in axis_arrays]):
+                # 1a-i:
+                if len(axis_arrays) == 1:
+                    # Is it sorted
+                    if is_increasing(axis_arrays[0]):
+                        return GridAbscissa(axis_data)
+                    else:
+                        return ScatterAbscissa(axis_data)
+                # 1a-ii
+                else:
+                    return ScatterAbscissa(axis_data)
+            # 1b
+            elif all([len(axis.shape) == len(axis_arrays) for axis in axis_arrays]):
+                return MeshgridAbscissa(axis_data)
+
+            else:
+                raise InterpretationError(Abscissa._determine_error_message(axis_arrays, ordinate_shape))
+
+        elif (
+            all([len(axis.shape) == 1 for axis in axis_arrays])
+            and tuple([axis.shape[0] for axis in axis_arrays]) == ordinate_shape
+        ):
+            # Require that they are sorted
+            if all([is_increasing(axis) for axis in axis_arrays]):
+                return GridAbscissa(axis_data)
+
+            else:
+                raise InterpretationError("Grid axes are not sorted")
+
+        else:
+            raise InterpretationError(Abscissa._determine_error_message(axis_arrays, ordinate_shape))
+
+
+class GridAbscissa(Abscissa):
+    @property
+    def is_grid(self):
+        return True
+
+
+class MeshgridAbscissa(Abscissa):
+    @property
+    def is_grid(self):
+        return True
+
+
+class ScatterAbscissa(Abscissa):
+    @property
+    def is_grid(self):
+        return False
diff --git a/sasdata/data.py b/sasdata/data.py
index 1a7101f9..a9cd9648 100644
--- a/sasdata/data.py
+++ b/sasdata/data.py
@@ -7,12 +7,66 @@
 from h5py._hl.group import Group as HDF5Group
 
 from sasdata import dataset_types
+from sasdata.abscissa import Abscissa
 from sasdata.dataset_types import DatasetType
-from sasdata.metadata import Metadata, MetadataEncoder
+from sasdata.metadata import DerivedMetadata, Metadata, MetadataEncoder
 from sasdata.quantities.quantity import Quantity
 
 
 class SasData:
+    """General object containing data in the SasView ecosystem"""
+
+    def __init__(
+        self,
+        name: str,
+        ordinate: Quantity,
+        abscissae: Abscissa,
+        mask: Quantity,
+        dependents: list["SasData"],
+        metadata: Metadata,
+    ):
+        self.name = name
+        self._ordinate = ordinate
+        self._abscissae = abscissae
+        self._mask = mask
+        self.dependents = dependents
+        self.metadata = metadata
+
+    @property
+    def ordinate(self) -> Quantity:
+        return self._ordinate
+
+    @property
+    def abscissae(self) -> Abscissa:
+        return self._abscissae
+
+    @property
+    def mask(self) -> Quantity:
+        return self._mask
+
+    def scatter_data(self):
+        """Return data in the coordinate/value form [(x1, x2, x3, y)...]"""
+
+
+class SasDerivedMeasurement(SasData):
+    """General object sas measurement that has not come directly from a file,
+    for example, the difference between two datasets"""
+
+    def __init__(
+        self,
+        name: str,
+        ordinate: Quantity,
+        abscissae: Abscissa,
+        mask: Quantity,
+        dependents: list["SasData"],
+        metadata: DerivedMetadata,
+    ):
+        super().__init__(
+            name=name, ordinate=ordinate, abscissae=abscissae, mask=mask, dependents=dependents, metadata=metadata
+        )
+
+
+class SasMeasurement(SasData):
     def __init__(
         self,
         name: str,
@@ -34,16 +88,14 @@ def __init__(
         self.dataset_type: DatasetType = dataset_type
 
         # Components that need to be organised after creation
-        self.mask = None  # TODO: fill out
+        self._mask = None  # TODO: fill out
         self.model_requirements = None  # TODO: fill out
 
     # TODO: Handle the other data types.
     @property
     def ordinate(self) -> Quantity:
         match self.dataset_type:
-            case (dataset_types.one_dim |
-                  dataset_types.two_dim |
-                  dataset_types.angle_dim):
+            case dataset_types.one_dim | dataset_types.two_dim | dataset_types.three_dim | dataset_types.angle_dim:
                 return self._data_contents["I"]
             case dataset_types.sesans:
                 return self._data_contents["Depolarisation"]
@@ -51,35 +103,22 @@ def ordinate(self) -> Quantity:
                 return None
 
     @property
-    def abscissae(self) -> Quantity:
+    def abscissae(self) -> Abscissa:
         match self.dataset_type:
             case dataset_types.one_dim:
-                return self._data_contents["Q"]
+                return Abscissa.determine([self._data_contents["Q"]], self.ordinate)
             case dataset_types.two_dim:
-                # Type hinting is a bit lacking. Assume each part of the zip is a scalar value.
-                data_contents = np.array(
-                    list(
-                        zip(
-                            self._data_contents["Qx"].value,
-                            self._data_contents["Qy"].value,
-                        )
-                    )
-                )
-                # Use this value to extract units etc. Assume they will be the same for Qy.
-                reference_data_content = self._data_contents["Qx"]
-                # TODO: If this is a derived quantity then we are going to lose that
-                # information.
-                #
-                # TODO: Won't work when there's errors involved. On reflection, we
-                # probably want to avoid creating a new Quantity but at the moment I
-                # can't see a way around it.
-                return Quantity(data_contents, reference_data_content.units, name=self._data_contents["Qx"].name, id_header=self._data_contents["Qx"]._id_header)
+                return Abscissa.determine([self._data_contents["Qx"], self._data_contents["Qy"]], self.ordinate)
             case dataset_types.angle_dim:
-                return self._data_contents["Phi"]
+                return Abscissa.determine([self._data_contents["Phi"]], self.ordinate)
+            case dataset_types.three_dim:
+                return Abscissa.determine(
+                    [self._data_contents["Qx"], self._data_contents["Qy"], self._data_contents["Qz"]], self.ordinate
+                )
             case dataset_types.sesans:
-                return self._data_contents["SpinEchoLength"]
+                return Abscissa.determine([self._data_contents["SpinEchoLength"]], self.ordinate)
             case _:
-                None
+                return None
 
     def __getitem__(self, item: str):
         return self._data_contents[item]
@@ -98,7 +137,7 @@ def summary(self, indent="  "):
 
     @staticmethod
     def from_json(obj):
-        return SasData(
+        return SasMeasurement(
             name=obj["name"],
             dataset_type=DatasetType(
                 name=obj["type"]["name"],
@@ -121,12 +160,11 @@ def _save_h5(self, sasentry: HDF5Group):
         for idx, (key, sasdata) in enumerate(self._data_contents.items()):
             sasdata.as_h5(group, key)
 
-
     @staticmethod
     def save_h5(data: dict[str, typing.Self], path: str | typing.BinaryIO):
         with h5py.File(path, "w") as f:
             for idx, (key, data) in enumerate(data.items()):
-                sasentry = f.create_group(f"sasentry{idx+1:02d}")
+                sasentry = f.create_group(f"sasentry{idx + 1:02d}")
                 if not key.startswith("sasentry"):
                     sasentry.attrs["sasview_key"] = key
                 data._save_h5(sasentry)
@@ -140,7 +178,7 @@ def deserialise(data: str) -> "SasData":
     def deserialise_json(json_data: dict) -> "SasData":
         name = json_data["name"]
         data_contents = {}
-        dataset_type = json_data["dataset_type"] # TODO: update when DatasetType is more finalized
+        dataset_type = json_data["dataset_type"]  # TODO: update when DatasetType is more finalized
         metadata = json_data["metadata"].deserialise_json()
         for quantity in json_data["data_contents"]:
             data_contents[quantity["label"]] = Quantity.deserialise_json(quantity)
@@ -159,11 +197,11 @@ def _serialise_json(self) -> dict[str, Any]:
         return {
             "name": self.name,
             "data_contents": data,
-            "dataset_type": None, # TODO: update when DatasetType is more finalized
+            "dataset_type": None,  # TODO: update when DatasetType is more finalized
             "verbose": self._verbose,
             "metadata": self.metadata.serialise_json(),
             "mask": {},
-            "model_requirements": {}
+            "model_requirements": {},
         }
 
 
@@ -210,16 +248,16 @@ def sasdata_reader2D_converter(data2d: SasData | None = None) -> SasData:
     qx_data = new_x.flatten()
     qy_data = new_y.flatten()
     err_data = np.sqrt(data2d._data_contents["I"].variance.value)
-    if not data2d._data_contents["I"].has_variance or np.any(err_data <= 0):
+    if not data2d._data_contents["I"].has_error or np.any(err_data <= 0):
         new_err_data = np.sqrt(np.abs(new_data))
     else:
         new_err_data = err_data.flatten()
     mask = np.ones(len(new_data), dtype=bool)
 
     data2d._data_contents["I"].value = new_data
-    data2d._data_contents["I"].variance.value = new_err_data ** 2
+    data2d._data_contents["I"].variance.value = new_err_data**2
     data2d._data_contents["Qx"].value = qx_data
     data2d._data_contents["Qy"].value = qy_data
-    data2d.mask = mask
+    data2d._mask = mask
 
     return data2d
diff --git a/sasdata/data_util/averaging.py b/sasdata/data_util/averaging.py
index 6a1e4bed..98ba9647 100644
--- a/sasdata/data_util/averaging.py
+++ b/sasdata/data_util/averaging.py
@@ -5,17 +5,17 @@
 import numpy as np
 import numpy.typing as npt
 
-from sasdata.data import SasData
+from sasdata.abscissa import Abscissa
+from sasdata.data import SasData, SasMeasurement
 from sasdata.data_util.binning import DirectionalAverage
 from sasdata.data_util.interval import IntervalType
 from sasdata.data_util.roi import CartesianROI, PolarROI
-from sasdata.dataset_types import angle_dim, one_dim
 from sasdata.quantities.constants import Pi, TwoPi
 from sasdata.quantities.quantity import Quantity
 from sasdata.quantities.units import radians
 
 
-def get_dq_data(data2d: SasData) -> npt.NDArray[np.floating]:
+def get_dq_data(data2d: SasMeasurement) -> npt.NDArray[np.floating]:
     """
     Get the dq for resolution averaging
     The pinholes and det. pix contribution present
@@ -25,13 +25,13 @@ def get_dq_data(data2d: SasData) -> npt.NDArray[np.floating]:
     Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average.
     """
 
-    q_data = np.sqrt(data2d._data_contents["Qx"].value**2 + data2d._data_contents["Qy"].value**2)
+    q_data = np.sqrt(data2d.abscissae.axes[0].value**2 + data2d.abscissae.axes[1].value**2)
 
     z_max = np.max(q_data)
     z_min = np.min(q_data)
 
-    dqx_data = np.sqrt(data2d._data_contents["Qx"].variance.value)
-    dqy_data = np.sqrt(data2d._data_contents["Qy"].variance.value)
+    dqx_data = np.sqrt(data2d.abscissae.axes[0].variance.value)
+    dqy_data = np.sqrt(data2d.abscissae.axes[1].variance.value)
 
     dqx_at_z_max = dqx_data[np.argmax(q_data)]
     dqx_at_z_min = dqx_data[np.argmin(q_data)]
@@ -78,11 +78,11 @@ def __init__(self, qx_range: tuple[float, float] = (0.0, 0.0), qy_range: tuple[f
         """
         super().__init__(qx_range=qx_range, qy_range=qy_range)
 
-    def __call__(self, data2d: SasData | None = None) -> tuple[float, float, float]:
+    def __call__(self, data2d: SasMeasurement | None = None) -> tuple[float, float, float]:
         """
         Coordinate data processing operations and return the results.
 
-        :param data2d: The SasData object for which the sum is calculated.
+        :param data2d: The SasMeasurement object for which the sum is calculated.
         """
         self.validate_and_assign_data(data2d)
         total_sum, error, count = self._sum()
@@ -126,11 +126,11 @@ def __init__(self, qx_range: tuple[float, float] = (0.0, 0.0), qy_range: tuple[f
         """
         super().__init__(qx_range=qx_range, qy_range=qy_range)
 
-    def __call__(self, data2d: SasData) -> tuple[float, float]:
+    def __call__(self, data2d: SasMeasurement) -> tuple[float, float]:
         """
         Coordinate data processing operations and return the results.
 
-        :param data2d: The SasData object for which the average is calculated.
+        :param data2d: The SasMeasurement object for which the average is calculated.
         """
         self.validate_and_assign_data(data2d)
         total_sum, error, count = super()._sum()
@@ -143,7 +143,7 @@ class SlabX(CartesianROI):
     Average I(Q_x, Q_y) along the y direction (within a ROI), giving I(Q_x).
 
     This class is initialised by specifying the boundaries of the ROI and is
-    called by supplying a SasData object. It returns a SasData object.
+    called by supplying a SasMeasurement object. It returns a SasData object.
     The averaging process can also be thought of as projecting 2D -> 1D.
 
     There also exists the option to "fold" the ROI, where Q data on opposite
@@ -173,11 +173,11 @@ def __init__(
         self.fold: bool = fold
         self.base: float | None = base
 
-    def __call__(self, data2d: SasData | None = None) -> SasData:
+    def __call__(self, data2d: SasMeasurement | None = None) -> SasData:
         """
         Compute the 1D average of 2D data, projecting along the Q_x axis.
 
-        :param data2d: The SasData object for which the average is computed.
+        :param data2d: The SasMeasurement object for which the average is computed.
         :return: SasData object for plotting.
         """
         self.validate_and_assign_data(data2d)
@@ -204,11 +204,17 @@ def __call__(self, data2d: SasData | None = None) -> SasData:
 
         qx_data, intensity, error = directional_average(data=self.data, err_data=self.err_data)
 
-        data_contents = {
-            "Q": Quantity(qx_data, data2d._data_contents["Qx"].units, None),
-            "I": Quantity(intensity, data2d.ordinate.units, error),
-        }
-        return SasData(f"{data2d.name}: Slab X Average", data_contents, one_dim, data2d.metadata)
+        ordinate = Quantity(intensity, data2d.ordinate.units, error)
+        abscissa = Abscissa.determine([Quantity(qx_data, data2d.abscissae.axes[0].units, None)], ordinate)
+
+        return SasData(
+            name=f"{data2d.name}: Slab X Average",
+            ordinate=ordinate,
+            abscissae=abscissa,
+            mask=data2d.mask,
+            dependents=[data2d],
+            metadata=data2d.metadata
+        )
 
 
 class SlabY(CartesianROI):
@@ -216,7 +222,7 @@ class SlabY(CartesianROI):
     Average I(Q_x, Q_y) along the x direction (within a ROI), giving I(Q_y).
 
     This class is initialised by specifying the boundaries of the ROI and is
-    called by supplying a SasData object. It returns a SasData object.
+    called by supplying a SasMeasurement object. It returns a SasData object.
     The averaging process can also be thought of as projecting 2D -> 1D.
 
     There also exists the option to "fold" the ROI, where Q data on opposite
@@ -247,11 +253,11 @@ def __init__(
         self.fold: bool = fold
         self.base: float | None = base
 
-    def __call__(self, data2d: SasData | None = None) -> SasData:
+    def __call__(self, data2d: SasMeasurement | None = None) -> SasData:
         """
         Compute the 1D average of 2D data, projecting along the Q_y axis.
 
-        :param data2d: The SasData object for which the average is computed.
+        :param data2d: The SasMeasurement object for which the average is computed.
         :return: SasData object for plotting.
         """
         self.validate_and_assign_data(data2d)
@@ -277,11 +283,17 @@ def __call__(self, data2d: SasData | None = None) -> SasData:
         )
         qy_data, intensity, error = directional_average(data=self.data, err_data=self.err_data)
 
-        data_contents = {
-            "Q": Quantity(qy_data, data2d._data_contents["Qy"].units, None),
-            "I": Quantity(intensity, data2d.ordinate.units, error),
-        }
-        return SasData(f"{data2d.name}: Slab Y Average", data_contents, one_dim, data2d.metadata)
+        ordinate = Quantity(intensity, data2d.ordinate.units, error)
+        abscissa = Abscissa.determine([Quantity(qy_data, data2d.abscissae.axes[1].units, None)], ordinate)
+
+        return SasData(
+            name=f"{data2d.name}: Slab Y Average",
+            ordinate=ordinate,
+            abscissae=abscissa,
+            mask=data2d.mask,
+            dependents=[data2d],
+            metadata=data2d.metadata
+        )
 
 
 class CircularAverage(PolarROI):
@@ -291,8 +303,8 @@ class CircularAverage(PolarROI):
     This class is initialised by specifying lower and upper limits on the
     magnitude of Q values to consider during the averaging, though currently
     SasView always calls this class using the full range of data. When called,
-    this class is supplied with a SasData object. It returns a SasData object
-    where intensity is given as a function of Q only.
+    this class is supplied with a SasMeasurement object.
+    It returns a SasData object where intensity is given as a function of Q only.
     """
 
     def __init__(
@@ -313,13 +325,13 @@ def __init__(
         self.nbins: int = nbins
         self.base: float | None = base
 
-    def __call__(self, data2D: SasData, ismask: bool = False) -> SasData:
+    def __call__(self, data2D: SasMeasurement, ismask: bool = False) -> SasData:
         """
         Perform circular averaging on the data. Uses DirectionalAverage for
         bin construction and weights, and computes dx (d_q) using get_dq_data
         averaged with those weights so behavior matches the legacy implementation.
 
-        :param data2D: SasData object
+        :param data2D: SasMeasurement object
         :param ismask: If True, respect data2D.mask (skip masked points). If False, ignore mask.
         :return: SasData object with x (bin centers), y (intensity), dy and dx (if available)
         """
@@ -329,8 +341,8 @@ def __call__(self, data2D: SasData, ismask: bool = False) -> SasData:
             raise RuntimeError(f"Circular averaging: invalid q_data: {data2D.q_data}")
 
         data = data2D.ordinate.value[finite_mask]
-        qx = data2D._data_contents["Qx"].value[finite_mask]
-        qy = data2D._data_contents["Qy"].value[finite_mask]
+        qx = data2D.abscissae.axes[0].value[finite_mask]
+        qy = data2D.abscissae.axes[1].value[finite_mask]
         q = np.sqrt(qx**2 + qy**2)
         err = np.sqrt(data2D.ordinate.variance.value)[finite_mask]
         mask = (data2D.mask if data2D.mask is not None else np.ones_like(data2D.ordinate.value, dtype=bool))[finite_mask]
@@ -349,7 +361,7 @@ def __call__(self, data2D: SasData, ismask: bool = False) -> SasData:
 
         # Prepare dq_data if available, aligned to the finite mask and selection
         dq_vals = None
-        if data2D._data_contents["Qx"].has_variance and data2D._data_contents["Qy"].has_variance:
+        if data2D.abscissae.axes[0].has_error and data2D.abscissae.axes[1].has_error:
             dq_full = get_dq_data(data2D)  # already uses np.isfinite(data2D.data)
             dq_vals = dq_full[sel]
 
@@ -386,12 +398,17 @@ def __call__(self, data2D: SasData, ismask: bool = False) -> SasData:
         else:
             dQ = None
 
-        data_contents = {
-            "Q": Quantity(x, data2D._data_contents["Qx"].units, dQ),
-            "I": Quantity(intensity, data2D.ordinate.units, error),
-        }
-        return SasData(f"{data2D.name}: Circular Average", data_contents, one_dim, data2D.metadata)
+        ordinate = Quantity(intensity, data2D.ordinate.units, error)
+        abscissa = Abscissa.determine([Quantity(x, data2D.abscissae.axes[0].units, dQ)], ordinate)
 
+        return SasData(
+            name=f"{data2D.name}: Circular Average",
+            ordinate=ordinate,
+            abscissae=abscissa,
+            mask=data2D.mask,
+            dependents=[data2D],
+            metadata=data2D.metadata
+        )
 
 
 class Ring(PolarROI):
@@ -400,8 +417,9 @@ class Ring(PolarROI):
 
     This class is initialised by specifying lower and upper limits on the
     magnitude of Q values to consider during the averaging. When called,
-    this class is supplied with a SasData object. It returns a SasData object
-    which gives intensity as a function of the angle from the positive x-axis, φ, only.
+    this class is supplied with a SasMeasurement object.
+    It returns a SasData object which gives intensity as a
+    function of the angle from the positive x-axis, φ, only.
     """
 
     def __init__(
@@ -430,16 +448,15 @@ def __call__(self, data2D: SasData) -> SasData:
         Apply the ring to the data set.
         Returns the angular distribution for a given q range
 
-        :param data2D: SasData object
+        :param data2D: SasMeasurement object
 
         :return: SasData object
         """
         if not isinstance(data2D, SasData):
             msg = "Data supplied for ring averaging must be of type SasData."
             raise RuntimeError(msg)
-        if not ("Qx" in data2D._data_contents and
-                "Qy" in data2D._data_contents):
-            msg = "SasData object for ring averaging must contain 'Qx' and 'Qy' data."
+        if data2D.abscissae.dimensionality < 2:
+            msg = "SasData object for ring averaging must contain at least two dimensions in Q."
             raise RuntimeError(msg)
 
         # Get data
@@ -447,8 +464,8 @@ def __call__(self, data2D: SasData) -> SasData:
 
         data = data2D.ordinate.value[valid_data]
         err_data = np.sqrt(data2D.ordinate.variance.value)[valid_data]
-        qx_data = data2D._data_contents["Qx"].value[valid_data]
-        qy_data = data2D._data_contents["Qy"].value[valid_data]
+        qx_data = data2D.abscissae.axes[0].value[valid_data]
+        qy_data = data2D.abscissae.axes[1].value[valid_data]
         q_data = np.sqrt(qx_data ** 2 + qy_data ** 2)
         mask_data = (data2D.mask if data2D.mask is not None else np.ones_like(data2D.ordinate.value, dtype=bool))[valid_data]
 
@@ -505,11 +522,17 @@ def __call__(self, data2D: SasData) -> SasData:
             msg = "Average Error: No points inside ROI to average..."
             raise ValueError(msg)
 
-        data_contents = {
-            "Phi": Quantity(phi_values[idx], radians, None),
-            "I": Quantity(phi_bins[idx], data2D.ordinate.units, phi_err[idx]),
-        }
-        return SasData(f"{data2D.name}: Ring Average", data_contents, angle_dim, data2D.metadata)
+        ordinate = Quantity(phi_bins[idx], data2D.ordinate.units, phi_err[idx])
+        abscissa = Abscissa.determine([Quantity(phi_values[idx], radians, None)], ordinate)
+
+        return SasData(
+            name=f"{data2D.name}: Ring Average",
+            ordinate=ordinate,
+            abscissae=abscissa,
+            mask=data2D.mask,
+            dependents=[data2D],
+            metadata=data2D.metadata
+        )
 
 
 class SectorQ(PolarROI):
@@ -530,8 +553,8 @@ class SectorQ(PolarROI):
     the data from the two regions are graphed separeately, with the secondary
     ROI data labelled using negative Q values.
 
-    When called, this class is supplied with a SasData object. It returns a
-    SasData object where intensity is given as a function of Q only.
+    When called, this class is supplied with a SasMeasurement object.
+    It returns a SasData object where intensity is given as a function of Q only.
     """
 
     def __init__(
@@ -559,16 +582,16 @@ def __init__(
         self.fold: bool = fold
         self.base: float | None = base
 
-    def __call__(self, data2d: SasData | None = None) -> SasData:
+    def __call__(self, data2d: SasMeasurement | None = None) -> SasData:
         """
         Compute the 1D average of 2D data, projecting along the Q_y axis.
 
-        :param data2d: The SasData object for which the average is computed.
+        :param data2d: The SasMeasurement object for which the average is computed.
         :return: SasData object for plotting.
         """
         self.validate_and_assign_data(data2d)
 
-        # Detect legacy phi convention (atan2 + pi -> values in [0, 2pi))
+        # Detect legacy phi convention (atan2 + pi -> SasDataalues in [0, 2pi))
         try:
             min_phi = np.nanmin(self.phi_data)
         except Exception:
@@ -641,22 +664,28 @@ def __call__(self, data2d: SasData | None = None) -> SasData:
 
             finite = np.isfinite(average_intensity)
 
-            data_contents = {
-                "Q": Quantity(combined_q[finite], data2d._data_contents["Qx"].units, None),
-                "I": Quantity(average_intensity[finite], data2d.ordinate.units, combined_err[finite]),
-            }
+            average_q = Quantity(combined_q[finite], data2d.abscissae.axes[0].units, None)
+            average_I = Quantity(average_intensity[finite], data2d.ordinate.units, combined_err[finite])
         else:
             # The secondary ROI is labelled with negative Q values.
             combined_q = np.append(np.flip(-1 * secondary_q), primary_q)
             combined_intensity = np.append(np.flip(secondary_I), primary_I)
             combined_error = np.append(np.flip(secondary_err), primary_err)
 
-            data_contents = {
-                "Q": Quantity(combined_q, data2d._data_contents["Qx"].units, None),
-                "I": Quantity(combined_intensity, data2d.ordinate.units, combined_error),
-            }
+            average_q = Quantity(combined_q, data2d.abscissae.axes[0].units, None)
+            average_I = Quantity(combined_intensity, data2d.ordinate.units, combined_error)
 
-        return SasData(f"{data2d.name}:SectorQ Average", data_contents, one_dim, data2d.metadata)
+        ordinate = average_I
+        abscissa = Abscissa.determine([average_q], ordinate)
+
+        return SasData(
+            name=f"{data2d.name}: SectorQ Average",
+            ordinate=ordinate,
+            abscissae=abscissa,
+            mask=data2d.mask,
+            dependents=[data2d],
+            metadata=data2d.metadata
+        )
 
 
 class WedgeQ(PolarROI):
@@ -669,8 +698,8 @@ class WedgeQ(PolarROI):
 
     This class is initialised by specifying lower and upper limits on both the
     magnitude of Q and the angle φ. When called, this class is supplied with a
-    SasData object. It returns a sasData object where intensity is given as a
-    function of Q only.
+    SasMeasurement object. It returns a SasData object where 
+    intensity is given as a function of Q only.
     """
 
     def __init__(
@@ -693,16 +722,16 @@ def __init__(
         self.nbins: int = nbins
         self.base: float | None = base
 
-    def __call__(self, data2d: SasData | None = None) -> SasData:
+    def __call__(self, data2d: SasMeasurement | None = None) -> SasData:
         """
         Compute the 1D average of 2D data, projecting along the Q_y axis.
 
-        :param data2d: The SasData object for which the average is computed.
+        :param data2d: The SasMeasurement object for which the average is computed.
         :return: SasData object for plotting.
         """
         self.validate_and_assign_data(data2d)
 
-        # Detect legacy phi convention (atan2 + pi -> values in [0, 2pi))
+        # Detect legacy phi convention (atan2 + pi -> SasDataalues in [0, 2pi))
         try:
             min_phi = np.nanmin(self.phi_data)
         except Exception:
@@ -738,11 +767,17 @@ def __call__(self, data2d: SasData | None = None) -> SasData:
 
         q_data, intensity, error = directional_average(data=self.data, err_data=self.err_data)
 
-        data_contents = {
-            "Q": Quantity(q_data, data2d._data_contents["Qx"].units, None),
-            "I": Quantity(intensity, data2d.ordinate.units, error),
-        }
-        return SasData(f"{data2d.name}: Wedge Q Average", data_contents, one_dim, data2d.metadata)
+        ordinate = Quantity(intensity, data2d.ordinate.units, error)
+        abscissa = Abscissa.determine([Quantity(q_data, data2d.abscissae.axes[0].units, None)], ordinate)
+
+        return SasData(
+            name=f"{data2d.name}: Wedge Q Average",
+            ordinate=ordinate,
+            abscissae=abscissa,
+            mask=data2d.mask,
+            dependents=[data2d],
+            metadata=data2d.metadata
+        )
 
 
 class WedgePhi(PolarROI):
@@ -754,8 +789,9 @@ class WedgePhi(PolarROI):
 
     This class is initialised by specifying lower and upper limits on both the
     magnitude of Q and the angle φ, measured anticlockwise from the positive
-    x-axis. When called, this class is supplied with a SasData object. It returns
-    a SasData object where intensity is given as a function of φ only.
+    x-axis. When called, this class is supplied with a SasMeasurement object.
+    It returns a SasData object where intensity is given as a 
+    function of φ only.
     """
 
     def __init__(
@@ -779,16 +815,16 @@ def __init__(
         self.nbins: int = nbins
         self.base: float | None = base
 
-    def __call__(self, data2d: SasData | None = None) -> SasData:
+    def __call__(self, data2d: SasMeasurement | None = None) -> SasData:
         """
         Compute the 1D average of 2D data, projecting along the Q_y axis.
 
-        :param data2d: The SasData object for which the average is computed.
+        :param data2d: The SasMeasurement object for which the average is computed.
         :return: SasData object for plotting.
         """
         self.validate_and_assign_data(data2d)
 
-        # Detect legacy phi convention (atan2 + pi -> values in [0, 2pi))
+        # Detect legacy phi convention (atan2 + pi -> SasDataalues in [0, 2pi))
         try:
             min_phi = np.nanmin(self.phi_data)
         except Exception:
@@ -846,11 +882,17 @@ def __call__(self, data2d: SasData | None = None) -> SasData:
         phi_centers = full_phi[populated] + directional_average.bin_widths[populated] / 2.0
 
         # intensity and error returned by DirectionalAverage are already filtered to the populated/finite bins
-        data_contents = {
-            "Phi": Quantity(phi_centers, radians, None),
-            "I": Quantity(intensity, data2d.ordinate.units, error),
-        }
-        return SasData(f"{data2d.name}: Wedge Phi Average", data_contents, angle_dim, data2d.metadata)
+        ordinate = Quantity(intensity, data2d.ordinate.units, error)
+        abscissa = Abscissa.determine([Quantity(phi_centers, radians, None)], ordinate)
+
+        return SasData(
+            name=f"{data2d.name}: Wedge Phi Average",
+            ordinate=ordinate,
+            abscissae=abscissa,
+            mask=data2d.mask,
+            dependents=[data2d],
+            metadata=data2d.metadata
+        )
 
 
 class SectorPhi(WedgePhi):
diff --git a/sasdata/data_util/manipulations.py b/sasdata/data_util/manipulations.py
index 491ebf45..074f1b95 100644
--- a/sasdata/data_util/manipulations.py
+++ b/sasdata/data_util/manipulations.py
@@ -23,6 +23,7 @@
 
 import numpy as np
 
+from sasdata.abscissa import Abscissa
 from sasdata.data import SasData
 
 ################################################################################
@@ -48,7 +49,6 @@
 )
 from sasdata.dataloader.data_info import Data2D
 from sasdata.dataloader.data_info import reader2D_converter as _di_reader2D_converter
-from sasdata.dataset_types import one_dim
 from sasdata.quantities.constants import Pi, TwoPi
 from sasdata.quantities.quantity import Quantity
 
@@ -56,23 +56,6 @@
      "sasdata.data_util.averaging.", DeprecationWarning, stacklevel=2)
 
 
-def deduce_qz(qx: float, qy: float, wavelength: float) -> float:
-    """
-    If you know qx, qy, and the wavelength, you can derive qz
-
-    :param qx: qx [inverse length]
-    :param dy: qy [inverse length]
-    :param wavelength: neutron wavelength [length]
-
-    :return: qz
-    """
-
-    k0 = 2*np.pi/wavelength
-    twotheta = np.arcsin((qx**2 + qy**2) / k0)
-    qz = (1 - np.cos(twotheta)) * k0
-    return qz
-
-
 def position_and_wavelength_to_q(dx: float, dy: float, detector_distance: float, wavelength: float) -> float:
     """
     :param dx: x-distance from beam center
@@ -393,12 +376,12 @@ def __call__(self, data2D, ismask=False):
         # Get data W/ finite values
         finite_mask = np.isfinite(data2D.ordinate.value)
         data = data2D.ordinate.value[finite_mask]
-        q_data = np.sqrt(data2D._data_contents["Qx"].value**2 + data2D._data_contents["Qy"].value**2)[finite_mask]
+        q_data = np.sqrt(data2D.abscissae.axes[0].value**2 + data2D.abscissae.axes[1].value**2)[finite_mask]
         err_data = np.sqrt(data2D.ordinate.variance.value)[finite_mask]
         mask_data = (data2D.mask if data2D.mask is not None else np.ones_like(data2D.ordinate.value, dtype=bool))[finite_mask]
 
         dq_data = None
-        if data2D._data_contents["Qx"].has_variance and data2D._data_contents["Qy"].has_variance:
+        if data2D.abscissae.axes[0].has_error and data2D.abscissae.axes[1].has_error:
             dq_data = get_dq_data(data2D)
 
         if len(q_data) == 0:
@@ -480,11 +463,17 @@ def __call__(self, data2D, ismask=False):
             msg = "Average Error: No points inside ROI to average..."
             raise ValueError(msg)
 
-        data_contents = {
-            "Q": Quantity(x[idx], data2D._data_contents["Qx"].units, dQ),
-            "I": Quantity(y[idx], data2D.ordinate.units, err_y[idx]),
-        }
-        return SasData("Circular Average", data_contents, one_dim, data2D.metadata)
+        ordinate = Quantity(y[idx], data2D.ordinate.units, err_y[idx])
+        abscissa = Abscissa.determine([Quantity(x[idx], data2D.abscissae.axes[0].units, dQ)], ordinate)
+
+        return SasData(
+            name="Circular Average",
+            ordinate=ordinate,
+            abscissae=abscissa,
+            mask=data2D.mask,
+            dependents=[data2D],
+            metadata=data2D.metadata
+        )
 
 ################################################################################
 
@@ -541,23 +530,20 @@ def _agv(self, data2D, run='phi'):
 
         :return: SasData object
         """
-        if not ("Qx" in data2D._data_contents and
-                "Qy" in data2D._data_contents):
-            raise RuntimeError("For averaging the SasData object must contain 'Qx' and 'Qy' data.")
+        if data2D.abscissae.dimensionality < 2:
+            raise RuntimeError("For averaging the SasData object must contain at least two dimensions in Q.")
 
         # Get all the data & info
         finite_mask = np.isfinite(data2D.ordinate.value)
         data = data2D.ordinate.value[finite_mask]
         err_data = np.sqrt(data2D.ordinate.variance.value)[finite_mask]
-        qx_data = data2D._data_contents["Qx"].value[finite_mask]
-        qy_data = data2D._data_contents["Qy"].value[finite_mask]
-        q_data = np.sqrt(data2D._data_contents["Qx"].value ** 2 +
-                         data2D._data_contents["Qy"].value ** 2
-                         )[finite_mask]
+        qx_data = data2D.abscissae.axes[0].value[finite_mask]
+        qy_data = data2D.abscissae.axes[1].value[finite_mask]
+        q_data = np.sqrt(data2D.abscissae.axes[0].value**2 + data2D.abscissae.axes[1].value**2)[finite_mask]
         mask_data = (data2D.mask if data2D.mask is not None else np.ones_like(data2D.ordinate.value, dtype=bool))[finite_mask]
 
         dq_data = None
-        if data2D._data_contents["Qx"].has_variance and data2D._data_contents["Qy"].has_variance:
+        if data2D.abscissae.axes[0].has_error and data2D.abscissae.axes[1].has_error:
             dq_data = get_dq_data(data2D)
 
         # set space for 1d outputs
@@ -722,12 +708,17 @@ def _agv(self, data2D, run='phi'):
             msg = "Average Error: No points inside sector of ROI to average..."
             raise ValueError(msg)
 
-        data_contents = {
-            "Q": Quantity(x[idx], data2D._data_contents["Qx"].units, dQ),
-            "I": Quantity(y[idx], data2D.ordinate.units, y_err[idx]),
-        }
-        return SasData("agv", data_contents, one_dim, data2D.metadata)
-
+        ordinate = Quantity(y[idx], data2D.ordinate.units, y_err[idx])
+        abscissa = Abscissa.determine([Quantity(x[idx], data2D.abscissae.axes[0].units, dQ)], ordinate)
+
+        return SasData(
+            name="agv",
+            ordinate=ordinate,
+            abscissae=abscissa,
+            mask=data2D.mask,
+            dependents=[data2D],
+            metadata=data2D.metadata
+        )
 
 
 class SectorPhi(_Sector):
diff --git a/sasdata/data_util/roi.py b/sasdata/data_util/roi.py
index 3c386aac..9e2f032c 100644
--- a/sasdata/data_util/roi.py
+++ b/sasdata/data_util/roi.py
@@ -38,13 +38,14 @@ def validate_and_assign_data(self, data2d: SasData | None = None) -> None:
         if not isinstance(data2d, SasData):
             msg = "Data supplied must be of type SasData."
             raise TypeError(msg)
-        if not ("Qx" in data2d._data_contents and
-                "Qy" in data2d._data_contents):
-            msg = "SasData object must contain 'Qx' and 'Qy' data."
+        if data2d.abscissae.dimensionality < 2:
+            msg = "SasData object must contain at least two dimensions in Q."
             raise TypeError(msg)
         if len(data2d.metadata.instrument.detector) > 1:
-            msg = (f"Invalid number of detectors: {len(data2d.metadata.instrument.detector)}."
-                   "Cannot have more than 1 detector.")
+            msg = (
+                f"Invalid number of detectors: {len(data2d.metadata.instrument.detector)}. "
+                "Cannot have more than 1 detector."
+            )
             raise ValueError(msg)
 
         # Only use data which is finite and not masked off
@@ -58,9 +59,10 @@ def validate_and_assign_data(self, data2d: SasData | None = None) -> None:
         self.data = data2d.ordinate.value[valid_data]
         self.err_data = np.sqrt(data2d.ordinate.variance.value)[valid_data]
 
-        self.qx_data = data2d._data_contents["Qx"].value[valid_data] - self.center_x
-        self.qy_data = data2d._data_contents["Qy"].value[valid_data] - self.center_y
-        self.q_data = np.sqrt(self.qx_data ** 2 + self.qy_data ** 2)
+        # We take the first two dimensions of the abscissae as the Qx and Qy data.
+        self.qx_data = data2d.abscissae.axes[0].value[valid_data] - self.center_x
+        self.qy_data = data2d.abscissae.axes[1].value[valid_data] - self.center_y
+        self.q_data = np.sqrt(self.qx_data**2 + self.qy_data**2)
 
         # Compute phi in the legacy convention: atan2(qy,qx) + pi
         # (legacy code used this origin; keeping it here makes all polar
@@ -72,6 +74,7 @@ def validate_and_assign_data(self, data2d: SasData | None = None) -> None:
         # previous functionality. It's a bit dodgy, so feel free to remove.
         self.err_data[self.err_data == 0] = np.sqrt(np.abs(self.data[self.err_data == 0]))
 
+
 class CartesianROI(GenericROI):
     """Base class for data manipulators with a Cartesian (rectangular) ROI."""
 
@@ -91,14 +94,16 @@ def __init__(self, qx_range: tuple[float, float] = (0.0, 0.0), qy_range: tuple[f
         self.qy_min = qy_min
         self.qy_max = qy_max
 
+
 class PolarROI(GenericROI):
     """Base class for data manipulators with a polar ROI."""
 
-    def __init__(self,
-                 r_range: tuple[float, float],
-                 phi_range: tuple[float, float] = (0.0, TwoPi),
-                 center: tuple[float, float] = (0.0, 0.0)
-                 ) -> None:
+    def __init__(
+        self,
+        r_range: tuple[float, float],
+        phi_range: tuple[float, float] = (0.0, TwoPi),
+        center: tuple[float, float] = (0.0, 0.0),
+    ) -> None:
         """
         Assign the variables used to label the properties of the SasData object.
         Also establish the upper and lower bounds defining the ROI.
@@ -108,7 +113,7 @@ def __init__(self,
 
         Note that Phi is measured anti-clockwise from the positive x-axis.
         """
-        super().__init__(center = center)
+        super().__init__(center=center)
 
         self.phi_data = None
 
diff --git a/sasdata/exceptions.py b/sasdata/exceptions.py
new file mode 100644
index 00000000..77c70a9e
--- /dev/null
+++ b/sasdata/exceptions.py
@@ -0,0 +1,2 @@
+class InterpretationError(Exception):
+    """Error interpreting data"""
diff --git a/sasdata/metadata.py b/sasdata/metadata.py
index 8c665d26..2a14b39c 100644
--- a/sasdata/metadata.py
+++ b/sasdata/metadata.py
@@ -890,3 +890,6 @@ def collect_tags(objs: list[dataclass]) -> TagCollection:
             result.variable.add(term)
 
     return result
+
+class DerivedMetadata(Metadata):
+    pass
diff --git a/sasdata/postprocess.py b/sasdata/postprocess.py
index 1d786609..799a767f 100644
--- a/sasdata/postprocess.py
+++ b/sasdata/postprocess.py
@@ -13,7 +13,6 @@ def fix_mantid_units_error(data: SasData) -> SasData:
     pass
 
 
-
 def apply_fixes(data: SasData, mantid_unit_error=True):
     if mantid_unit_error:
         data = fix_mantid_units_error(data)
@@ -24,34 +23,22 @@ def apply_fixes(data: SasData, mantid_unit_error=True):
 def deduce_qz(data: SasData):
     """Calculates and appends Qz to SasData if Qx, Qy, and wavelength are all present"""
     # if Qz is not already in the dataset, but Qx and Qy are
-    if 'Qz' not in data._data_contents and 'Qx' in data._data_contents and 'Qy' in data._data_contents:
+    if data.abscissae.dimensionality == 2:
         # we start by making the approximation that qz=0
-        data._data_contents['Qz'] = 0*data._data_contents['Qx']
+        qx = data.abscissae.axes[0]
+        qy = data.abscissae.axes[1]
+        qz = np.zeros_like(qx.value)
 
         # now check if metadata has wavelength information
         wavelength = getattr(
-            getattr(
-                getattr(
-                    getattr(data, "metadata", None),
-                    "instrument",
-                    None
-                ),
-                "source",
-                None
-            ),
-            "wavelength",
-            None
+            getattr(getattr(getattr(data, "metadata", None), "instrument", None), "source", None), "wavelength", None
         )
 
         if wavelength is not None:
             # we can deduce the value of qz from qx and qy
             # if we have the wavelength
-            qx = data._data_contents['Qx']
-            qy = data._data_contents['Qy']
-
             # this is how you convert qx, qy, and wavelength to qz
-            k0 = 2*np.pi/wavelength
-            qz = k0-(k0**2-qx**2-qy**2)**(0.5)
-
-            data._data_contents['Qz'] = qz
+            k0 = 2.0 * np.pi / wavelength
+            qz = k0 - (k0**2 - qx**2 - qy**2) ** 0.5
 
+        return qz
diff --git a/sasdata/temp_ascii_reader.py b/sasdata/temp_ascii_reader.py
index 96e8634b..938031e7 100644
--- a/sasdata/temp_ascii_reader.py
+++ b/sasdata/temp_ascii_reader.py
@@ -11,7 +11,7 @@
     bidirectional_pairings,
     pairings,
 )
-from sasdata.data import SasData
+from sasdata.data import SasMeasurement
 from sasdata.dataset_types import DatasetType, one_dim, unit_kinds
 from sasdata.default_units import get_default_unit
 from sasdata.guess import (
@@ -188,11 +188,13 @@ def merge_uncertainties(quantities: dict[str, Quantity]) -> dict[str, Quantity]:
     return new_quantities
 
 
-def load_data(params: AsciiReaderParams) -> list[SasData]:
-    """This loads a series of SasData objects based on the params. The amount of
-    SasData objects loaded will depend on how many filenames are present in the
-    list contained in the params."""
-    loaded_data: list[SasData] = []
+def load_data(params: AsciiReaderParams) -> list[SasMeasurement]:
+    """
+    This loads a series of SasMeasurement objects based on the params.
+    The amount of SasMeasurement objects loaded will depend on how many
+    filenames are present in the list contained in the params.
+    """
+    loaded_data: list[SasMeasurement] = []
     for filename in params.filenames:
         raw_metadata = import_metadata(
             params.metadata.all_file_metadata(path.basename(filename))
@@ -207,7 +209,7 @@ def load_data(params: AsciiReaderParams) -> list[SasData]:
             raw=raw_metadata,
         )
         quantities = load_quantities(params, filename, metadata)
-        data = SasData(
+        data = SasMeasurement(
             path.basename(filename),
             merge_uncertainties(quantities),
             params.dataset_type,
@@ -217,6 +219,6 @@ def load_data(params: AsciiReaderParams) -> list[SasData]:
     return loaded_data
 
 
-def load_data_default_params(filename: str) -> list[SasData]:
+def load_data_default_params(filename: str) -> list[SasMeasurement]:
     params = guess_params_from_filename(filename, guess_dataset_type(filename))
     return load_data(params)
diff --git a/sasdata/temp_hdf5_reader.py b/sasdata/temp_hdf5_reader.py
index e439486a..fc883519 100644
--- a/sasdata/temp_hdf5_reader.py
+++ b/sasdata/temp_hdf5_reader.py
@@ -6,7 +6,7 @@
 from h5py._hl.dataset import Dataset as HDF5Dataset
 from h5py._hl.group import Group as HDF5Group
 
-from sasdata.data import SasData
+from sasdata.data import SasMeasurement
 from sasdata.data_backing import Dataset as SASDataDataset
 from sasdata.data_backing import Group as SASDataGroup
 from sasdata.dataset_types import one_dim, three_dim, two_dim
@@ -404,9 +404,9 @@ def parse_metadata(node : HDF5Group) -> Metadata:
                     raw=raw)
 
 
-def load_data(filename: str) -> dict[str, SasData]:
+def load_data(filename: str) -> dict[str, SasMeasurement]:
     with h5py.File(filename, "r") as f:
-        loaded_data: dict[str, SasData] = {}
+        loaded_data: dict[str, SasMeasurement] = {}
 
         for root_key in f.keys():
             entry = f[root_key]
@@ -439,7 +439,7 @@ def load_data(filename: str) -> dict[str, SasData]:
 
             entry_key = entry.attrs["sasview_key"] if "sasview_key" in entry.attrs else root_key
 
-            loaded_data[entry_key] = SasData(
+            loaded_data[entry_key] = SasMeasurement(
                     name=root_key,
                     dataset_type=dataset_type,
                     data_contents=data_contents,
diff --git a/sasdata/temp_sesans_reader.py b/sasdata/temp_sesans_reader.py
index 23d56159..00094958 100644
--- a/sasdata/temp_sesans_reader.py
+++ b/sasdata/temp_sesans_reader.py
@@ -8,7 +8,7 @@
 
 import numpy as np
 
-from sasdata.data import SasData
+from sasdata.data import SasMeasurement
 from sasdata.data_util.loader_exceptions import FileContentsException
 from sasdata.dataset_types import sesans
 from sasdata.metadata import (
@@ -191,11 +191,11 @@ def parse_data(lines: list[str], kvs: dict[str, str]) -> dict[str, Quantity]:
     return data_contents
 
 
-def parse_sesans(lines: list[str]) -> SasData:
+def parse_sesans(lines: list[str]) -> SasMeasurement:
     version, lines = parse_version(lines)
     metadata, kvs, lines = parse_metadata(lines)
     data_contents = parse_data(lines, kvs)
-    return SasData(
+    return SasMeasurement(
         name="Sesans",
         dataset_type=sesans,
         data_contents=data_contents,
@@ -204,7 +204,7 @@ def parse_sesans(lines: list[str]) -> SasData:
     )
 
 
-def load_data(filename) -> SasData:
+def load_data(filename) -> SasMeasurement:
     with open(filename) as infile:
         lines = infile.readlines()
     return parse_sesans(lines)
diff --git a/sasdata/temp_xml_reader.py b/sasdata/temp_xml_reader.py
index 174fe738..8fce64ed 100644
--- a/sasdata/temp_xml_reader.py
+++ b/sasdata/temp_xml_reader.py
@@ -5,7 +5,7 @@
 from lxml import etree
 
 import sasdata.quantities.unit_parser as unit_parser
-from sasdata.data import SasData
+from sasdata.data import SasMeasurement
 from sasdata.dataset_types import one_dim
 from sasdata.metadata import (
     Aperture,
@@ -282,9 +282,9 @@ def load_raw(node: etree._Element, version: str) -> MetaNode:
     return MetaNode(name=etree.QName(node).localname, attrs=attrib, contents=contents)
 
 
-def load_data(filename: str) -> dict[str, SasData]:
+def load_data(filename: str) -> dict[str, SasMeasurement]:
     """Load scattering data from an XML file"""
-    loaded_data: dict[str, SasData] = {}
+    loaded_data: dict[str, SasMeasurement] = {}
     tree = etree.parse(filename)
     root = tree.getroot()
 
@@ -327,7 +327,7 @@ def load_data(filename: str) -> dict[str, SasData]:
             data = data_set
             break
 
-        loaded_data[name] = SasData(
+        loaded_data[name] = SasMeasurement(
             name=name,
             dataset_type=one_dim,
             data_contents=data,
diff --git a/sasdata/trend.py b/sasdata/trend.py
index 9b1a371a..394814dd 100644
--- a/sasdata/trend.py
+++ b/sasdata/trend.py
@@ -2,13 +2,14 @@
 
 import numpy as np
 
-from sasdata.data import SasData
+from sasdata.data import SasData, SasMeasurement
 from sasdata.data_backing import Dataset, Group
 from sasdata.quantities.quantity import Quantity
 from sasdata.transforms.rebinning import calculate_interpolation_matrix_1d
 
 # Axis strs refer to the name of their associated NamedQuantity.
 
+
 # TODO: This probably shouldn't be here but will keep it here for now.
 # TODO: Not sure how to type hint the return.
 def get_metadatum_from_path(data: SasData, metadata_path: list[str]):
@@ -16,12 +17,12 @@ def get_metadatum_from_path(data: SasData, metadata_path: list[str]):
     for path_item in metadata_path:
         current_item = current_group.children.get(path_item, None)
         if current_item is None or (isinstance(current_item, Dataset) and path_item != metadata_path[-1]):
-            raise ValueError('Path does not lead to valid a metadatum.')
+            raise ValueError("Path does not lead to valid a metadatum.")
         elif isinstance(current_item, Group):
             current_group = current_item
         else:
             return current_item.data
-    raise ValueError('End of path without finding a dataset.')
+    raise ValueError("End of path without finding a dataset.")
 
 
 @dataclass
@@ -40,6 +41,7 @@ def __getitem__(self, item) -> SasData:
             if metadatum == item:
                 return datum
         raise KeyError()
+
     @property
     def trend_axes(self) -> list[float]:
         return [get_metadatum_from_path(datum, self.trend_axis) for datum in self.data]
@@ -77,13 +79,12 @@ def interpolate(self, axis: str) -> "Trend":
                     continue
                 new_quantities[name] = quantity @ mat
 
-            new_datum = SasData(
+            new_datum = SasMeasurement(
                 name=datum.name,
                 data_contents=new_quantities,
                 dataset_type=datum.dataset_type,
                 metadata=datum.metadata,
             )
             new_data.append(new_datum)
-        new_trend = Trend(new_data,
-                          self.trend_axis)
+        new_trend = Trend(new_data, self.trend_axis)
         return new_trend
diff --git a/sasdata/util.py b/sasdata/util.py
index 8decc68b..6e64b257 100644
--- a/sasdata/util.py
+++ b/sasdata/util.py
@@ -1,10 +1,13 @@
 from collections.abc import Callable
 from typing import TypeVar
 
+import numpy as np
+
 T = TypeVar("T")
 
+
 def cache[T](fun: Callable[[], T]):
-    """ Decorator to store values """
+    """Decorator to store values"""
 
     cache_state = [False, None]
 
@@ -16,3 +19,18 @@ def wrapper() -> T:
         return cache_state[1]
 
     return wrapper
+
+
+def is_increasing(data: np.ndarray):
+    """Check if a 1D array is sorted in strictly increasing order"""
+    return np.all(data[1:] > data[:-1])
+
+
+def is_decreasing(data: np.ndarray):
+    """Check if a 1D array is sorted in strictly decreasing order"""
+    return np.all(data[1:] < data[:-1])
+
+
+def is_sorted(data: np.ndarray):
+    """Check if a 1D array is strictly sorted"""
+    return is_increasing(data) or is_decreasing(data)
diff --git a/test/sasdataloader/utest_data_names.py b/test/sasdataloader/utest_data_names.py
index 11dcb677..668d9122 100644
--- a/test/sasdataloader/utest_data_names.py
+++ b/test/sasdataloader/utest_data_names.py
@@ -28,7 +28,7 @@ def local_load(path: str) -> SasData:
     ("ascii_test_1", "::Q:3KrS58TPgclJ1rgyr0VQp3"),
     ("ISIS_1_1", "TK49 c10_SANS:79680:Q:4TghWEoJi6xxhyeDXhS751"),
     ("cansas1d", "Test title:1234:Q:440tNBqdx9jvci6CgjmrmD"),
-    ("MAR07232_rest", "MAR07232_rest_out.dat:2:/sasentry01/sasdata01/Qx:2Y0qTTb054KSJnJaJv0rFl"),
+    ("MAR07232_rest", "MAR07232_rest_out.dat:2:/sasentry01/sasdata01/Qx:37t0tPj1o8oQcQEB3DVUlw"),
     ("simpleexamplefile", "::/sasentry01/sasdata01/Q:uoHMeB8mukElC1uLCy7Sd"),
 ]
 
@@ -39,5 +39,5 @@ def test_quantity_name(x):
     (f, expected) = x
     data = [v for v in local_load(f"data/{f}")][0]
     if data.metadata.title is not None:
-        assert data.abscissae.unique_id.startswith(data.metadata.title)
-    assert data.abscissae.unique_id == expected
+        assert data.abscissae.axes[0].unique_id.startswith(data.metadata.title)
+    assert data.abscissae.axes[0].unique_id == expected
diff --git a/test/sasdataloader/utest_sasdataload.py b/test/sasdataloader/utest_sasdataload.py
index 97cf14e0..4af8b09c 100644
--- a/test/sasdataloader/utest_sasdataload.py
+++ b/test/sasdataloader/utest_sasdataload.py
@@ -12,7 +12,7 @@
 import pytest
 
 import sasdata.quantities.units as units
-from sasdata.data import SasData, SasDataEncoder
+from sasdata.data import SasData, SasDataEncoder, SasMeasurement
 from sasdata.dataset_types import one_dim
 from sasdata.guess import guess_columns
 from sasdata.quantities.quantity import Quantity
@@ -424,7 +424,7 @@ def test_load_file(test_case: BaseTestCase):
             raw = json.loads("".join(infile.readlines()))
             parsed = {}
             for k in raw:
-                parsed[k] = SasData.from_json(raw[k])
+                parsed[k] = SasMeasurement.from_json(raw[k])
 
         for k in combined_data:
             expect = combined_data[k]
@@ -437,7 +437,7 @@ def test_load_file(test_case: BaseTestCase):
 
     if test_case.round_trip:
         bio = io.BytesIO()
-        SasData.save_h5(combined_data, bio)
+        SasMeasurement.save_h5(combined_data, bio)
         bio.seek(0)
 
         result = hdf_load_data(bio)
diff --git a/test/sasmanipulations/helper.py b/test/sasmanipulations/helper.py
index cc1c8025..e21ee8ff 100644
--- a/test/sasmanipulations/helper.py
+++ b/test/sasmanipulations/helper.py
@@ -4,8 +4,8 @@
 import numpy as np
 from scipy import integrate
 
+from sasdata.abscissa import Abscissa
 from sasdata.data import SasData
-from sasdata.dataset_types import two_dim
 from sasdata.metadata import Instrument, Metadata, Source
 from sasdata.quantities.constants import TwoPi
 from sasdata.quantities.quantity import Quantity
@@ -33,12 +33,10 @@ def integrate_1d_output(output, method="simpson"):
     - If output is a SasData-like object with "Q" and "I" -> integrate I(Q)
     - If output is a tuple (result, error[, npoints]) -> return numeric result
     """
-    if (hasattr(output, "_data_contents") and
-        "Q" in output._data_contents and
-        "I" in output._data_contents):
+    if isinstance(output, SasData) and output.abscissae.dimensionality == 1:
         if method == "trapezoid":
-            return integrate.trapezoid(output._data_contents["I"].value, output._data_contents["Q"].value)
-        return integrate.simpson(output._data_contents["I"].value, output._data_contents["Q"].value)
+            return integrate.trapezoid(output.ordinate.value, output.abscissae.axes[0].value)
+        return integrate.simpson(output.ordinate.value, output.abscissae.axes[0].value)
     if isinstance(output, tuple) and len(output) >= 1:
         return output[0]
     raise TypeError("Unsupported averager output type: %r" % type(output))
@@ -113,14 +111,10 @@ def __init__(self, data2d=None, err_data=None):
         self.qmax = 1
 
         # Create a SasData object to use for testing the averagers.
-        data_contents = {
-            "Qx": Quantity(qx_data, per_angstrom),
-            "Qy": Quantity(qy_data, per_angstrom),
-            "I": Quantity(data_flat, per_centimeter),
-            "dI": Quantity(err_flat, per_centimeter)
-        }
-
-        wavelength = Quantity(1., angstroms)
+        ordinate = Quantity(data_flat, per_centimeter, err_flat)
+        abscissa = Abscissa.determine([Quantity(qx_data, per_angstrom), Quantity(qy_data, per_angstrom)], ordinate)
+
+        wavelength = Quantity(1.0, angstroms)
         source = Source(radiation=None,
                         beam_shape=None,
                         beam_size=None,
@@ -131,15 +125,22 @@ def __init__(self, data2d=None, err_data=None):
         instrument = Instrument(collimations=[],
                                 source=source,
                                 detector=[])
-        metadata=Metadata(title=None,
-                        run=[],
-                        definition=None,
-                        process=[],
-                        sample=None,
-                        instrument=instrument,
-                        raw=None)
-
-        self.data = SasData("Matrix Data", data_contents, two_dim, metadata)
+        basic_metadata = Metadata(title=None,
+                                run=[],
+                                definition=None,
+                                process=[],
+                                sample=None,
+                                instrument=instrument,
+                                raw=None)
+
+        self.data = SasData(
+            name="Matrix Data",
+            ordinate=ordinate,
+            abscissae=abscissa,
+            mask=None,
+            dependents=None,
+            metadata=basic_metadata
+            )
 
     def _validate_and_convert_inputs(self, data2d, err_data):
         """Validate inputs and coerce to numpy arrays. Returns (matrix, err_data_or_None)."""
diff --git a/test/sasmanipulations/utest_averaging.py b/test/sasmanipulations/utest_averaging.py
index 5b6b02ea..0b425a34 100644
--- a/test/sasmanipulations/utest_averaging.py
+++ b/test/sasmanipulations/utest_averaging.py
@@ -3,7 +3,7 @@
 
 import numpy as np
 
-from sasdata.data import SasData, sasdata_reader2D_converter
+from sasdata.data import SasMeasurement, sasdata_reader2D_converter
 from sasdata.data_util.manipulations import (
     Boxavg,
     Boxsum,
@@ -72,7 +72,7 @@ def setUp(self):
                           instrument=instrument,
                           raw=None)
 
-        self.data = SasData("Test Averaging", data_contents, two_dim, metadata)
+        self.data = SasMeasurement("Test Averaging", data_contents, two_dim, metadata)
 
         # get_q(dx, dy, det_dist, wavelength) where units are mm,mm,mm,and A
         # respectively.
@@ -102,7 +102,7 @@ def test_ring_flat_distribution(self):
 
         o = r(self.data)
         for i in range(20):
-            self.assertEqual(o._data_contents["I"].value[i], 1.0)
+            self.assertEqual(o.ordinate.value[i], 1.0)
 
     def test_sectorphi_full(self):
         """Test sector averaging."""
@@ -111,7 +111,7 @@ def test_sectorphi_full(self):
         r.nbins_phi = 20
         o = r(self.data)
         for i in range(7):
-            self.assertEqual(o._data_contents["I"].value[i], 1.0)
+            self.assertEqual(o.ordinate.value[i], 1.0)
 
     def test_sectorphi_partial(self):
         """Test sector averaging."""
@@ -123,7 +123,7 @@ def test_sectorphi_partial(self):
         o = r(self.data)
         self.assertEqual(r.phi_max, phi_max)
         for i in range(17):
-            self.assertEqual(o._data_contents["I"].value[i], 1.0)
+            self.assertEqual(o.ordinate.value[i], 1.0)
 
 
 class DataInfoTests(unittest.TestCase):
@@ -156,9 +156,9 @@ def test_ring(self):
         self.assertEqual(len(answer_list), 1)
         for i in range(r.nbins_phi - 1):
             # Current ascii reader implementation assumes file data is "one_dim"
-            self.assertAlmostEqual(o._data_contents["Phi"].value[i], answer._data_contents["Q"].value[i], 4)
-            self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], 4)
-            self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], 4)
+            self.assertAlmostEqual(o.abscissae.axes[0].value[i], answer.abscissae.axes[0].value[i], 4)
+            self.assertAlmostEqual(o.ordinate.value[i], answer.ordinate.value[i], 4)
+            self.assertAlmostEqual(o.ordinate.variance.value[i], answer.ordinate.variance.value[i], 4)
 
     def test_circularavg(self):
         """
@@ -174,9 +174,9 @@ def test_circularavg(self):
         filepath = find('avg_testdata.txt')
         answer = ascii_load_data(filepath)[0]
         for i in range(r.nbins_phi):
-            self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i], delta=1e-4)
-            self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], delta=1e-4)
-            self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], delta=1e-4)
+            self.assertAlmostEqual(o.abscissae.axes[0].value[i], answer.abscissae.axes[0].value[i], delta=1e-4)
+            self.assertAlmostEqual(o.ordinate.value[i], answer.ordinate.value[i], delta=1e-4)
+            self.assertAlmostEqual(o.ordinate.variance.value[i], answer.ordinate.variance.value[i], delta=1e-4)
 
     def test_box(self):
         """
@@ -208,10 +208,10 @@ def test_slabX(self):
 
         filepath = find('slabx_testdata.txt')
         answer = ascii_load_data(filepath)[0]
-        for i in range(len(o._data_contents["Q"].value)):
-            self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i], 4)
-            self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], 4)
-            self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], 4)
+        for i in range(len(o.abscissae.axes[0].value)):
+            self.assertAlmostEqual(o.abscissae.axes[0].value[i], answer.abscissae.axes[0].value[i], 4)
+            self.assertAlmostEqual(o.ordinate.value[i], answer.ordinate.value[i], 4)
+            self.assertAlmostEqual(o.ordinate.variance.value[i], answer.ordinate.variance.value[i], 4)
 
     def test_slabY(self):
         """
@@ -226,10 +226,10 @@ def test_slabY(self):
 
         filepath = find('slaby_testdata.txt')
         answer = ascii_load_data(filepath)[0]
-        for i in range(len(o._data_contents["Q"].value)):
-            self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i], 4)
-            self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], 4)
-            self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], 4)
+        for i in range(len(o.abscissae.axes[0].value)):
+            self.assertAlmostEqual(o.abscissae.axes[0].value[i], answer.abscissae.axes[0].value[i], 4)
+            self.assertAlmostEqual(o.ordinate.value[i], answer.ordinate.value[i], 4)
+            self.assertAlmostEqual(o.ordinate.variance.value[i], answer.ordinate.variance.value[i], 4)
 
     def test_sectorphi_full(self):
         """
@@ -252,10 +252,10 @@ def test_sectorphi_full(self):
 
         filepath = find('ring_testdata.txt')
         answer = ascii_load_data(filepath)[0]
-        for i in range(len(o._data_contents["Q"].value)-1):
-            self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i+1], 4)
-            self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i+1], 4)
-            self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i+1], 4)
+        for i in range(len(o.abscissae.axes[0].value)-1):
+            self.assertAlmostEqual(o.abscissae.axes[0].value[i], answer.abscissae.axes[0].value[i+1], 4)
+            self.assertAlmostEqual(o.ordinate.value[i], answer.ordinate.value[i+1], 4)
+            self.assertAlmostEqual(o.ordinate.variance.value[i], answer.ordinate.variance.value[i+1], 4)
 
     def test_sectorphi_quarter(self):
         """
@@ -269,10 +269,10 @@ def test_sectorphi_quarter(self):
 
         filepath = find('sectorphi_testdata.txt')
         answer = ascii_load_data(filepath)[0]
-        for i in range(len(o._data_contents["Q"].value)):
-            self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i], 4)
-            self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], 4)
-            self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], 4)
+        for i in range(len(o.abscissae.axes[0].value)):
+            self.assertAlmostEqual(o.abscissae.axes[0].value[i], answer.abscissae.axes[0].value[i], 4)
+            self.assertAlmostEqual(o.ordinate.value[i], answer.ordinate.value[i], 4)
+            self.assertAlmostEqual(o.ordinate.variance.value[i], answer.ordinate.variance.value[i], 4)
 
     def test_sectorq_full(self):
         """
@@ -286,10 +286,10 @@ def test_sectorq_full(self):
 
         filepath = find('sectorq_testdata.txt')
         answer = ascii_load_data(filepath)[0]
-        for i in range(len(o._data_contents["Q"].value)):
-            self.assertAlmostEqual(o._data_contents["Q"].value[i], answer._data_contents["Q"].value[i], 4)
-            self.assertAlmostEqual(o._data_contents["I"].value[i], answer._data_contents["I"].value[i], 4)
-            self.assertAlmostEqual(o._data_contents["I"].variance.value[i], answer._data_contents["I"].variance.value[i], 4)
+        for i in range(len(o.abscissae.axes[0].value)):
+            self.assertAlmostEqual(o.abscissae.axes[0].value[i], answer.abscissae.axes[0].value[i], 4)
+            self.assertAlmostEqual(o.ordinate.value[i], answer.ordinate.value[i], 4)
+            self.assertAlmostEqual(o.ordinate.variance.value[i], answer.ordinate.variance.value[i], 4)
 
     def test_sectorq_log(self):
         """
@@ -302,8 +302,8 @@ def test_sectorq_log(self):
         o = r(self.data)
 
         expected_binning = np.logspace(np.log10(0.005), np.log10(0.01), 20, base=10)
-        for i in range(len(o._data_contents["Q"].value)):
-            self.assertAlmostEqual(o._data_contents["Q"].value[i], expected_binning[i], 3)
+        for i in range(len(o.abscissae.axes[0].value)):
+            self.assertAlmostEqual(o.abscissae.axes[0].value[i], expected_binning[i], 3)
 
         # TODO: Test for Y values (o.y)
         # print len(self.data.x_bins)
diff --git a/test/sasmanipulations/utest_averaging_circle.py b/test/sasmanipulations/utest_averaging_circle.py
index adcecb2a..cd4b4be3 100644
--- a/test/sasmanipulations/utest_averaging_circle.py
+++ b/test/sasmanipulations/utest_averaging_circle.py
@@ -99,7 +99,7 @@ def test_circularaverage_averages_circularly(self):
         data1d = circ_object(averager_data.data)
 
         expected_area = test_data.area_under_region(r_min=r_min, r_max=r_max)
-        actual_area = integrate.trapezoid(data1d._data_contents["I"].value, data1d._data_contents["Q"].value)
+        actual_area = integrate.trapezoid(data1d.ordinate.value, data1d.abscissae.axes[0].value)
 
         # This used to be able to pass with a precision of 3 d.p. with the old
         # manipulations.py - I'm not sure why it doesn't anymore.
@@ -169,7 +169,7 @@ def test_ring_averages_azimuthally(self):
         data1d = ring_object(averager_data.data)
 
         expected_area = test_data.area_under_region(r_min=r_min, r_max=r_max)
-        actual_area = integrate.simpson(data1d._data_contents["I"].value, data1d._data_contents["Phi"].value)
+        actual_area = integrate.simpson(data1d.ordinate.value, data1d.abscissae.axes[0].value)
 
         self.assertAlmostEqual(actual_area, expected_area, 1)
 
@@ -241,7 +241,7 @@ def test_sectorq_averaging_without_fold(self):
         expected_area += test_data.area_under_region(r_min=r_min, r_max=r_max,
                                                      phi_min=phi_min+Pi,
                                                      phi_max=phi_max+Pi)
-        actual_area = integrate.simpson(data1d._data_contents["I"].value, data1d._data_contents["Q"].value)
+        actual_area = integrate.simpson(data1d.ordinate.value, data1d.abscissae.axes[0].value)
 
         self.assertAlmostEqual(actual_area, expected_area, 1)
 
@@ -275,7 +275,7 @@ def test_sectorq_averaging_with_fold(self):
                                                      phi_min=phi_min+Pi,
                                                      phi_max=phi_max+Pi)
         expected_area /= 2
-        actual_area = integrate.simpson(data1d._data_contents["I"].value, data1d._data_contents["Q"].value)
+        actual_area = integrate.simpson(data1d.ordinate.value, data1d.abscissae.axes[0].value)
 
         self.assertAlmostEqual(actual_area, expected_area, 1)
 
@@ -326,7 +326,7 @@ def test_wedgeq_averaging(self):
 
         expected_area = test_data.area_under_region(r_min=r_min, r_max=r_max,
                                                     phi_min=phi_min, phi_max=phi_max)
-        actual_area = integrate.simpson(data1d._data_contents["I"].value, data1d._data_contents["Q"].value)
+        actual_area = integrate.simpson(data1d.ordinate.value, data1d.abscissae.axes[0].value)
 
         self.assertAlmostEqual(actual_area, expected_area, 1)
 
@@ -387,7 +387,7 @@ def test_wedgephi_averaging(self):
         expected_area = test_data.area_under_region(r_min=r_min, r_max=r_max,
                                                     phi_min=phi_min, phi_max=phi_max)
 
-        actual_area = integrate.simpson(data1d._data_contents["I"].value, data1d._data_contents["Phi"].value)
+        actual_area = integrate.simpson(data1d.ordinate.value, data1d.abscissae.axes[0].value)
 
         self.assertAlmostEqual(actual_area, expected_area, 1)
 
diff --git a/test/sasmanipulations/utest_averaging_directional.py b/test/sasmanipulations/utest_averaging_directional.py
index 7fbb5607..864151b0 100644
--- a/test/sasmanipulations/utest_averaging_directional.py
+++ b/test/sasmanipulations/utest_averaging_directional.py
@@ -85,8 +85,8 @@ def setUp(self):
         self.bin_width = (self.lims[1] - self.lims[0]) / self.nbins
 
         self.directional_average = \
-            DirectionalAverage(major_axis=self.data2d.data._data_contents["Qx"].value,
-                               minor_axis=self.data2d.data._data_contents["Qy"].value,
+            DirectionalAverage(major_axis=self.data2d.data.abscissae.axes[0].value,
+                               minor_axis=self.data2d.data.abscissae.axes[1].value,
                                lims=(self.lims,self.lims),
                                nbins=self.nbins)
 
@@ -141,8 +141,8 @@ def test_directional_averaging(self):
         the bins.
         """
         x_axis_values, intensity, errors = \
-            self.directional_average(data=self.data2d.data._data_contents["I"].value,
-                                     err_data=self.data2d.data._data_contents["dI"].value)
+            self.directional_average(data=self.data2d.data.ordinate.value,
+                                     err_data=self.data2d.data.ordinate.standard_error.value)
 
         expected_x = self.qx_data[self.in_roi]
         expected_intensity = np.mean(self.qy_data[self.in_roi]) * expected_x**2
@@ -159,8 +159,8 @@ def test_no_points_in_roi(self):
         self.directional_average.minor_lims = (2, 3)
         self.assertRaises(ValueError,
                  self.directional_average,
-                          self.data2d.data._data_contents["I"].value,
-                          self.data2d.data._data_contents["dI"].value)
+                          self.data2d.data.ordinate.value,
+                          self.data2d.data.ordinate.standard_error.value)
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/test/utest_abscissa.py b/test/utest_abscissa.py
new file mode 100644
index 00000000..2de50f1e
--- /dev/null
+++ b/test/utest_abscissa.py
@@ -0,0 +1,77 @@
+import numpy as np
+import pytest
+
+from sasdata.abscissa import Abscissa, GridAbscissa, MeshgridAbscissa, ScatterAbscissa
+from sasdata.exceptions import InterpretationError
+from sasdata.quantities.quantity import Quantity
+from sasdata.quantities.units import none
+
+
+def test_deterimine_1d_grid():
+    """Test that 1D ordered data is grid type"""
+    q = Quantity(np.arange(10), units=none)
+
+    determined = Abscissa.determine([q], q)
+
+    assert isinstance(determined, GridAbscissa)
+
+
+def test_deterimine_1d_scatter():
+    """Test that 1D unordered data is scatter type"""
+    a = Quantity(np.array([1, 2, 3, 4, 5, 0, 9, 8, 7, 6]), units=none)
+    d = Quantity(np.arange(10), units=none)
+
+    determined = Abscissa.determine([a], d)
+
+    assert isinstance(determined, ScatterAbscissa)
+
+
+def test_2D_scatter():
+    """Test the nD scatter case with 2D data"""
+    q = Quantity(np.arange(10), units=none)
+
+    determined = Abscissa.determine([q, q], q)
+
+    assert isinstance(determined, ScatterAbscissa)
+
+
+def test_2D_meshgrid():
+    """Test the meshgrid case with 2x5"""
+    q = Quantity(np.arange(10).reshape(2, 5), units=none)
+
+    determined = Abscissa.determine([q, q], q)
+
+    assert isinstance(determined, MeshgridAbscissa)
+
+
+def test_2D_grid():
+    """Test the nD grid case with 2x5"""
+    a1 = Quantity(np.arange(2), units=none)
+    a2 = Quantity(np.arange(5), units=none)
+
+    d = Quantity(np.arange(10).reshape(2, 5), units=none)
+
+    determined = Abscissa.determine([a1, a2], d)
+
+    assert isinstance(determined, GridAbscissa)
+
+
+def test_2D_grid_axis_error():
+    """Test the nD grid case with bad axes"""
+
+    a1 = Quantity(np.arange(2), units=none)
+    a2 = Quantity(np.arange(5), units=none)
+
+    d = Quantity(np.arange(10).reshape(5, 2, 1), units=none)
+
+    with pytest.raises(InterpretationError):
+        Abscissa.determine([a1, a2], d)
+
+
+def test_2D_meshgrid_error_mismatched_dimensionality():
+    """Test the nD meshgrid case with bad axes"""
+
+    q = Quantity(np.arange(10).reshape(5, 2), units=none)
+
+    with pytest.raises(InterpretationError):
+        Abscissa.determine([q, q, q], q)  # three axes, each 2D
diff --git a/test/utest_new_sasdata.py b/test/utest_new_sasdata.py
deleted file mode 100644
index 5ee073b3..00000000
--- a/test/utest_new_sasdata.py
+++ /dev/null
@@ -1,125 +0,0 @@
-import numpy as np
-
-from sasdata.data import SasData
-from sasdata.data_backing import Group
-from sasdata.dataset_types import angle_dim, one_dim, three_dim, two_dim
-from sasdata.metadata import Instrument, Metadata, Source
-from sasdata.postprocess import deduce_qz
-from sasdata.quantities.constants import Pi
-from sasdata.quantities.quantity import Quantity
-from sasdata.quantities.units import angstroms, per_angstrom, per_centimeter, radians
-
-
-def test_1d():
-    q = [1, 2, 3, 4, 5]
-    i = [5, 4, 3, 2, 1]
-
-    q_quantity = Quantity(np.array(q), per_angstrom)
-    i_quantity = Quantity(np.array(i), per_centimeter)
-
-    data_contents = {
-        'Q': q_quantity,
-        'I': i_quantity
-    }
-
-    data = SasData('TestData', data_contents, one_dim, Group('root', {}), True)
-
-    assert all(data.abscissae.value == np.array(q))
-    assert all(data.ordinate.value == np.array(i))
-
-
-def test_2d():
-    # This could be autogenerated but I am hard coding to reduce the logic in
-    # the test.
-    qx = [1, 1, 1, 2, 2, 2, 3, 3, 3]
-    qy = [1, 2, 3, 1, 2, 3, 1, 2, 3]
-    i = [1, 2, 3]
-
-    qx_quantity = Quantity(np.array(qx), per_angstrom)
-    qy_quantity = Quantity(np.array(qy), per_angstrom)
-    i_quantity = Quantity(np.array(i), per_centimeter)
-
-    data_contents = {
-        'Qx': qx_quantity,
-        'Qy': qy_quantity,
-        'I': i_quantity
-    }
-
-    data = SasData('TestData', data_contents, two_dim, Group('root', {}), True)
-
-    assert all(data.ordinate.value == np.array(i))
-    assert (data.abscissae.value == np.array([[1, 1], [1, 2], [1, 3], [2, 1], [2, 2], [2, 3], [3, 1], [3, 2], [3, 3]])).all().all()
-
-def test_3d():
-    # test base 3D class
-    qx = [1, 1, 1, 2, 2, 2, 3, 3, 3]
-    qy = [1, 2, 3, 1, 2, 3, 1, 2, 3]
-    qz = [0, 1, 0, 1, 0, 1, 0, 1, 0]
-    i = [1, 2, 3]
-
-    qx_quantity = Quantity(np.array(qx), per_angstrom)
-    qy_quantity = Quantity(np.array(qy), per_angstrom)
-    qz_quantity = Quantity(np.array(qz), per_angstrom)
-    i_quantity = Quantity(np.array(i), per_centimeter)
-
-    data_contents = {
-        'Qx': qx_quantity,
-        'Qy': qy_quantity,
-        'Qz': qz_quantity,
-        'I': i_quantity
-    }
-
-    data = SasData('TestData', data_contents, three_dim, Group('root', {}), True)
-
-    assert (data._data_contents['Qx'].value == np.array(qx)).all()
-
-
-
-    # test autogenerated qz from qx, qy, and wavelength
-    wavelength = Quantity(1., angstroms)
-    source = Source(radiation=None,
-                    beam_shape=None,
-                    beam_size=None,
-                    wavelength=wavelength,
-                    wavelength_max=None,
-                    wavelength_min=None,
-                    wavelength_spread=None)
-    instrument = Instrument(collimations=[],
-                            source=source,
-                            detector=[])
-    metadata=Metadata(title=None,
-                      run=[],
-                      definition=None,
-                      process=[],
-                      sample=None,
-                      instrument=instrument,
-                      raw=None)
-
-    data_contents = {
-        'Qx': qx_quantity,
-        'Qy': qy_quantity,
-        'I': i_quantity
-    }
-
-    data = SasData('TestData', data_contents, two_dim, metadata, True)
-
-    deduce_qz(data)
-
-    assert (data._data_contents['Qz'].value != (0*data._data_contents['Qx'].value)).all()
-
-def test_angle():
-    phi = [0.4 * Pi, 0.8 * Pi, 1.2 * Pi, 1.6 * Pi, 2 * Pi]
-    i = [5, 4, 3, 2, 1]
-
-    phi_quantity = Quantity(np.array(phi), radians)
-    i_quantity = Quantity(np.array(i), per_centimeter)
-
-    data_contents = {
-        'Phi': phi_quantity,
-        'I': i_quantity
-    }
-
-    data = SasData('TestData', data_contents, angle_dim, Group('root', {}), True)
-
-    assert all(data.abscissae.value == np.array(phi))
-    assert all(data.ordinate.value == np.array(i))
diff --git a/test/utest_sasmeasurement.py b/test/utest_sasmeasurement.py
new file mode 100644
index 00000000..232ebbd4
--- /dev/null
+++ b/test/utest_sasmeasurement.py
@@ -0,0 +1,135 @@
+import numpy as np
+import pytest
+
+from sasdata.data import SasMeasurement
+from sasdata.dataset_types import angle_dim, one_dim, three_dim, two_dim
+from sasdata.postprocess import deduce_qz
+from sasdata.quantities.constants import Pi
+from sasdata.quantities.quantity import Quantity
+from sasdata.quantities.units import none as unitless
+from sasdata.quantities.units import per_angstrom, per_centimeter, radians
+
+
+@pytest.mark.parametrize(
+    "q, i",
+    [(np.array([1, 2, 3, 4, 5]), np.array([5, 4, 3, 2, 1])), (np.array([2, 4, 1, 5, 3]), np.array([1, 2, 3, 4, 5]))],
+)
+def test_1d(q, i, basic_metadata):
+    q_quantity = Quantity(q, per_angstrom)
+    i_quantity = Quantity(i, per_centimeter)
+
+    data_contents = {"Q": q_quantity, "I": i_quantity}
+
+    data = SasMeasurement("TestData", data_contents, one_dim, basic_metadata, True)
+
+    assert data.abscissae.dimensionality == 1
+    assert all(data.abscissae.axes[0].value == q)
+    assert all(data.ordinate == i_quantity)
+
+
+# This could be autogenerated but I am hard coding to reduce the logic in the test.
+@pytest.mark.parametrize(
+    "qx, qy, i",
+    [
+        (
+            np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3]]),
+            np.array([[1, 2, 3], [1, 2, 3], [1, 2, 3]]),
+            np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+        ),
+        (np.array([1, 2]), np.array([1, 2, 3, 4, 5]), np.array([[1, 0, 1, 0, 1], [0, 1, 0, 1, 0]])),
+    ],
+)
+def test_2d(qx, qy, i, basic_metadata):
+    qx_quantity = Quantity(qx, per_angstrom)
+    qy_quantity = Quantity(qy, per_angstrom)
+    i_quantity = Quantity(i, per_centimeter)
+
+    data_contents = {"Qx": qx_quantity, "Qy": qy_quantity, "I": i_quantity}
+
+    data = SasMeasurement("TestData", data_contents, two_dim, basic_metadata, True)
+
+    assert data.abscissae.dimensionality == 2
+    assert (data.abscissae.axes[0].value == qx).all()
+    assert (data.abscissae.axes[1].value == qy).all()
+    assert (data.ordinate == i_quantity).all()
+
+
+@pytest.mark.parametrize(
+    "qx, qy, qz, i",
+    [
+        (
+            np.array([[[1, 1], [1, 1]], [[2, 2], [2, 2]]]),
+            np.array([[[1, 1], [2, 2]], [[1, 1], [2, 2]]]),
+            np.array([[[1, 2], [1, 2]], [[1, 2], [1, 2]]]),
+            np.array([[[1, 0], [1, 0]], [[0, 1], [0, 1]]]),
+        ),
+        (
+            np.array([1, 2]),
+            np.array([1, 2]),
+            np.array([1, 2, 3, 4, 5]),
+            np.array([[[1, 0, 1, 0, 1], [1, 0, 1, 0, 1]], [[0, 1, 0, 1, 0], [0, 1, 0, 1, 0]]]),
+        ),
+    ],
+)
+def test_3d(qx, qy, qz, i, basic_metadata):
+    qx_quantity = Quantity(qx, per_angstrom)
+    qy_quantity = Quantity(qy, per_angstrom)
+    qz_quantity = Quantity(qz, per_angstrom)
+    i_quantity = Quantity(i, per_centimeter)
+
+    data_contents = {"Qx": qx_quantity, "Qy": qy_quantity, "Qz": qz_quantity, "I": i_quantity}
+
+    data = SasMeasurement("TestData", data_contents, three_dim, basic_metadata, True)
+
+    assert data.abscissae.dimensionality == 3
+    assert (data.abscissae.axes[0].value == qx).all()
+    assert (data.abscissae.axes[1].value == qy).all()
+    assert (data.abscissae.axes[2].value == qz).all()
+    assert (data.ordinate == i_quantity).all()
+
+
+def test_deduce_qz(basic_metadata):
+    root_two_pi = np.sqrt(2) * np.pi
+    qx = np.array(
+        [
+            [root_two_pi, root_two_pi, root_two_pi],
+            [root_two_pi, root_two_pi, root_two_pi],
+            [root_two_pi, root_two_pi, root_two_pi],
+        ]
+    )
+    qy = np.array(
+        [
+            [root_two_pi, root_two_pi, root_two_pi],
+            [root_two_pi, root_two_pi, root_two_pi],
+            [root_two_pi, root_two_pi, root_two_pi],
+        ]
+    )
+    i = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+
+    qx_quantity = Quantity(qx, unitless) / basic_metadata.instrument.source.wavelength
+    qy_quantity = Quantity(qy, unitless) / basic_metadata.instrument.source.wavelength
+    i_quantity = Quantity(i, per_centimeter)
+
+    data_contents = {"Qx": qx_quantity, "Qy": qy_quantity, "I": i_quantity}
+
+    data = SasMeasurement("TestData", data_contents, two_dim, basic_metadata, True)
+    qz = deduce_qz(data)
+
+    assert qz.value.shape == data.abscissae.axes[0].value.shape
+    assert (qz.value == 2 * np.pi).all()
+
+
+def test_angle(basic_metadata):
+    phi = np.array([0.4 * Pi, 0.8 * Pi, 1.2 * Pi, 1.6 * Pi, 2 * Pi])
+    i = np.array([5, 4, 3, 2, 1])
+
+    phi_quantity = Quantity(phi, radians)
+    i_quantity = Quantity(i, per_centimeter)
+
+    data_contents = {"Phi": phi_quantity, "I": i_quantity}
+
+    data = SasMeasurement("TestData", data_contents, angle_dim, basic_metadata, True)
+
+    assert data.abscissae.dimensionality == 1
+    assert all(data.abscissae.axes[0].value == phi)
+    assert all(data.ordinate == i_quantity)