BUG: calculate default depth for off axis projections

pyproject.toml

@@ -200,6 +200,7 @@ test = [
     "pytest-mpl>=0.16.1",
     "sympy!=1.10,!=1.9", # see https://github.com/sympy/sympy/issues/22241
     "imageio!=2.35.0", # see https://github.com/yt-project/yt/issues/4966
+    "contourpy",
 ]
 
 [project.scripts]

yt/visualization/fixed_resolution.py

@@ -639,7 +639,6 @@ def _generate_image_and_mask(self, item) -> None:
                 self.bounds[5] - self.bounds[4],
             )
         )
-        depth = dd.depth[0] if dd.depth is not None else None
         buff = off_axis_projection(
             dd.dd,
             dd.center,
@@ -652,7 +651,7 @@ def _generate_image_and_mask(self, item) -> None:
             no_ghost=dd.no_ghost,
             interpolated=dd.interpolated,
             north_vector=dd.north_vector,
-            depth=depth,
+            depth=dd.depth,
             method=dd.method,
         )
         if self.data_source.moment == 2:

yt/visualization/plot_window.py

@@ -83,7 +83,11 @@ def get_window_parameters(axis, center, width, ds):
 
 
 def get_oblique_window_parameters(
-    normal, center, width, ds, depth=None, get3bounds=False
+    normal,
+    center,
+    width,
+    ds,
+    depth=None,
 ):
     center, display_center = ds.coordinates.sanitize_center(center, axis=None)
     width = ds.coordinates.sanitize_width(normal, width, depth)
@@ -101,15 +105,7 @@ def get_oblique_window_parameters(
 
     w = tuple(el.in_units("code_length") for el in width)
     bounds = tuple(((2 * (i % 2)) - 1) * w[i // 2] / 2 for i in range(len(w) * 2))
-    if get3bounds and depth is None:
-        # off-axis projection, depth not specified
-        # -> set 'large enough' depth using half the box diagonal + margin
-        d2 = ds.domain_width[0].in_units("code_length") ** 2
-        d2 += ds.domain_width[1].in_units("code_length") ** 2
-        d2 += ds.domain_width[2].in_units("code_length") ** 2
-        diag = np.sqrt(d2)
-        bounds = bounds + (-0.51 * diag, 0.51 * diag)
-    return (bounds, center)
+    return bounds, center
 
 
 def get_axes_unit(width, ds):
@@ -2389,7 +2385,8 @@ class OffAxisProjectionPlot(ProjectionPlot, PWViewerMPL):
     depth : A tuple or a float
         A tuple containing the depth to project through and the string
         key of the unit: (width, 'unit'). If set to a float, code units
-        are assumed
+        are assumed. In not set, then a depth equal to the diagonal of
+        the domain width plus a small margin will be used.
     weight_field : string
         The name of the weighting field.  Set to None for no weight.
     max_level: int
@@ -2466,6 +2463,13 @@ def __init__(
                 "currently supported geometries:"
                 f" {self._supported_geometries!r}"
             )
+
+        if depth is None:
+            # off-axis projection, depth not specified
+            # -> set 'large enough' depth using half the box diagonal + margin
+            depth = np.linalg.norm(ds.domain_width.in_units("code_length")) * 1.0001
+        depth = ds.coordinates.sanitize_depth(depth)[0]
+
         # center_rot normalizes the center to (0,0),
         # units match bounds
         # for SPH data, we want to input the original center
@@ -2479,7 +2483,6 @@ def __init__(
             width,
             ds,
             depth=depth,
-            get3bounds=True,
         )
         # will probably fail if you try to project an SPH and non-SPH
         # field in a single call

yt/visualization/tests/test_offaxisprojection_pytestonly.py

@@ -1,14 +1,18 @@
 import numpy as np
 import pytest
 import unyt
+from numpy.testing import assert_allclose
 
 from yt.testing import (
     assert_rel_equal,
     cubicspline_python,
+    fake_amr_ds,
     fake_sph_flexible_grid_ds,
+    fake_sph_grid_ds,
     integrate_kernel,
+    requires_module_pytest,
 )
-from yt.visualization.api import ProjectionPlot
+from yt.visualization.api import OffAxisProjectionPlot, ProjectionPlot
 
 
 @pytest.mark.parametrize("weighted", [True, False])
@@ -165,3 +169,77 @@ def makemasses(i, j, k):
     # print("expected:\n", expected_out)
     # print("recovered:\n", img.v)
     assert_rel_equal(expected_out, img.v, 4)
+
+
+@pytest.fixture
+def fame_amr_ds_for_proj():
+    return fake_amr_ds()
+
+
+_diag_dist = np.sqrt(3.0)  # diagonal distance of a cube with length 1.
+# each case is depth, center, expected integrated distance
+_cases_to_test = [
+    (_diag_dist / 3.0, "domain_left_edge", _diag_dist / 3.0 / 2.0),
+    (_diag_dist * 2.0, "domain_left_edge", _diag_dist),
+    (_diag_dist * 4.0, "domain_left_edge", _diag_dist),
+    (None, "domain_center", _diag_dist),
+]
+
+
+@pytest.mark.parametrize("depth,proj_center,expected", _cases_to_test)
+def test_offaxisprojection_depth(depth, proj_center, expected, fame_amr_ds_for_proj):
+    # this checks that the depth keyword argument works as expected.
+    # in all cases, it integrates the (index, ones) field for a normal
+    # pointing to the right edge corner of the domain.
+    #
+    # For the tests where the projection is centered on the left edge,
+    # the integrate distance will scale as depth / 2.0. When centered
+    # on the origin, it will scale with depth. The integrated distance
+    # should max out at the diagonal distance of the cube (when the depth
+    # exceeds the cube diagonal distance).
+    #
+    # Also note that the accuracy will depend on the buffer dimensions:
+    # using the default (800,800) results in accuracy of about 1 percent
+
+    ds = fame_amr_ds_for_proj
+
+    n = [1.0, 1.0, 1.0]
+    c = getattr(ds, proj_center)
+    field = ("index", "ones")
+
+    p = ProjectionPlot(ds, n, field, depth=depth, weight_field=None, center=c)
+
+    maxval = p.frb[field].max().d
+    assert_allclose(expected, maxval, atol=1e-2)
+
+
+_sph_test_cases = [
+    ([1.0, 1.0, 0.7], 27),
+    ([1.0, 1.0, 1], 19),
+    ([0.0, 0.0, 1], 9),
+]
+
+
+@requires_module_pytest("contourpy")
+@pytest.mark.parametrize("normal_vec, n_particles", _sph_test_cases)
+def test_offaxisprojection_sph_defaultdepth(normal_vec, n_particles):
+    # checks that particles are picked up as expected for a range of
+    # depths and normal vectors. Certain viewing angles will result
+    # in overlapping particles (since fake_sph_grid_ds aligns particles
+    # on a grid): n_particles is the expected number of circles in the
+    # resulting image for the given normal vector. Circle counts are
+    # calculated here using a contour generator.
+    from contourpy import contour_generator
+
+    ds = fake_sph_grid_ds()
+    c = ds.domain_center
+    diag_dist = np.linalg.norm(ds.domain_width)
+    field = ("gas", "mass")
+    p = OffAxisProjectionPlot(
+        ds, normal_vec, field, weight_field=None, center=c, width=diag_dist
+    )
+    p.render()
+
+    # get the number of circles in the plot
+    cg = contour_generator(z=p.frb[("gas", "mass")].d)
+    assert n_particles == len(cg.lines(1.0))