search

yt-project / yt

Main yt repository
http://yt-project.org
Other
468 stars 279 forks source link

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

Open evanmayer opened 9 months ago

evanmayer commented 9 months ago

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: ```python 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..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... : Sources: source_00: :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: : 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_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

yt is installed from source.

welcome[bot] commented 9 months ago

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.