This example creates a fake in-memory particle dataset and then loads it as a yt dataset using the load_particles function.
Our “fake” dataset will be numpy arrays filled with normally distributed randoml particle positions and uniform particle masses. Since real data is often scaled, I arbitrarily multiply by 1e6 to show how to deal with scaled data.
[1]:
import numpy as np
n_particles = 5000000
ppx, ppy, ppz = 1e6 * np.random.normal(size=[3, n_particles])
ppm = np.ones(n_particles)
The load_particles function accepts a dictionary populated with particle data fields loaded in memory as numpy arrays or python lists:
[2]:
data = {
"particle_position_x": ppx,
"particle_position_y": ppy,
"particle_position_z": ppz,
"particle_mass": ppm,
}
To hook up with yt’s internal field system, the dictionary keys must be ‘particle_position_x’, ‘particle_position_y’, ‘particle_position_z’, and ‘particle_mass’, as well as any other particle field provided by one of the particle frontends.
The load_particles function transforms the data dictionary into an in-memory yt Dataset object, providing an interface for further analysis with yt. The example below illustrates how to load the data dictionary we created above.
[3]:
import yt
from yt.units import Msun, parsec
bbox = 1.1 * np.array(
[[min(ppx), max(ppx)], [min(ppy), max(ppy)], [min(ppz), max(ppz)]]
)
ds = yt.load_particles(data, length_unit=1.0 * parsec, mass_unit=1e8 * Msun, bbox=bbox)
The length_unit and mass_unit are the conversion from the units used in the data dictionary to CGS. I’ve arbitrarily chosen one parsec and 10^8 Msun for this example.
The n_ref parameter controls how many particle it takes to accumulate in an oct-tree cell to trigger refinement. Larger n_ref will decrease poisson noise at the cost of resolution in the octree.
Finally, the bbox parameter is a bounding box in the units of the dataset that contains all of the particles. This is used to set the size of the base octree block.
This new dataset acts like any other yt Dataset object, and can be used to create data objects and query for yt fields. This example shows how to access “deposit” fields:
[4]:
ad = ds.all_data()
# This is generated with "cloud-in-cell" interpolation.
cic_density = ad["deposit", "all_cic"]
# These three are based on nearest-neighbor cell deposition
nn_density = ad["deposit", "all_density"]
nn_deposited_mass = ad["deposit", "all_mass"]
particle_count_per_cell = ad["deposit", "all_count"]
---------------------------------------------------------------------------
NeedsGridType Traceback (most recent call last)
File /tmp/yt/yt/data_objects/data_containers.py:287, in YTDataContainer._generate_fluid_field(self, field)
286 try:
--> 287 finfo.check_available(gen_obj)
288 except NeedsGridType as ngt_exception:
File /tmp/yt/yt/fields/derived_field.py:235, in DerivedField.check_available(self, data)
234 for validator in self.validators:
--> 235 validator(data)
236 # If we don't get an exception, we're good to go
File /tmp/yt/yt/fields/derived_field.py:620, in ValidateSpatial.__call__(self, data)
619 if not getattr(data, "_spatial", False):
--> 620 raise NeedsGridType(self.ghost_zones, self.fields)
621 if self.ghost_zones <= data._num_ghost_zones:
NeedsGridType: fields None require 0 ghost zones.
During handling of the above exception, another exception occurred:
YTNonIndexedDataContainer Traceback (most recent call last)
Cell In[4], line 4
1 ad = ds.all_data()
3 # This is generated with "cloud-in-cell" interpolation.
----> 4 cic_density = ad["deposit", "all_cic"]
6 # These three are based on nearest-neighbor cell deposition
7 nn_density = ad["deposit", "all_density"]
File /tmp/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 /tmp/yt/yt/data_objects/selection_objects/data_selection_objects.py:220, in YTSelectionContainer.get_data(self, fields)
217 self.field_data[f].convert_to_units(finfos[f].output_units)
219 fields_to_generate += gen_fluids + gen_particles
--> 220 self._generate_fields(fields_to_generate)
221 for field in list(self.field_data.keys()):
222 if field not in ofields:
File /tmp/yt/yt/data_objects/selection_objects/data_selection_objects.py:252, in YTSelectionContainer._generate_fields(self, fields_to_generate)
250 fi = self.ds._get_field_info(field)
251 try:
--> 252 fd = self._generate_field(field)
253 if hasattr(fd, "units"):
254 fd.units.registry = self.ds.unit_registry
File /tmp/yt/yt/data_objects/data_containers.py:273, in YTDataContainer._generate_field(self, field)
271 tr = self._generate_particle_field(field)
272 else:
--> 273 tr = self._generate_fluid_field(field)
274 if tr is None:
275 raise YTCouldNotGenerateField(field, self.ds)
File /tmp/yt/yt/data_objects/data_containers.py:289, in YTDataContainer._generate_fluid_field(self, field)
287 finfo.check_available(gen_obj)
288 except NeedsGridType as ngt_exception:
--> 289 rv = self._generate_spatial_fluid(field, ngt_exception.ghost_zones)
290 else:
291 rv = finfo(gen_obj)
File /tmp/yt/yt/data_objects/data_containers.py:315, in YTDataContainer._generate_spatial_fluid(self, field, ngz)
313 o = self._current_chunk.objs[0]
314 if accumulate:
--> 315 rv = self.ds.arr(np.empty(o.ires.size, dtype="float64"), units)
316 outputs.append(rv)
317 ind = 0 # Does this work with mesh?
File /tmp/yt/yt/data_objects/selection_objects/data_selection_objects.py:451, in YTSelectionContainer.ires(self)
449 if self._current_chunk is None:
450 self.index._identify_base_chunk(self)
--> 451 return self._current_chunk.ires
File /tmp/yt/yt/geometry/geometry_handler.py:272, in cacheable_property.<locals>.cacheable_func(self)
270 return getattr(self, n)
271 if self.data_size is None:
--> 272 tr = self._accumulate_values(n[1:])
273 else:
274 tr = func(self)
File /tmp/yt/yt/geometry/geometry_handler.py:308, in YTDataChunk._accumulate_values(self, method)
306 for obj in self._fast_index or self.objs:
307 f = getattr(obj, mname)
--> 308 arrs.append(f(self.dobj))
309 if method == "dtcoords":
310 arrs = [arr[0] for arr in arrs]
File /tmp/yt/yt/data_objects/index_subobjects/particle_container.py:17, in _non_indexed.<locals>._func_non_indexed(self, *args, **kwargs)
16 def _func_non_indexed(self, *args, **kwargs):
---> 17 raise YTNonIndexedDataContainer(self)
YTNonIndexedDataContainer: The data container type (ParticleContainer) is an unindexed type. Operations such as ires, icoords, fcoords and fwidth will not work on it.
Did you just attempt to perform an off-axis operation ? Be sure to consult the latest documentation to see whether the operation you tried is actually supported for your data type.
[5]:
ds.field_list
[5]:
[('all', 'particle_mass'),
('all', 'particle_position_x'),
('all', 'particle_position_y'),
('all', 'particle_position_z'),
('io', 'particle_mass'),
('io', 'particle_position_x'),
('io', 'particle_position_y'),
('io', 'particle_position_z'),
('nbody', 'particle_mass'),
('nbody', 'particle_position_x'),
('nbody', 'particle_position_y'),
('nbody', 'particle_position_z')]
[6]:
ds.derived_field_list
[6]:
[('all', 'particle_mass'),
('all', 'particle_ones'),
('all', 'particle_position'),
('all', 'particle_position_cylindrical_radius'),
('all', 'particle_position_cylindrical_theta'),
('all', 'particle_position_cylindrical_z'),
('all', 'particle_position_relative_x'),
('all', 'particle_position_relative_y'),
('all', 'particle_position_relative_z'),
('all', 'particle_position_spherical_phi'),
('all', 'particle_position_spherical_radius'),
('all', 'particle_position_spherical_theta'),
('all', 'particle_position_x'),
('all', 'particle_position_y'),
('all', 'particle_position_z'),
('all', 'particle_radius'),
('all', 'relative_particle_position'),
('all', 'relative_particle_position_x'),
('all', 'relative_particle_position_y'),
('all', 'relative_particle_position_z'),
('deposit', 'all_cic'),
('deposit', 'all_count'),
('deposit', 'all_density'),
('deposit', 'all_mass'),
('deposit', 'io_cic'),
('deposit', 'io_count'),
('deposit', 'io_density'),
('deposit', 'io_mass'),
('deposit', 'nbody_cic'),
('deposit', 'nbody_count'),
('deposit', 'nbody_density'),
('deposit', 'nbody_mass'),
('gas', 'cell_volume'),
('gas', 'dx'),
('gas', 'dy'),
('gas', 'dz'),
('gas', 'path_element_x'),
('gas', 'path_element_y'),
('gas', 'path_element_z'),
('gas', 'volume'),
('gas', 'x'),
('gas', 'y'),
('gas', 'z'),
('index', 'cell_volume'),
('index', 'cylindrical_radius'),
('index', 'cylindrical_theta'),
('index', 'cylindrical_z'),
('index', 'dx'),
('index', 'dy'),
('index', 'dz'),
('index', 'grid_indices'),
('index', 'grid_level'),
('index', 'morton_index'),
('index', 'ones'),
('index', 'ones_over_dx'),
('index', 'path_element_x'),
('index', 'path_element_y'),
('index', 'path_element_z'),
('index', 'radius'),
('index', 'spherical_phi'),
('index', 'spherical_radius'),
('index', 'spherical_theta'),
('index', 'virial_radius_fraction'),
('index', 'volume'),
('index', 'x'),
('index', 'y'),
('index', 'z'),
('index', 'zeros'),
('io', 'particle_mass'),
('io', 'particle_ones'),
('io', 'particle_position'),
('io', 'particle_position_cylindrical_radius'),
('io', 'particle_position_cylindrical_theta'),
('io', 'particle_position_cylindrical_z'),
('io', 'particle_position_relative_x'),
('io', 'particle_position_relative_y'),
('io', 'particle_position_relative_z'),
('io', 'particle_position_spherical_phi'),
('io', 'particle_position_spherical_radius'),
('io', 'particle_position_spherical_theta'),
('io', 'particle_position_x'),
('io', 'particle_position_y'),
('io', 'particle_position_z'),
('io', 'particle_radius'),
('io', 'relative_particle_position'),
('io', 'relative_particle_position_x'),
('io', 'relative_particle_position_y'),
('io', 'relative_particle_position_z'),
('nbody', 'particle_mass'),
('nbody', 'particle_ones'),
('nbody', 'particle_position'),
('nbody', 'particle_position_cylindrical_radius'),
('nbody', 'particle_position_cylindrical_theta'),
('nbody', 'particle_position_cylindrical_z'),
('nbody', 'particle_position_relative_x'),
('nbody', 'particle_position_relative_y'),
('nbody', 'particle_position_relative_z'),
('nbody', 'particle_position_spherical_phi'),
('nbody', 'particle_position_spherical_radius'),
('nbody', 'particle_position_spherical_theta'),
('nbody', 'particle_position_x'),
('nbody', 'particle_position_y'),
('nbody', 'particle_position_z'),
('nbody', 'particle_radius'),
('nbody', 'relative_particle_position'),
('nbody', 'relative_particle_position_x'),
('nbody', 'relative_particle_position_y'),
('nbody', 'relative_particle_position_z'),
('stream', 'cell_volume'),
('stream', 'dx'),
('stream', 'dy'),
('stream', 'dz'),
('stream', 'path_element_x'),
('stream', 'path_element_y'),
('stream', 'path_element_z'),
('stream', 'volume'),
('stream', 'x'),
('stream', 'y'),
('stream', 'z')]
[7]:
slc = yt.SlicePlot(ds, 2, ("deposit", "all_cic"))
slc.set_width((8, "Mpc"))
---------------------------------------------------------------------------
NeedsGridType Traceback (most recent call last)
File /tmp/yt/yt/data_objects/data_containers.py:287, in YTDataContainer._generate_fluid_field(self, field)
286 try:
--> 287 finfo.check_available(gen_obj)
288 except NeedsGridType as ngt_exception:
File /tmp/yt/yt/fields/derived_field.py:235, in DerivedField.check_available(self, data)
234 for validator in self.validators:
--> 235 validator(data)
236 # If we don't get an exception, we're good to go
File /tmp/yt/yt/fields/derived_field.py:620, in ValidateSpatial.__call__(self, data)
619 if not getattr(data, "_spatial", False):
--> 620 raise NeedsGridType(self.ghost_zones, self.fields)
621 if self.ghost_zones <= data._num_ghost_zones:
NeedsGridType: fields None require 0 ghost zones.
During handling of the above exception, another exception occurred:
YTNonIndexedDataContainer Traceback (most recent call last)
Cell In[7], line 1
----> 1 slc = yt.SlicePlot(ds, 2, ("deposit", "all_cic"))
2 slc.set_width((8, "Mpc"))
File /tmp/yt/yt/visualization/plot_window.py:1839, in AxisAlignedSlicePlot.__init__(self, ds, normal, fields, center, width, axes_unit, origin, fontsize, field_parameters, window_size, aspect, data_source, buff_size, north_vector)
1831 else:
1832 slc = ds.slice(
1833 axis,
1834 center[axis],
(...)
1837 data_source=data_source,
1838 )
-> 1839 slc.get_data(fields)
1840 validate_mesh_fields(slc, fields)
1841 PWViewerMPL.__init__(
1842 self,
1843 slc,
(...)
1851 geometry=ds.geometry,
1852 )
File /tmp/yt/yt/data_objects/selection_objects/data_selection_objects.py:220, in YTSelectionContainer.get_data(self, fields)
217 self.field_data[f].convert_to_units(finfos[f].output_units)
219 fields_to_generate += gen_fluids + gen_particles
--> 220 self._generate_fields(fields_to_generate)
221 for field in list(self.field_data.keys()):
222 if field not in ofields:
File /tmp/yt/yt/data_objects/selection_objects/data_selection_objects.py:252, in YTSelectionContainer._generate_fields(self, fields_to_generate)
250 fi = self.ds._get_field_info(field)
251 try:
--> 252 fd = self._generate_field(field)
253 if hasattr(fd, "units"):
254 fd.units.registry = self.ds.unit_registry
File /tmp/yt/yt/data_objects/data_containers.py:273, in YTDataContainer._generate_field(self, field)
271 tr = self._generate_particle_field(field)
272 else:
--> 273 tr = self._generate_fluid_field(field)
274 if tr is None:
275 raise YTCouldNotGenerateField(field, self.ds)
File /tmp/yt/yt/data_objects/data_containers.py:289, in YTDataContainer._generate_fluid_field(self, field)
287 finfo.check_available(gen_obj)
288 except NeedsGridType as ngt_exception:
--> 289 rv = self._generate_spatial_fluid(field, ngt_exception.ghost_zones)
290 else:
291 rv = finfo(gen_obj)
File /tmp/yt/yt/data_objects/data_containers.py:315, in YTDataContainer._generate_spatial_fluid(self, field, ngz)
313 o = self._current_chunk.objs[0]
314 if accumulate:
--> 315 rv = self.ds.arr(np.empty(o.ires.size, dtype="float64"), units)
316 outputs.append(rv)
317 ind = 0 # Does this work with mesh?
File /tmp/yt/yt/data_objects/selection_objects/data_selection_objects.py:451, in YTSelectionContainer.ires(self)
449 if self._current_chunk is None:
450 self.index._identify_base_chunk(self)
--> 451 return self._current_chunk.ires
File /tmp/yt/yt/geometry/geometry_handler.py:272, in cacheable_property.<locals>.cacheable_func(self)
270 return getattr(self, n)
271 if self.data_size is None:
--> 272 tr = self._accumulate_values(n[1:])
273 else:
274 tr = func(self)
File /tmp/yt/yt/geometry/geometry_handler.py:308, in YTDataChunk._accumulate_values(self, method)
306 for obj in self._fast_index or self.objs:
307 f = getattr(obj, mname)
--> 308 arrs.append(f(self.dobj))
309 if method == "dtcoords":
310 arrs = [arr[0] for arr in arrs]
File /tmp/yt/yt/data_objects/index_subobjects/particle_container.py:17, in _non_indexed.<locals>._func_non_indexed(self, *args, **kwargs)
16 def _func_non_indexed(self, *args, **kwargs):
---> 17 raise YTNonIndexedDataContainer(self)
YTNonIndexedDataContainer: The data container type (ParticleContainer) is an unindexed type. Operations such as ires, icoords, fcoords and fwidth will not work on it.
Did you just attempt to perform an off-axis operation ? Be sure to consult the latest documentation to see whether the operation you tried is actually supported for your data type.
Finally, one can specify multiple particle types in the data directory by setting the field names to be field tuples (the default field type for particles is "io") if one is not specified:
[8]:
n_star_particles = 1000000
n_dm_particles = 2000000
ppxd, ppyd, ppzd = 1e6 * np.random.normal(size=[3, n_dm_particles])
ppmd = np.ones(n_dm_particles)
ppxs, ppys, ppzs = 5e5 * np.random.normal(size=[3, n_star_particles])
ppms = 0.1 * np.ones(n_star_particles)
data2 = {
("dm", "particle_position_x"): ppxd,
("dm", "particle_position_y"): ppyd,
("dm", "particle_position_z"): ppzd,
("dm", "particle_mass"): ppmd,
("star", "particle_position_x"): ppxs,
("star", "particle_position_y"): ppys,
("star", "particle_position_z"): ppzs,
("star", "particle_mass"): ppms,
}
ds2 = yt.load_particles(
data2, length_unit=1.0 * parsec, mass_unit=1e8 * Msun, n_ref=256, bbox=bbox
)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[8], line 21
8 ppms = 0.1 * np.ones(n_star_particles)
10 data2 = {
11 ("dm", "particle_position_x"): ppxd,
12 ("dm", "particle_position_y"): ppyd,
(...)
18 ("star", "particle_mass"): ppms,
19 }
---> 21 ds2 = yt.load_particles(
22 data2, length_unit=1.0 * parsec, mass_unit=1e8 * Msun, n_ref=256, bbox=bbox
23 )
TypeError: load_particles() got an unexpected keyword argument 'n_ref'
We now have separate "dm" and "star" particles, as well as their deposited fields:
[9]:
slc = yt.SlicePlot(ds2, 2, [("deposit", "dm_cic"), ("deposit", "star_cic")])
slc.set_width((8, "Mpc"))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 1
----> 1 slc = yt.SlicePlot(ds2, 2, [("deposit", "dm_cic"), ("deposit", "star_cic")])
2 slc.set_width((8, "Mpc"))
NameError: name 'ds2' is not defined