diff --git a/yt/visualization/plot_window.py b/yt/visualization/plot_window.py
index 78f96f14a28..5ce9bd790ea 100644
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -26,6 +26,7 @@
     validate_moment,
 )
 from yt.geometry.api import Geometry
+from yt.geometry.oct_geometry_handler import OctreeIndex
 from yt.units.unit_object import Unit  # type: ignore
 from yt.units.unit_registry import UnitParseError  # type: ignore
 from yt.units.yt_array import YTArray, YTQuantity
@@ -2494,7 +2495,12 @@ def __init__(
         is_sph_field = finfo.is_sph_field
         particle_datasets = (ParticleDataset, StreamParticlesDataset)
 
-        if isinstance(data_source.ds, particle_datasets) and is_sph_field:
+        dom_width = data_source.ds.domain_width
+        cubic_domain = dom_width.max() == dom_width.min()
+
+        if (isinstance(data_source.ds, particle_datasets) and is_sph_field) or (
+            isinstance(data_source.ds.index, OctreeIndex) and cubic_domain
+        ):
             center_use = parse_center_array(center, ds=data_source.ds, axis=None)
         else:
             center_use = center_rot
diff --git a/yt/visualization/tests/test_offaxisprojection.py b/yt/visualization/tests/test_offaxisprojection.py
index b9839745d7a..46dccf860ae 100644
--- a/yt/visualization/tests/test_offaxisprojection.py
+++ b/yt/visualization/tests/test_offaxisprojection.py
@@ -16,6 +16,7 @@
 from yt.visualization.api import (
     OffAxisProjectionPlot,
     OffAxisSlicePlot,
+    ProjectionPlot,
 )
 from yt.visualization.image_writer import write_projection
 from yt.visualization.volume_rendering.api import off_axis_projection
@@ -210,6 +211,35 @@ def test_field_cut_off_axis_octree():
     assert_equal(np.nanmin(p4rho[p4rho > 0.0]) >= 0.5, True)
 
 
+def test_off_axis_octree():
+    ds = fake_octree_ds()
+    p1 = ProjectionPlot(
+        ds,
+        "x",
+        ("gas", "density"),
+        center=[0.6] * 3,
+        width=0.8,
+        weight_field=("gas", "density"),
+    )
+    p2 = OffAxisProjectionPlot(
+        ds,
+        [1, 0, 0],
+        ("gas", "density"),
+        center=[0.6] * 3,
+        width=0.8,
+        weight_field=("gas", "density"),
+    )
+
+    # Note: due to our implementation, the off-axis projection will have a
+    # slightly blurred cell edges so we can't do an exact comparison
+    v1, v2 = p1.frb["gas", "density"], p2.frb["gas", "density"]
+    diff = (v1 - v2) / (v1 + v2) * 2
+
+    # Make sure the difference is zero-centered with a small standard deviation
+    assert np.mean(diff).max() < 1e-3  # 0.1%: very little bias
+    assert np.std(diff) < 0.05  # <2% error on average
+
+
 def test_offaxis_moment():
     ds = fake_random_ds(64)
 
diff --git a/yt/visualization/volume_rendering/off_axis_projection.py b/yt/visualization/volume_rendering/off_axis_projection.py
index 3c2359a3f80..38c5677ef55 100644
--- a/yt/visualization/volume_rendering/off_axis_projection.py
+++ b/yt/visualization/volume_rendering/off_axis_projection.py
@@ -444,21 +444,23 @@ def temp_weightfield(field, data):
         # We need the width of the plot window in projected coordinates,
         # i.e. we ignore the z-component
         wmax = width[:2].max()
-
-        # Normalize the positions & dx so that they are in the range [-0.5, 0.5]
-        xyz = np.stack(
-            [
-                ((data_source["index", k] - center[i]) / wmax).to("1").d
-                for i, k in enumerate("xyz")
-            ],
-            axis=-1,
+        xyz = data_source.ds.arr(
+            np.zeros((len(data_source[vol.field]), 3)), "code_length"
         )
 
         for idim, periodic in enumerate(data_source.ds.periodicity):
+            axis = data_source.ds.coordinates.axis_order[idim]
+            # Recenter positions w.r.t. center of the plot window
+            xyz[..., idim] = data_source["index", axis] - center[idim]
             if not periodic:
                 continue
-            # Wrap into [-0.5, +0.5]
-            xyz[..., idim] = (xyz[..., idim] + 0.5) % 1 - 0.5
+            # If we have periodic boundaries, we need to wrap the corresponding
+            # coordinates into [-w/2, +w/2]
+            w = data_source.ds.domain_width[idim]
+            xyz[..., idim] = (xyz[..., idim] + w / 2) % w - w / 2
+
+        # Rescale to [-0.5, +0.5]
+        xyz = (xyz / wmax).to("1").d
 
         dx = (data_source["index", "dx"] / wmax).to("1").d