BUG: Broken sample outputs of off-axis projection plots in docs

[xshaokun]

Bug report

Bug summary

The outputs of two samples are broken in the docs
https://yt-project.org/docs/dev/cookbook/simple_plots.html#off-axis-projection
https://yt-project.org/docs/dev/visualizing/plots.html#off-axis-projections

The same code produces the following outputs when run with yt-4.3.1

import yt

ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
L = [1, 1, 0]  # vector normal to cutting plane
north_vector = [-1, 1, 0]
prj = yt.ProjectionPlot(
    ds, L, ("gas", "density"), width=(25, "kpc"), north_vector=north_vector
)
prj.save()

import yt

# Load the dataset.
ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")

# Create a 15 kpc radius sphere, centered on the center of the sim volume
sp = ds.sphere("center", (15.0, "kpc"))

# Get the angular momentum vector for the sphere.
L = sp.quantities.angular_momentum_vector()

print(f"Angular momentum vector: {L}")

# Create an off-axis ProjectionPlot of density centered on the object with the L
# vector as its normal and a width of 25 kpc on a side
p = yt.ProjectionPlot(
    ds, L, fields=("gas", "density"), center=sp.center, width=(25, "kpc")
)
p.save()

[xshaokun]

The color map is different because I have locally customized it in yt.toml

[neutrinoceros]

Thanks for reporting. I don't have time to handle it myself now but I trust you are correct in that is a regression in yt 4.4.0 so it should be addressed in 4.4.1.

[chrishavlin]

Not sure of the source yet, but I did narrow the regression to the big sph backend changes (#4939).

Checking out the last commit (hash 77714f7) just before the sph backend changes, you get the correct expected behavior. Checking out the sph backend changes (hash cadddb6), you get the incorrect result. Shouldn't be too hard to figure out the issue, since there shouldn't be too much overlap in parts of the code that were changed that would affect both grid-based and sph-based off axis projections... I'll keep looking for a bit, but @nastasha-w if you're available, mind taking a look as well in case you find it quicker?

Also, i'm a bit baffled that we don't have an answer test to catch this degree of a regression...

[chrishavlin]

well i've narrowed it a bit more -- has to do with the updated _get_oblique_window_parameters and default depth.

yt/yt/visualization/plot_window.py

Lines 85 to 112 in e6e5f1f

	def get_oblique_window_parameters(
	normal, center, width, ds, depth=None, get3bounds=False
	):
	center, display_center = ds.coordinates.sanitize_center(center, axis=None)
	width = ds.coordinates.sanitize_width(normal, width, depth)
	if len(width) == 2:
	# Transforming to the cutting plane coordinate system
	# the original dimensionless center messes up off-axis
	# SPH projections though -> don't use this center there
	center = (
	(center - ds.domain_left_edge) / ds.domain_width - 0.5
	) * ds.domain_width
	(normal, perp1, perp2) = ortho_find(normal)
	mat = np.transpose(np.column_stack((perp1, perp2, normal)))
	center = np.dot(mat, center)
	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)

If I apply the following diff that adds back the original _get_oblique_window_parameters and changes the default depth back to (1,""), I get the expected result (not saying that it is the fix, but it is the source of the issue):

diff --git a/yt/visualization/plot_window.py b/yt/visualization/plot_window.py
index 4a52e803f..ee1e890d5 100644
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -83,8 +83,23 @@ def get_window_parameters(axis, center, width, ds):
     )
     return (bounds, center, display_center)
 
+def get_oblique_window_parameters(normal, center, width, ds, depth=None):
+    center, display_center = ds.coordinates.sanitize_center(center, axis=None)
+    width = ds.coordinates.sanitize_width(normal, width, depth)
 
-def get_oblique_window_parameters(
+    if len(width) == 2:
+        # Transforming to the cutting plane coordinate system
+        center = (center - ds.domain_left_edge) / ds.domain_width - 0.5
+        (normal, perp1, perp2) = ortho_find(normal)
+        mat = np.transpose(np.column_stack((perp1, perp2, normal)))
+        center = np.dot(mat, center)
+
+    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))
+
+    return (bounds, center)
+            
+def get_oblique_window_parameters_sph(
     normal, center, width, ds, depth=None, get3bounds=False
 ):
     center, display_center = ds.coordinates.sanitize_center(center, axis=None)
@@ -2444,7 +2459,7 @@ class OffAxisProjectionPlot(ProjectionPlot, PWViewerMPL):
         fields,
         center="center",
         width=None,
-        depth=None,
+        depth=(1, "1"),
         axes_unit=None,
         weight_field=None,
         max_level=None,
@@ -2474,14 +2489,7 @@ class OffAxisProjectionPlot(ProjectionPlot, PWViewerMPL):
         # rotation.
         # get3bounds gets a depth 0.5 * diagonal + margin in the
         # depth=None case.
-        (bounds, center_rot) = get_oblique_window_parameters(
-            normal,
-            center,
-            width,
-            ds,
-            depth=depth,
-            get3bounds=True,
-        )
+
         # will probably fail if you try to project an SPH and non-SPH
         # field in a single call
         # checks for SPH fields copied from the
@@ -2497,8 +2505,23 @@ class OffAxisProjectionPlot(ProjectionPlot, PWViewerMPL):
         particle_datasets = (ParticleDataset, StreamParticlesDataset)
 
         if isinstance(data_source.ds, particle_datasets) and is_sph_field:
+            (bounds, center_rot) = get_oblique_window_parameters_sph(
+            normal,
+            center,
+            width,
+            ds,
+            depth=depth,
+            get3bounds=True,
+            )
             center_use = parse_center_array(center, ds=data_source.ds, axis=None)
         else:
+            (bounds, center_rot) = get_oblique_window_parameters(
+            normal,
+            center,
+            width,
+            ds,
+            depth=depth,            
+            )
             center_use = center_rot
         fields = list(iter_fields(fields))[:]
         # oap_width = ds.arr(

[chrishavlin]

aaand one more update: has to do mainly with how the depth is handled: supplying a depth works as expected. e.g., the following returns the correct image

import yt
ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
L = [1, 1, 0]  # vector normal to cutting plane
north_vector = [-1, 1, 0]
prj = yt.ProjectionPlot(
    ds, L, ("gas", "density"),  width=(25, 'kpc'), north_vector=north_vector, depth = (1, 'code_length')
)
prj.show()

[chrishavlin]

fix incoming...

[xshaokun]

Update another broken case in docs, not sure if it is caused by the same source, but is also related to the method ProjectionPlot
https://yt-project.org/docs/dev/analyzing/particle_filter.html#Filtering-Particle-Data

[neutrinoceros]

I think this is #3980, but worth inspecting if it's solved by Chris' patch too

[nastasha-w]

@xshaokun thanks for finding this! I'm sorry my SPH changes broke this.

[chrishavlin]

confirming that my fix won't fix the particle filtering issue -- definitely a separate issue.