Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
chrishavlin opened this issue Oct 24, 2023 · 4 comments · May be fixed by #4712
Open

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

chrishavlin opened this issue Oct 24, 2023 · 4 comments · May be fixed by #4712

Comments

@chrishavlin
Copy link
Contributor

chrishavlin commented Oct 24, 2023

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)

bounds = [x_min, x_max, y_min, y_max, z_min, z_max]
if weight is None:
for chunk in data_source.chunks([], "io"):
off_axis_projection_SPH(
chunk[ptype, ppos[0]].to("code_length").d,
chunk[ptype, ppos[1]].to("code_length").d,
chunk[ptype, ppos[2]].to("code_length").d,
chunk[ptype, "mass"].to("code_mass").d,
chunk[ptype, "density"].to("code_density").d,
chunk[ptype, "smoothing_length"].to("code_length").d,
bounds,
center.to("code_length").d,
width.to("code_length").d,
chunk[item].in_units(ounits),
buf,
mask,
normal_vector,
north,

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

px_rotated, py_rotated, \
rot_bounds_x0, rot_bounds_x1, \
rot_bounds_y0, rot_bounds_y1 = rotate_particle_coord(px, py, pz,
center, width, normal_vector, north_vector)
pixelize_sph_kernel_projection(projection_array,
mask,
px_rotated,
py_rotated,
smoothing_lengths,
particle_masses,
particle_densities,
quantity_to_smooth,
[rot_bounds_x0, rot_bounds_x1,
rot_bounds_y0, rot_bounds_y1],
weight_field=weight_field,
check_period=0)

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
Copy link
Contributor Author

chrishavlin commented Oct 24, 2023

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

@matthewturk
Copy link
Member

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.

@chrishavlin chrishavlin linked a pull request Oct 25, 2023 that will close this issue
5 tasks
@nastasha-w
Copy link

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
Copy link
Contributor Author

chrishavlin commented Jan 29, 2024

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants