diff --git a/tests/test_spatial.py b/tests/test_spatial.py
index 1a486c5b..fe0361cd 100644
--- a/tests/test_spatial.py
+++ b/tests/test_spatial.py
@@ -180,6 +180,24 @@ def test_spatial_average_for_lat_region(self):
 
         assert result.identical(expected)
 
+    def test_spatial_average_for_domain_wrapping_p_meridian_non_cf_conventions(
+        self,
+    ):
+        ds = self.ds.copy()
+
+        # get spatial average for original dataset
+        ref = ds.spatial.average("ts").ts
+
+        # change first bound from -0.9375 to 359.0625
+        lon_bnds = ds.lon_bnds.copy()
+        lon_bnds[0, 0] = 359.0625
+        ds["lon_bnds"] = lon_bnds
+
+        # check spatial average with new (bad) bound
+        result = ds.spatial.average("ts").ts
+
+        assert result.identical(ref)
+
     @requires_dask
     def test_spatial_average_for_lat_region_and_keep_weights_with_dask(self):
         ds = self.ds.copy().chunk(2)
@@ -269,26 +287,17 @@ def setup(self):
             decode_times=True, cf_compliant=False, has_bounds=True
         )
 
-    def test_bounds_reordered_when_upper_indexed_first(self):
+    def test_value_error_thrown_for_multiple_out_of_order_lon_bounds(self):
         domain_bounds = xr.DataArray(
             name="lon_bnds",
-            data=np.array(
-                [[-89.375, -90], [0.0, -89.375], [0.0, 89.375], [89.375, 90]]
-            ),
+            data=np.array([[3, 1], [5, 3], [5, 7], [7, 9]]),
             coords={"lat": self.ds.lat},
             dims=["lat", "bnds"],
         )
-        result = self.ds.spatial._force_domain_order_low_to_high(domain_bounds)
-
-        expected_domain_bounds = xr.DataArray(
-            name="lon_bnds",
-            data=np.array(
-                [[-90, -89.375], [-89.375, 0.0], [0.0, 89.375], [89.375, 90]]
-            ),
-            coords={"lat": self.ds.lat},
-            dims=["lat", "bnds"],
-        )
-        assert result.identical(expected_domain_bounds)
+        # Check _get_longitude_weights raises error when there are
+        # > 1 out-of-order bounds for the dataset.
+        with pytest.raises(ValueError):
+            self.ds.spatial._get_longitude_weights(domain_bounds, region_bounds=None)
 
     def test_raises_error_if_dataset_has_multiple_bounds_variables_for_an_axis(self):
         ds = self.ds.copy()
diff --git a/xcdat/spatial.py b/xcdat/spatial.py
index c51af518..2c50595a 100644
--- a/xcdat/spatial.py
+++ b/xcdat/spatial.py
@@ -281,12 +281,6 @@ def get_weights(
                     "to reference a specific data variable's axis bounds."
                 )
 
-            # The logic for generating longitude weights depends on the
-            # bounds being ordered such that d_bounds[:, 0] < d_bounds[:, 1].
-            # They are re-ordered (if need be) for the purpose of creating
-            # weights.
-            d_bounds = self._force_domain_order_low_to_high(d_bounds)
-
             r_bounds = axis_bounds[key]["region"]
 
             weights = axis_bounds[key]["weights_method"](d_bounds, r_bounds)
@@ -324,33 +318,6 @@ def _validate_axis_arg(self, axis: List[SpatialAxis]):
             # Check the axis coordinate variable exists in the Dataset.
             get_dim_coords(self._dataset, key)
 
-    def _force_domain_order_low_to_high(self, domain_bounds: xr.DataArray):
-        """Reorders the ``domain_bounds`` low-to-high.
-
-        This method ensures all lower bound values are less than the upper bound
-        values (``domain_bounds[:, 1] < domain_bounds[:, 1]``).
-
-        Parameters
-        ----------
-        domain_bounds: xr.DataArray
-            The bounds of an axis.
-
-        Returns
-        ------
-        xr.DataArray
-            The bounds of an axis (re-ordered if applicable).
-        """
-        index_bad_cells = np.where(domain_bounds[:, 1] - domain_bounds[:, 0] < 0)[0]
-
-        if len(index_bad_cells) > 0:
-            new_domain_bounds = domain_bounds.copy()
-            new_domain_bounds[index_bad_cells, 0] = domain_bounds[index_bad_cells, 1]
-            new_domain_bounds[index_bad_cells, 1] = domain_bounds[index_bad_cells, 0]
-
-            return new_domain_bounds
-
-        return domain_bounds
-
     def _validate_region_bounds(self, axis: SpatialAxis, bounds: RegionAxisBounds):
         """Validates the ``bounds`` arg based on a set of criteria.
 
@@ -441,12 +408,29 @@ def _get_longitude_weights(
         -------
         xr.DataArray
             The longitude axis weights.
+
+        Raises
+        ------
+        ValueError
+            If the there are multiple instances in which the
+            domain_bounds[:, 0] > domain_bounds[:, 1]
         """
         p_meridian_index: Optional[np.ndarray] = None
         d_bounds = domain_bounds.copy()
 
+        pm_cells = np.where(domain_bounds[:, 1] - domain_bounds[:, 0] < 0)[0]
+        if len(pm_cells) > 1:
+            raise ValueError(
+                "More than one longitude bound is out of order. Only one bound "
+                "value spanning the prime meridian is permitted in data on "
+                "a rectilinear grid."
+            )
+        d_bounds: xr.DataArray = self._swap_lon_axis(d_bounds, to=360)  # type: ignore
+        p_meridian_index = _get_prime_meridian_index(d_bounds)
+        if p_meridian_index is not None:
+            d_bounds = _align_lon_bounds_to_360(d_bounds, p_meridian_index)
+
         if region_bounds is not None:
-            d_bounds: xr.DataArray = self._swap_lon_axis(d_bounds, to=360)  # type: ignore
             r_bounds: np.ndarray = self._swap_lon_axis(
                 region_bounds, to=360
             )  # type:ignore
@@ -455,10 +439,6 @@ def _get_longitude_weights(
             if is_region_circular:
                 r_bounds = np.array([0.0, 360.0])
 
-            p_meridian_index = _get_prime_meridian_index(d_bounds)
-            if p_meridian_index is not None:
-                d_bounds = _align_lon_bounds_to_360(d_bounds, p_meridian_index)
-
             d_bounds = self._scale_domain_to_region(d_bounds, r_bounds)
 
         weights = self._calculate_weights(d_bounds)