BUG: wrong pixelization for "r-normal" projections in spherical coordinates ?

[neutrinoceros]

Bug report

Bug summary
originally reported as #3529

Code for reproduction

import yt

ds = yt.load_sample("bw_spherical_2d", unit_system="code")
p = yt.SlicePlot(ds, "r", "density")
p.save(f"/tmp/test.png")

Actual outcome

Expected outcome
Seemingly, working examples of this usage were documented in #1791
but it also looks like the important code here has really not changed much since then (except for formatting), and the PR included tests, so I'm surprised to see it failing now. I'll try to inspect this deeper later.

[neutrinoceros]

Here's another example that doesn't require yt.load_sample or the AMRVAC frontend, that I was able to build from @matthewturk's notebook embedded in #1791 (thank you so much for that Matt !)

It should be pretty clear here that the theta coordinate should go up to 2*PI

import numpy as np
import yt

shape = (72,  721, 1152)
prng = np.random.RandomState(0)
data = {("gas", "density"): prng.random_sample(shape)}
bbox = np.array([[1, 72], [0, 2.0*np.pi], [0, np.pi] ])

ds = yt.load_uniform_grid(
    data,
    shape,
    1.0,
    geometry=("spherical", ('r', 'theta', 'phi')),
    bbox=bbox
)

p = yt.SlicePlot(ds, "r", "density")
p.save(f"/tmp/test.png")

[neutrinoceros]

Using this patch

diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx
index 72803050c..ab7e896c7 100644
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -687,10 +687,10 @@ def pixelize_aitoff(np.float64_t[:] theta,
                 phi0 = math.asin(z*y)
                 # Now we just need to figure out which pixel contributes.
                 # We do not have a fast search.
-                if not (theta_p - dtheta_p <= theta0 <= theta_p + dtheta_p):
-                    continue
-                if not (phi_p - dphi_p <= phi0 <= phi_p + dphi_p):
-                    continue
+                #if not (theta_p - dtheta_p <= theta0 <= theta_p + dtheta_p):
+                #    continue
+                #if not (phi_p - dphi_p <= phi0 <= phi_p + dphi_p):
+                #    continue
                 img[i, j] = field[fi]
     return img

recompiling and running the example above gives

so it seems clear that the pixels are correctly selected, but apparently half of them are not correctly initialised

[neutrinoceros]

My example above is broken on the user side... theta is the angle that varies over a PI-wide interval, while phi varies over 2PI. The code dealing with this in yt is known to be inconsistent in naming these (#1796), but it should be ok on the user side. Here's a corrected version of my test

import numpy as np
import yt

shape = (72,  721, 1152)
prng = np.random.RandomState(0)
data = {("gas", "density"): prng.random_sample(shape)}
bbox = np.array([[1, 72], [0, np.pi], [0, 2*np.pi]])

ds = yt.load_uniform_grid(
    data,
    shape,
    1.0,
    geometry=("spherical", ('r', 'theta', 'phi')),
    bbox=bbox
)

p = yt.SlicePlot(ds, "r", "density")
p.save(f"/tmp/test1.png")


# transposition
bbox = np.array([[1, 72], [0, 2*np.pi], [0, np.pi]])

ds = yt.load_uniform_grid(
    data,
    shape,
    1.0,
    geometry=("spherical", ('r', 'phi', 'theta')),
    bbox=bbox
)

p.save(f"/tmp/test2.png")

both example produce this image

which I think is incorrect ? I'd except x/y axes to be flipped, now let's try fix that.

[matthewturk]

A few minor things -- before you get too far in, can you check what happens if you use the cartopy projections?

[neutrinoceros]

I'd love to, how do I do that ?

[neutrinoceros]

(I might be unable to actually, cartopy is notoriously hard to install outside of conda land)

[neutrinoceros]

yet another way this is broken:

import yt

ds = yt.load_sample("KeplerianDisk")
p = yt.SlicePlot(ds, "r", "density")
p.save(f"/tmp/disk_slice_r.png")

unroll !

likely a crucial difference with the previous example is that this data set has

>>> ds.coordinates.axis_order
('r', 'theta', 'phi')

while the previous one had

>>> ds.coordinates.axis_order
('r', 'phi', 'theta')

[neutrinoceros]

I think I'll need to adress #1796 thoroughly before I can fix this one