[BUG] Querying a reloaded sphere dataset fails when chunked

Hello, unfortunately, the following script fails when run on a reloaded sphere dataset from an Enzo simulation:

import yt

sp_ds = yt.load("sphere_dataset.h5")  # saved intermediate sphere
sp_ds.data["grid", "density"]

The following issue is raised:

TypeError: Boolean indexing array has incompatible shape

The culprit line is at https://github.com/yt-project/yt/blob/main/yt/frontends/ytdata/io.py#L233-L235 and follows from the else branch a few lines above. It was introduced in 6405200.

In the code a selector object is used to return what I am assuming is a boolean mask. However, in this instance, the branch sizes the input array to selector from a slice of the data file (a chunk subset). In my case the reloaded sphere is large and split into 173 chunks (according to the index attrs). The mask is therefore sized to these chunks. The file which is loaded with h5py is the full data file, meaning when this mask is applied to f in L234, we get the boolean indexing shape error. I presume this was an oversight in how index_mask was implemented for chunked datasets.

By contrast, if I do:

sp_ds.all_data()["grid", "density"]

No error is encountered.

Versions:
yt: 4.4.2
numpy: 2.5
h5py: 3.16
Python: 3.14
OS: Linux

I am not intimately familiar with the source code but would be happy to work on a PR for this if it is a relatively straightforward fix as I believe it should be. The sphere dataset I'm using is relatively large (~10GB). My cluster just went down so I can't access the full traceback but if you need that I'd be happy to grab it when it's back up.

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.

thanks for the report! Unfortunately the re-loading from saved objects has been a tad flaky lately. I don't think it would matter too much, but mind showing how you saved the sphere? the dataset type might have an influence on behavior. a full reproducible example would be helpful too if possible, e.g., along the lines of

from yt.testing import fake_amr_ds

ds = fake_amr_ds()
<--- save and reload sphere --_.

@chrishavlin

Hello @chrishavlin, thank you for your reply. Here is a reproducible example following the test dataset you specified:

import yt; from yt.testing import fake_amr_ds

ds = fake_amr_ds()
sp = ds.sphere(ds.domain_center, (1.0, "kpc"))
sp.save_as_dataset("test_sphere.h5", fields=ds.field_list)

sp_ds = yt.load("test_sphere.h5")
print(sp_ds.data[("grid", "Density")])

Produces:

Traceback (most recent call last):
  File "<python-input-7>", line 1, in <module>
    print(sp_ds.data[("grid", "Density")])
          ~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^
  File ".../python3.14/site-packages/yt/data_objects/data_containers.py", line 248, in __getitem__
    self.get_data(f)
    ~~~~~~~~~~~~~^^^
  File ".../python3.14/site-packages/yt/data_objects/selection_objects/data_selection_objects.py", line 211, in get_data
    read_particles, gen_particles = self.index._read_particle_fields(
                                    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
        particles, self, self._current_chunk
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File ".../python3.14/site-packages/yt/geometry/geometry_handler.py", line 227, in _read_particle_fields
    fields_to_return = self.io._read_particle_selection(
        chunks, selector, fields_to_read
    )
  File ".../python3.14/site-packages/yt/utilities/io_handler.py", line 194, in _read_particle_selection
    for field_r, vals in self._read_particle_fields(chunks, ptf, selector):
                         ~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^
  File ".../python3.14/site-packages/yt/utilities/io_handler.py", line 222, in _read_particle_fields
    data_file_data = self._read_particle_data_file(data_file, ptf, selector)
  File ".../python3.14/site-packages/yt/frontends/ytdata/io.py", line 234, in _read_particle_data_file
    data = f[ptype][field][mask].astype("float64", copy=False)
           ~~~~~~~~~~~~~~~^^^^^^
  File "h5py/_objects.pyx", line 54, in h5py._objects.with_phil.wrapper
  File "h5py/_objects.pyx", line 55, in h5py._objects.with_phil.wrapper
  File ".../python3.14/site-packages/h5py/_hl/dataset.py", line 945, in __getitem__
    selection = sel.select(self.shape, args, dataset=self)
  File ".../python3.14/site-packages/h5py/_hl/selections.py", line 67, in select
    raise TypeError("Boolean indexing array has incompatible shape")
TypeError: Boolean indexing array has incompatible shape

My apologies for not including this originally.

Original script from my code

def create_star_sphere(
    row: dict[str, object], 
    outpath: Path, 
    radius_kpc: float, 
) -> Path:
    """
    Generates a sphere of specified radius_kpc to outpath, and stores the fields
    specified in the module-level SPHERE_FIELDS. This function should be called by workers.
    """
    yt.set_log_level("error")

    ds = yt.load(row["snapshot_path"])
    add_metallicity3(ds=ds)

    centre = ds.arr(row["positions_unitary"], "unitary")
    radius = ds.quan(radius_kpc, "kpc")

    sp = ds.sphere(centre, radius)
    sp.get_data(SPHERE_FIELDS)  # faster to use get_data so all queries are done in one pass then stored on the sphere
    sp.save_as_dataset(filename=str(outpath), fields=SPHERE_FIELDS)  # and then save

    del ds
    gc.collect()

    return outpath

thanks!

@jp-duminy

ok, thanks for the reproducer, @jp-duminy ! it was helpful. and i think you're right -- it was an oversight in how chunking is applied to the saved file (and a lack of tests).

The index_mask in that section of code should select the index range from the full file that is being handled in the current chunk. So in the else branch that you pointed out, the mask will be the size of the index range given in index_mask (smaller than the total range of the file) and the corresponding selection would need to apply that mask after the index selection. e.g., I think line 234 would become

data = f[ptype][field][index_mask][mask].astype("float64", copy=False)

where mask applies just to the index_mask subselection.

Would need to clean up the other branch of the if too so that it will select everything from the index range mask returns everything from that subset (probably just mask=slice(None)).

A PR would be welcome if you want! We're actually trying to get yt 4.4.3 out the door very soon, so it'd be good timing :)

Splendid, I will try to get this done for you at my soonest convenience.

Side note: with a cursory ctrl+f on the file for slices, I see this class https://github.com/yt-project/yt/blob/main/yt/frontends/ytdata/io.py#L264C1-L301C51 is doing a similar operation on the z-axis with the all_count[ptype] on the _read_particle_fields() function. It follows the same logic of constructing a boolean mask on a (potentially chunked) file which is then applied to the h5py file. The coords function looks fine though since the z-axis is sized to x directly. This bug may therefore be present elsewhere and require a more careful treatment: I see the chunking lines went in ~3 years ago.

on a quick read it looks to me like the other spots you highlighted are not quite the same issue: looks like they are reading and masking the whole file on every iteration rather than masking just the index range for the chunk. very likely it could have the same performance issue that the index_mask here was aiming to solve but probably wouldn't surface any errors. and sometimes the data file objects really do reference different hdf files (chunking can be index ranges within a file or chunks across files). certainly worth a closer look but I'd split that work out separately since I think it'd take signicantly more time to sort out.

Navigation Menu

Repository navigation

[BUG] Querying a reloaded sphere dataset fails when chunked #5576

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

[BUG] Querying a reloaded sphere dataset fails when chunked #5576

Description

Activity

welcome commented on Oct 8, 2026

chrishavlin commented on Oct 8, 2026

jp-duminy commented on Oct 8, 2026

chrishavlin commented on Oct 8, 2026

chrishavlin commented on Oct 8, 2026

jp-duminy commented on Oct 8, 2026

chrishavlin commented on Oct 8, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions