depth keyword in off-axis projections with sph particles has no effect?

[chrishavlin]

The depth keyword in OffAxisProjectionPlots has no effect with sph particles. Reading through the code it looks to me like while the depth will get used to set the bounds, when projecting with sph particles only the rotated x-y bounds and particle positions get used, so I'd guess that it has not been implemented fully but it'd be great to get confirmation from others who are more familiar with this part of the code before I attempt a fix.

reproduction

The following returns identical images:

import yt 
import numpy as np

ds = yt.load_sample("gizmo_64")  

z_hat_original = ds.arr([0.3, 0.3, 0.3], 'kpc')
y_hat_original = ds.arr([1., 0., 0.], 'kpc')

for d in np.arange(0, 10, 1):
    prj = yt.ProjectionPlot(ds, - z_hat_original, 
                                       center = ds.domain_center, 
                                       fields=('gas', 'density'),                                 
                                       depth = (d, 'kpc'), 
                                       north_vector = - y_hat_original, 
                                       )
    prj.save(f"depth_{str(d).zfill(2)}_proj.png")

I'd expect the images to vary as the depth parameter is changed.

what to do

If I'm not missing something, seems to me there are 2 options:

• short fix: raise a NotImplemented error when depth is used with particles
• better fix: implement it

code notes

By the time we get down to off_axis_projection, the bounds are still correct (accounting for the depth argument)

https://github.com/yt-project/yt/blob/9bd319d24ec3312ba8670efdc0cd5d5165f9d394/yt/visualization/volume_rendering/off_axis_projection.py#L211-L229

But when calculating the rotated coordinates and bounds, only the rotated x-y are returned and used in the kernel projection

https://github.com/yt-project/yt/blob/9bd319d24ec3312ba8670efdc0cd5d5165f9d394/yt/utilities/lib/pixelization_routines.pyx#L1934-L1950

So I think we'd want rotate_particle_coord to also return the rotated z vals and the rotated z bounds and then have pixelize_sph_kernel_projection also discard particles outsize those rotated z bounds.. I'm happy to attempt that, but wanted to double check that I'm not going down the wrong path here...

[chrishavlin]

note: thanks to @y-samuel-lu for raising the issue on slack!

[matthewturk]

I believe you're right. My recollection is that for on-axis projections, axis aligned bounding boxes (AABB) are used, but this was not implemented here. I think that it would be the appropriate thing to do to use the transformed z coordinates.

[nastasha-w]

Hi, I'm pretty new to YT, but I have been looking into SPH data projections for a different issue (#4788) that involves the pixelize_sph_kernel_projection function. The approach sounds reasonable enough to me. I'd suggest doing a simple cut on SPH particle positions in the z direction; the kernel interpolation table won't work for integrating through only partway through the smoothing-length-diameter sphere, and I figure this is one of those things that's good enough.

One issue I do want to point out involves periodic datasets. Now, I think we can reasonably get away with ignoring cases (or raising errors) if someone analyzing a periodic volume wants to make a slice through the whole volume that wraps back around, where the slice width is larger than the box size (e.g. ASCII sketch below).

+--------------------+
|       \.     \     |
|.       \.     \    |
|.        \.     \   |
|.   \     \         |
|.    \     \        |
+--------------------+

However, I can imagine a case where someone analyzing e.g., a cosmological volume wants to make an image of a galaxy overlapping the periodic boundary, aligned with the galaxy edge-on/face-on direction. To handle that, I'd suggest re-centering the output coordinates on zero, or half the box size, then applying any periodicity adjustments in the output arrays, before rotating them. (It would probably make the most sense to do that in the same loop as the rotation.) We can then output the rotated bounds as (0.5 * boxsize -0.5 * width, 0.5 * boxsize + 0.5* width, etc.) to match the adjusted coordinates. (This is possible in numpy using float modulo division, but I don't know how efficient that is.)

[chrishavlin]

Thanks for the suggestions! I need to look back at where my PR stalled...

One issue I do want to point out involves periodic datasets.

If I remember right, the current state does not handle periodcitiy at all for off-axis projections and particles but I may be remembering incorrectly. In any case, yes, it'd be good to properly handle periodicity.