Error while rendering SpectralCube with non-square face: IndexError while reading fields

[evanmayer]

Bug report

Bug summary

When calling sc.show() on a scene composed of a SpectralCube object with a non-square face, the call fails with an IndexError.

Code for reproduction

from astropy.io import fits
from spectral_cube import SpectralCube # pip install spectral-cube
import yt

fits_data = fits.open('./evanmayer_mwe.fits')
cube = SpectralCube.read(fits_data)
fits_data.close()

# https://spectral-cube.readthedocs.io/en/latest/yt_example.html
ytcube = cube.to_yt()
sc = yt.create_scene(ytcube.dataset, ('fits','flux'))
sc.show()

Unfortunately, I think this requires the file I'm using to reproduce (attached). I'm also not sure if the problem lies in yt or SpectralCube, since the yt object is created by the SpectralCube to_yt() method.

FITS file

I failed to reproduce this with just a normal data volume with a non-square face:

import yt
import numpy as np

N = 8
data = {"density": np.random.random((2*N, N, N))}

ds = yt.load_uniform_grid(
    data,
    [2*N, N, N],
    bbox=np.array([[0.0, 2.0], [0.0, 1.0], [0.0, 1.0]]),
)
sc = yt.create_scene(ds)
sc.show()

Maybe I didn't try hard enough, or this kind of MWE doesn't exercise the failing code?

Actual outcome

Error message:

yt : [WARNING  ] 2024-02-03 11:08:56,206 Cannot find time
yt : [INFO     ] 2024-02-03 11:08:56,209 Detected these axes: RA---SIN DEC--SIN FREQ 
yt : [WARNING  ] 2024-02-03 11:08:56,226 No length conversion provided. Assuming 1 = 1 cm.
yt : [WARNING  ] 2024-02-03 11:08:56,229 Assuming 1.0 = 1.0 s
yt : [WARNING  ] 2024-02-03 11:08:56,232 Assuming 1.0 = 1.0 g
yt : [INFO     ] 2024-02-03 11:08:56,288 Parameters: current_time              = 0.0
yt : [INFO     ] 2024-02-03 11:08:56,289 Parameters: domain_dimensions         = [120  24 188]
yt : [INFO     ] 2024-02-03 11:08:56,290 Parameters: domain_left_edge          = [0.5 0.5 0.5]
yt : [INFO     ] 2024-02-03 11:08:56,291 Parameters: domain_right_edge         = [120.5  24.5 188.5]
yt : [INFO     ] 2024-02-03 11:08:56,292 Parameters: cosmological_simulation   = 0
yt : [INFO     ] 2024-02-03 11:08:56,302 Adding field flux to the list of fields.
yt : [INFO     ] 2024-02-03 11:08:56,386 Rendering scene (Can take a while).
yt : [INFO     ] 2024-02-03 11:08:56,388 Creating volume

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
File ~/.local/lib/python3.10/site-packages/IPython/core/formatters.py:344, in BaseFormatter.__call__(self, obj)
    342     method = get_real_method(obj, self.print_method)
    343     if method is not None:
--> 344         return method()
    345     return None
    346 else:

File ~/github/yt/yt/visualization/volume_rendering/scene.py:987, in Scene._repr_png_(self)
    985 def _repr_png_(self):
    986     if self._last_render is None:
--> 987         self.render()
    988     png = self._last_render.write_png(
    989         filename=None, sigma_clip=self._sigma_clip, background="black"
    990     )
    991     self._sigma_clip = None

File ~/github/yt/yt/visualization/volume_rendering/scene.py:225, in Scene.render(self, camera)
    223 assert camera is not None
    224 self._validate()
--> 225 bmp = self.composite(camera=camera)
    226 self._last_render = bmp
    227 return bmp

File ~/github/yt/yt/visualization/volume_rendering/scene.py:598, in Scene.composite(self, camera)
    595     im = source.zbuffer.rgba
    597 for _, source in self.transparent_sources:
--> 598     im = source.render(camera, zbuffer=opaque)
    599     opaque.rgba = im
    601 # rotate image 180 degrees so orientation agrees with e.g.
    602 # a PlotWindow plot

File ~/github/yt/yt/visualization/volume_rendering/render_source.py:131, in validate_volume.<locals>.wrapper(*args, **kwargs)
    129     log_fields.append(obj.log_field)
    130 if not obj._volume_valid:
--> 131     obj.volume.set_fields(
    132         fields, log_fields, no_ghost=(not obj.use_ghost_zones)
    133     )
    134 obj._volume_valid = True
    135 return f(*args, **kwargs)

File ~/github/yt/yt/utilities/amr_kdtree/amr_kdtree.py:232, in AMRKDTree.set_fields(self, fields, log_fields, no_ghost, force)
    229 self.brick_dimensions = []
    230 bricks = []
--> 232 for b in self.traverse():
    233     list(map(_apply_log, b.my_data, flip_log, self.log_fields))
    234     bricks.append(b)

File ~/github/yt/yt/utilities/amr_kdtree/amr_kdtree.py:250, in AMRKDTree.traverse(self, viewpoint)
    248 def traverse(self, viewpoint=None):
    249     for node in self.tree.trunk.kd_traverse(viewpoint=viewpoint):
--> 250         yield self.get_brick_data(node)

File ~/github/yt/yt/utilities/amr_kdtree/amr_kdtree.py:339, in AMRKDTree.get_brick_data(self, node)
    337 else:
    338     dds = []
--> 339     vcd = grid.get_vertex_centered_data(
    340         self.fields, smoothed=True, no_ghost=self.no_ghost
    341     )
    342     for i, field in enumerate(self.fields):
    343         if self.log_fields[i]:

File ~/github/yt/yt/data_objects/index_subobjects/grid_patch.py:289, in AMRGridPatch.get_vertex_centered_data(self, fields, smoothed, no_ghost)
    285 if no_ghost:
    286     for field in fields:
    287         # Ensure we have the native endianness in this array.  Avoid making
    288         # a copy if possible.
--> 289         old_field = np.asarray(self[field], dtype="=f8")
    290         # We'll use the ghost zone routine, which will naturally
    291         # extrapolate here.
    292         input_left = np.array([0.5, 0.5, 0.5], dtype="float64")

File ~/github/yt/yt/data_objects/index_subobjects/grid_patch.py:76, in AMRGridPatch.__getitem__(self, key)
     75 def __getitem__(self, key):
---> 76     tr = super().__getitem__(key)
     77     try:
     78         fields = self._determine_fields(key)

File ~/github/yt/yt/data_objects/data_containers.py:235, in YTDataContainer.__getitem__(self, key)
    233         return self.field_data[f]
    234     else:
--> 235         self.get_data(f)
    236 # fi.units is the unit expression string. We depend on the registry
    237 # hanging off the dataset to define this unit object.
    238 # Note that this is less succinct so that we can account for the case
    239 # when there are, for example, no elements in the object.
    240 try:

File ~/github/yt/yt/data_objects/selection_objects/data_selection_objects.py:205, in YTSelectionContainer.get_data(self, fields)
    201         fluids.append(field_key)
    202 # The _read method will figure out which fields it needs to get from
    203 # disk, and return a dict of those fields along with the fields that
    204 # need to be generated.
--> 205 read_fluids, gen_fluids = self.index._read_fluid_fields(
...
    413     for i, sl in enumerate(slices):
--> 414         dest[offset : offset + count, i] = source[tuple(sl)][np.squeeze(mask)]
    415 return count

IndexError: boolean index did not match indexed array along dimension 2; dimension is 3 but corresponding boolean dimension is 28

Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...
<Scene Object>:
Sources: 
    source_00: <Volume Source>:YTRegion (InMemoryFITSFile_47550824a5de43229d56af5a0b742a3d): , center=[60.5 12.5 94.5] cm, left_edge=[0.5 0.5 0.5] cm, right_edge=[120.5  24.5 188.5] cm transfer_function:None
Camera: 
    <Camera Object>:
	position:[120.5  24.5 188.5] code_length
	focus:[60.5 12.5 94.5] code_length
	north_vector:[-0.82180894 -0.16436179  0.54554126]
	width:[282. 282. 282.] code_length
	light:None
	resolution:(512, 512)
Lens: <Lens Object>:
	lens_type:plane-parallel
	viewpoint:[-1.50855129e+08 -3.01710254e+07 -2.36339703e+08] code_length
yt : [INFO     ] 2024-02-03 11:01:17,440 Parameters: current_time              = 0.0
yt : [INFO     ] 2024-02-03 11:01:17,442 Parameters: domain_dimensions         = [16  8  8]
yt : [INFO     ] 2024-02-03 11:01:17,444 Parameters: domain_left_edge          = [0. 0. 0.]
yt : [INFO     ] 2024-02-03 11:01:17,445 Parameters: domain_right_edge         = [2. 1. 1.]
yt : [INFO     ] 2024-02-03 11:01:17,446 Parameters: cosmological_simulation   = 0
yt : [INFO     ] 2024-02-03 11:01:17,583 Setting default field to ('gas', 'density')
yt : [INFO     ] 2024-02-03 11:01:17,589 Rendering scene (Can take a while).
yt : [INFO     ] 2024-02-03 11:01:17,602 Creating volume
yt : [INFO     ] 2024-02-03 11:01:17,624 Creating transfer function
yt : [INFO     ] 2024-02-03 11:01:17,625 Calculating data bounds. This may take a while. Set the TransferFunctionHelper.bounds to avoid this.

YTRegion (UniformGridData): , center=[1.  0.5 0.5] cm, left_edge=[0. 0. 0.] cm, right_edge=[2. 1. 1.] cm

Expected outcome

yt renders the cube, just as it does for cubes with square faces.

Version Information

• Operating System: Ubuntu 22.04.3 LTS
• Python Version: 3.10.13
• yt version: 4.3.dev0
• Other Libraries (if applicable): astropy

yt is installed from source.

[welcome]

Hi, and welcome to yt! Thanks for opening your first issue. We have an issue template that helps us to gather relevant information to help diagnosing and fixing the issue.