diff --git a/sasdata/abscissa.py b/sasdata/abscissa.py
new file mode 100644
index 00000000..a3826309
--- /dev/null
+++ b/sasdata/abscissa.py
@@ -0,0 +1,124 @@
+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 73e64e61..059c9b1f 100644
--- a/sasdata/data.py
+++ b/sasdata/data.py
@@ -6,11 +6,68 @@
 
 from sasdata import dataset_types
 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: list[Quantity],
+                 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) -> list[Quantity]:
+        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: list[Quantity],
+                 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,
@@ -32,7 +89,7 @@ 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.
@@ -52,32 +109,15 @@ def ordinate(self) -> Quantity:
     def abscissae(self) -> Quantity:
         match self.dataset_type:
             case dataset_types.one_dim:
-                return self._data_contents["Q"]
+                return [self._data_contents["Q"]]
             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 [self._data_contents["Qx"], self._data_contents["Qy"]]
             case dataset_types.angle_dim:
-                return self._data_contents["Phi"]
+                return [self._data_contents["Phi"]]
             case dataset_types.sesans:
-                return self._data_contents["SpinEchoLength"]
+                return [self._data_contents["SpinEchoLength"]]
             case _:
-                None
+                return None
 
     def __getitem__(self, item: str):
         return self._data_contents[item]
@@ -96,7 +136,7 @@ def summary(self, indent="  "):
 
     @staticmethod
     def from_json(obj):
-        return SasData(
+        return SasMeasurement(
             name=obj["name"],
             dataset_type=DatasetType(
                 name=obj["type"]["name"],
@@ -119,7 +159,6 @@ 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:
@@ -130,7 +169,6 @@ def save_h5(data: dict[str, typing.Self], path: str | typing.BinaryIO):
                 data._save_h5(sasentry)
 
 
-
 class SasDataEncoder(MetadataEncoder):
     def default(self, obj):
         match obj:
@@ -174,7 +212,7 @@ 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()
@@ -184,6 +222,6 @@ def sasdata_reader2D_converter(data2d: SasData | None = None) -> SasData:
     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..f2a6ebfa 100644
--- a/sasdata/data_util/averaging.py
+++ b/sasdata/data_util/averaging.py
@@ -5,17 +5,16 @@
 import numpy as np
 import numpy.typing as npt
 
-from sasdata.data import SasData
+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 +24,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[0].value**2 + data2d.abscissae[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[0].variance.value)
+    dqy_data = np.sqrt(data2d.abscissae[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 +77,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 +125,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 +142,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 +172,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 +203,14 @@ 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)
+        return SasData(
+            name=f"{data2d.name}: Slab X Average",
+            ordinate=Quantity(intensity, data2d.ordinate.units, error),
+            abscissae=[Quantity(qx_data, data2d.abscissae[0].units, None)],
+            mask=data2d.mask,
+            dependents=[data2d],
+            metadata=data2d.metadata
+        )
 
 
 class SlabY(CartesianROI):
@@ -216,7 +218,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 +249,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 +279,14 @@ 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)
+        return SasData(
+            name=f"{data2d.name}: Slab Y Average",
+            ordinate=Quantity(intensity, data2d.ordinate.units, error),
+            abscissae=[Quantity(qy_data, data2d.abscissae[1].units, None)],
+            mask=data2d.mask,
+            dependents=[data2d],
+            metadata=data2d.metadata
+        )
 
 
 class CircularAverage(PolarROI):
@@ -291,8 +296,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 +318,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 +334,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[0].value[finite_mask]
+        qy = data2D.abscissae[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 +354,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[0].has_error and data2D.abscissae[1].has_error:
             dq_full = get_dq_data(data2D)  # already uses np.isfinite(data2D.data)
             dq_vals = dq_full[sel]
 
@@ -386,12 +391,14 @@ 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)
-
+        return SasData(
+            name=f"{data2D.name}: Circular Average",
+            ordinate=Quantity(intensity, data2D.ordinate.units, error),
+            abscissae=[Quantity(x, data2D.abscissae[0].units, dQ)],
+            mask=data2D.mask,
+            dependents=[data2D],
+            metadata=data2D.metadata
+        )
 
 
 class Ring(PolarROI):
@@ -400,8 +407,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 +438,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 len(data2D.abscissae) < 2:
+            msg = "SasData object for ring averaging must contain at least two dimensions in Q."
             raise RuntimeError(msg)
 
         # Get data
@@ -447,8 +454,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[0].value[valid_data]
+        qy_data = data2D.abscissae[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 +512,14 @@ 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)
+        return SasData(
+            name=f"{data2D.name}: Ring Average",
+            ordinate=Quantity(phi_bins[idx], data2D.ordinate.units, phi_err[idx]),
+            abscissae=[Quantity(phi_values[idx], radians, None)],
+            mask=data2D.mask,
+            dependents=[data2D],
+            metadata=data2D.metadata
+        )
 
 
 class SectorQ(PolarROI):
@@ -530,8 +540,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 +569,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 +651,25 @@ 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[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[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)
+        return SasData(
+            name=f"{data2d.name}: SectorQ Average",
+            ordinate=average_I,
+            abscissae=[average_q],
+            mask=data2d.mask,
+            dependents=[data2d],
+            metadata=data2d.metadata
+        )
 
 
 class WedgeQ(PolarROI):
@@ -669,8 +682,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 +706,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 +751,14 @@ 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)
+        return SasData(
+            name=f"{data2d.name}: Wedge Q Average",
+            ordinate=Quantity(intensity, data2d.ordinate.units, error),
+            abscissae=[Quantity(q_data, data2d.abscissae[0].units, None)],
+            mask=data2d.mask,
+            dependents=[data2d],
+            metadata=data2d.metadata
+        )
 
 
 class WedgePhi(PolarROI):
@@ -754,8 +770,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 +796,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 +863,14 @@ 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)
+        return SasData(
+            name=f"{data2d.name}: Wedge Phi Average",
+            ordinate=Quantity(intensity, data2d.ordinate.units, error),
+            abscissae=[Quantity(phi_centers, radians, None)],
+            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..01253c8a 100644
--- a/sasdata/data_util/manipulations.py
+++ b/sasdata/data_util/manipulations.py
@@ -48,7 +48,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
 
@@ -393,12 +392,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[0].value**2 + data2D.abscissae[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[0].has_error and data2D.abscissae[1].has_error:
             dq_data = get_dq_data(data2D)
 
         if len(q_data) == 0:
@@ -480,11 +479,14 @@ 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)
+        return SasData(
+            name="Circular Average",
+            ordinate=Quantity(y[idx], data2D.ordinate.units, err_y[idx]),
+            abscissae=[Quantity(x[idx], data2D.abscissae[0].units, dQ)],
+            mask=data2D.mask,
+            dependents=[data2D],
+            metadata=data2D.metadata
+        )
 
 ################################################################################
 
@@ -541,23 +543,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 len(data2D.abscissae) != 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[0].value[finite_mask]
+        qy_data = data2D.abscissae[1].value[finite_mask]
+        q_data = np.sqrt(data2D.abscissae[0].value**2 + data2D.abscissae[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[0].has_error and data2D.abscissae[1].has_error:
             dq_data = get_dq_data(data2D)
 
         # set space for 1d outputs
@@ -722,12 +721,14 @@ 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)
-
+        return SasData(
+            name="agv",
+            ordinate=Quantity(y[idx], data2D.ordinate.units, y_err[idx]),
+            abscissae=[Quantity(x[idx], data2D.abscissae[0].units, dQ)],
+            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..414aa557 100644
--- a/sasdata/data_util/roi.py
+++ b/sasdata/data_util/roi.py
@@ -38,9 +38,8 @@ 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 len(data2d.abscissae) < 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)}."
@@ -58,8 +57,9 @@ 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
+        # We take the first two dimensions of the abscissae as the Qx and Qy data.
+        self.qx_data = data2d.abscissae[0].value[valid_data] - self.center_x
+        self.qy_data = data2d.abscissae[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
diff --git a/sasdata/exceptions.py b/sasdata/exceptions.py
new file mode 100644
index 00000000..6fc3e545
--- /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 d53c3102..e20ba062 100644
--- a/sasdata/metadata.py
+++ b/sasdata/metadata.py
@@ -847,3 +847,6 @@ def collect_tags(objs: list[dataclass]) -> TagCollection:
             result.variable.add(term)
 
     return result
+
+class DerivedMetadata(Metadata):
+    pass
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..4ab10729 100644
--- a/sasdata/trend.py
+++ b/sasdata/trend.py
@@ -2,7 +2,7 @@
 
 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
@@ -77,7 +77,7 @@ 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,
diff --git a/sasdata/util.py b/sasdata/util.py
index 8decc68b..5dd0d404 100644
--- a/sasdata/util.py
+++ b/sasdata/util.py
@@ -1,6 +1,8 @@
 from collections.abc import Callable
 from typing import TypeVar
 
+import numpy as np
+
 T = TypeVar("T")
 
 def cache[T](fun: Callable[[], T]):
@@ -16,3 +18,15 @@ 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..2bcfded7 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[0].unique_id.startswith(data.metadata.title)
+    assert data.abscissae[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..844a9169 100644
--- a/test/sasmanipulations/helper.py
+++ b/test/sasmanipulations/helper.py
@@ -5,7 +5,6 @@
 from scipy import integrate
 
 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 +32,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 len(output.abscissae) == 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[0].value)
+        return integrate.simpson(output.ordinate.value, output.abscissae[0].value)
     if isinstance(output, tuple) and len(output) >= 1:
         return output[0]
     raise TypeError("Unsupported averager output type: %r" % type(output))
@@ -113,13 +110,6 @@ 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)
         source = Source(radiation=None,
                         beam_shape=None,
@@ -139,7 +129,14 @@ def __init__(self, data2d=None, err_data=None):
                         instrument=instrument,
                         raw=None)
 
-        self.data = SasData("Matrix Data", data_contents, two_dim, metadata)
+        self.data = SasData(
+            name="Matrix Data",
+            ordinate=Quantity(data_flat, per_centimeter, err_flat),
+            abscissae=[Quantity(qx_data, per_angstrom), Quantity(qy_data, per_angstrom)],
+            mask=None,
+            dependents=None,
+            metadata=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..27378c56 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[0].value[i], answer.abscissae[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[0].value[i], answer.abscissae[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[0].value)):
+            self.assertAlmostEqual(o.abscissae[0].value[i], answer.abscissae[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[0].value)):
+            self.assertAlmostEqual(o.abscissae[0].value[i], answer.abscissae[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[0].value)-1):
+            self.assertAlmostEqual(o.abscissae[0].value[i], answer.abscissae[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[0].value)):
+            self.assertAlmostEqual(o.abscissae[0].value[i], answer.abscissae[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[0].value)):
+            self.assertAlmostEqual(o.abscissae[0].value[i], answer.abscissae[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[0].value)):
+            self.assertAlmostEqual(o.abscissae[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..2209a8af 100644
--- a/test/sasmanipulations/utest_averaging_circle.py
+++ b/test/sasmanipulations/utest_averaging_circle.py
@@ -94,12 +94,13 @@ def test_circularaverage_averages_circularly(self):
         r_min = averager_data.qmax * 0.25
         r_max = averager_data.qmax * 0.75
 
+        print(averager_data.data)
         nbins = test_data.matrix_size
         circ_object = CircularAverage(r_range=(r_min, r_max), nbins=nbins)
         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[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 +170,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[0].value)
 
         self.assertAlmostEqual(actual_area, expected_area, 1)
 
@@ -241,7 +242,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[0].value)
 
         self.assertAlmostEqual(actual_area, expected_area, 1)
 
@@ -275,7 +276,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[0].value)
 
         self.assertAlmostEqual(actual_area, expected_area, 1)
 
@@ -326,7 +327,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[0].value)
 
         self.assertAlmostEqual(actual_area, expected_area, 1)
 
@@ -387,7 +388,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[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..eb266a8f 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[0].value,
+                               minor_axis=self.data2d.data.abscissae[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_asbscissa.py b/test/utest_asbscissa.py
new file mode 100644
index 00000000..3b166313
--- /dev/null
+++ b/test/utest_asbscissa.py
@@ -0,0 +1,72 @@
+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
index 5ee073b3..36ef3ab1 100644
--- a/test/utest_new_sasdata.py
+++ b/test/utest_new_sasdata.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-from sasdata.data import SasData
+from sasdata.data import SasMeasurement
 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
@@ -22,9 +22,9 @@ def test_1d():
         'I': i_quantity
     }
 
-    data = SasData('TestData', data_contents, one_dim, Group('root', {}), True)
+    data = SasMeasurement('TestData', data_contents, one_dim, Group('root', {}), True)
 
-    assert all(data.abscissae.value == np.array(q))
+    assert all(data.abscissae[0].value == np.array(q))
     assert all(data.ordinate.value == np.array(i))
 
 
@@ -45,10 +45,11 @@ def test_2d():
         'I': i_quantity
     }
 
-    data = SasData('TestData', data_contents, two_dim, Group('root', {}), True)
+    data = SasMeasurement('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()
+    assert (data.abscissae[0].value == np.array([1, 1, 1, 2, 2, 2, 3, 3, 3])).all()
+    assert (data.abscissae[1].value == np.array([1, 2, 3, 1, 2, 3, 1, 2, 3])).all()
 
 def test_3d():
     # test base 3D class
@@ -69,7 +70,7 @@ def test_3d():
         'I': i_quantity
     }
 
-    data = SasData('TestData', data_contents, three_dim, Group('root', {}), True)
+    data = SasMeasurement('TestData', data_contents, three_dim, Group('root', {}), True)
 
     assert (data._data_contents['Qx'].value == np.array(qx)).all()
 
@@ -101,7 +102,7 @@ def test_3d():
         'I': i_quantity
     }
 
-    data = SasData('TestData', data_contents, two_dim, metadata, True)
+    data = SasMeasurement('TestData', data_contents, two_dim, metadata, True)
 
     deduce_qz(data)
 
@@ -119,7 +120,7 @@ def test_angle():
         'I': i_quantity
     }
 
-    data = SasData('TestData', data_contents, angle_dim, Group('root', {}), True)
+    data = SasMeasurement('TestData', data_contents, angle_dim, Group('root', {}), True)
 
-    assert all(data.abscissae.value == np.array(phi))
+    assert all(data.abscissae[0].value == np.array(phi))
     assert all(data.ordinate.value == np.array(i))