Selecting Extreme Halos and Visualizing (with Particles)

In this tutorial, we are going to define a measure of halo “extremeness”, filter halos by this quantity, and visualize a subset of interesting halos using the Opencosmo Toolkit. This requires loading the halo properties and particles in a collection.

Data Halo Properties Halo Particles Linked Datasets

Tasks Query Add Columns Offset Visualization

import numpy as np
import matplotlib.pyplot as plt

import opencosmo as oc

ds = oc.open("halo_properties_and_particles.hdf5")

There are many different ways that we can classify “extremeness” for a halo. One way is to look at the distance between the center-of-mass position and the potential minimum of the halo. We could compute this quantity as a derived column using the evaluate function; however, the offset_3d function already exists to compute 3D offsets between different positions.

from opencosmo.column import offset_3d

# compute 3D offset, normalized by R200c
offset_column = offset_3d("fof_halo_center", "sod_halo_com") / oc.col("sod_halo_radius")

# insert 3D offset as a new column
ds = ds.with_new_columns(dataset="halo_properties", xoff=offset_column)

Let’s now sort the halos by this new quantity “xoff”, and select the 3 halos with the largest offset.

ds_extreme = ds.sort_by("xoff", invert=True).take(3, at="start")

Now we can fetch their halo IDs and visualize with halo_projection_array, which returns a matplotlib Figure object. :::{note} For technical reasons, you cannot use plt.show() to display images generated with visualize_halo and halo_projection_array. Instead, save the figure using fig.savefig() and display with IPython.display.Image.

from opencosmo.analysis import halo_projection_array

halo_ids = ds_extreme["halo_properties"].get_data()["unique_tag"].value
halo_projection_array(halo_ids, ds_extreme).savefig("extreme_halos_dm.png")

This created a 1x3 array of DM particle projections

from IPython.display import Image, display
display(Image(filename="extreme_halos_dm.png"))

To instead look at stars, we can set field = ("star", "particle_mass")

fig = halo_projection_array(halo_ids, ds_extreme, field=("star", "particle_mass")).savefig("extreme_halos_stars.png")

display(Image(filename="extreme_halos_stars.png"))

Now, let’s make a 3x3 image that shows dark matter, gas and stars for all 3 halos. We will also add length scale indicators using the length_scale parameter.

halo_ids_3x3 = ( halo_ids, halo_ids, halo_ids )

params = {
    "fields": (
        3*[("dm", "particle_mass")],
        3*[("gas", "particle_mass")],
        3*[("star", "particle_mass")]
    ),
    "labels": (
        ["Dark Matter", None, None],
        ["Gas", None, None],
        ["Stars", None, None]
    ),
    "cmaps": (
        3*["gray"],
        3*["cividis"], 
        3*["bone"]
    ),
}

fig = halo_projection_array(halo_ids_3x3, ds_extreme, params=params, length_scale="all top")
fig.savefig("extreme_halos.png")

display(Image(filename="extreme_halos.png"))

Say we only wanted to look at one halo. halo_projection_array can be used to look at one or multiple fields for the halo. Another option is the visualize_halo function, which plots dark matter, gas, stars, and gas temperature with minimal input.

from opencosmo.analysis import visualize_halo

visualize_halo(halo_ids[0], ds_extreme).savefig("visualize_halo_0.png")

display(Image(filename="visualize_halo_0.png"))

Visualizing a Custom Field¶

Now that we are familiar with the tools, let’s take it a step further. Let’s derive a new column for radial velocity and visualize it for the most massive halos in the dataset.

from matplotlib.colors import SymLogNorm

def radial_velocity(ds):
    '''
    Computes radial velocity for all particle types, relative to the halo's center of mass.
    '''
    
    keys = ["dm_particles", "gas_particles", "star_particles", "agn_particles"]

    halo_properties = ds["halo_properties"].get_data()
    for key in keys:
        ptype = ds[key].get_data()

        x, y, z = ptype["x"] - halo_properties["sod_halo_com_x"], ptype["y"] - halo_properties["sod_halo_com_y"],ptype["z"] - halo_properties["sod_halo_com_z"]
        vx, vy, vz = ptype["vx"] - halo_properties["sod_halo_com_vx"], ptype["vy"] - halo_properties["sod_halo_com_vy"], ptype["vz"] - halo_properties["sod_halo_com_vz"]

        # direction given by v dot r_unit
        denom = 1/np.sqrt(np.square(x) + np.square(y) + np.square(z))
        runit = (x*denom,y*denom,z*denom)

        vr = vx*runit[0] + vy*runit[1] + vz*runit[2]

        # add derived column for this particle type
        ds = ds.with_new_columns(dataset=key, vr=vr)

    return ds
    
# fetch data for the most massive halo
ds_most_massive = radial_velocity(ds.sort_by("sod_halo_mass").take(1, at="end"))

halo_ids = 3*[ ds_most_massive["halo_properties"].get_data()["unique_tag"] ]
params = {
    # the derived column is accessible through yt via the (<particle_type>, "vr") key
    "fields": (
        [("dm", "vr"), ("gas", "vr"), ("star", "vr")],
    ),
    "labels": (
        ["Dark Matter", "Gas", "Stars"],
    ),
    "cmaps": (
        3*["PiYG"],
    ),
    # use SymLogNorm because radial velocity can be negative
    "cmap_norms": (
        3*[SymLogNorm(1)],
    ),
}

fig = halo_projection_array(halo_ids, ds_most_massive, params=params, length_scale="all")
fig.savefig("radial_velocity_most_massive.png")

display(Image(filename="radial_velocity_most_massive.png"))

The above image shows radial velocity for dark matter, gas, and star particles for the most massive halo in the dataset, color-coded such that inflowing particles are pink and outflowing particles are green.

Viewing list of available yt fields¶

The halo_projection_array function uses yt under-the-hood, so input fields must follow yt’s naming conventions. To view the list of available fields to visualize (and their format), you can directly load a given halo’s particle data into yt using the create_yt_dataset function:

from opencosmo.analysis import create_yt_dataset

# create yt dataset
ds_yt = create_yt_dataset(next(ds_most_massive.halos()))

# print all yt fields, including additional derived fields
for f in ds_yt.derived_field_list: print(f)

('agn', 'age')
('agn', 'bhr')
('agn', 'ebh')
('agn', 'fof_halo_tag')
('agn', 'id')
('agn', 'mask')
('agn', 'mbh')
('agn', 'particle_mass')
('agn', 'particle_ones')
('agn', 'particle_position')
('agn', 'particle_position_cylindrical_radius')
('agn', 'particle_position_cylindrical_theta')
('agn', 'particle_position_cylindrical_z')
('agn', 'particle_position_relative_x')
('agn', 'particle_position_relative_y')
('agn', 'particle_position_relative_z')
('agn', 'particle_position_spherical_phi')
('agn', 'particle_position_spherical_radius')
('agn', 'particle_position_spherical_theta')
('agn', 'particle_position_x')
('agn', 'particle_position_y')
('agn', 'particle_position_z')
('agn', 'particle_radius')
('agn', 'phi')
('agn', 'relative_particle_position')
('agn', 'relative_particle_position_x')
('agn', 'relative_particle_position_y')
('agn', 'relative_particle_position_z')
('agn', 'vr')
('agn', 'vx')
('agn', 'vy')
('agn', 'vz')
('all', 'fof_halo_tag')
('all', 'id')
('all', 'mask')
('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', 'phi')
('all', 'relative_particle_position')
('all', 'relative_particle_position_x')
('all', 'relative_particle_position_y')
('all', 'relative_particle_position_z')
('all', 'vr')
('all', 'vx')
('all', 'vy')
('all', 'vz')
('deposit', 'agn_cic')
('deposit', 'agn_cic_age')
('deposit', 'agn_count')
('deposit', 'agn_density')
('deposit', 'agn_mass')
('deposit', 'agn_nn_age')
('deposit', 'all_cic')
('deposit', 'all_count')
('deposit', 'all_density')
('deposit', 'all_mass')
('deposit', 'dm_cic')
('deposit', 'dm_count')
('deposit', 'dm_density')
('deposit', 'dm_mass')
('deposit', 'gas_cic')
('deposit', 'gas_count')
('deposit', 'gas_density')
('deposit', 'gas_mass')
('deposit', 'nbody_cic')
('deposit', 'nbody_count')
('deposit', 'nbody_density')
('deposit', 'nbody_mass')
('deposit', 'star_cic')
('deposit', 'star_cic_age')
('deposit', 'star_count')
('deposit', 'star_density')
('deposit', 'star_mass')
('deposit', 'star_nn_age')
('dm', 'fof_halo_tag')
('dm', 'id')
('dm', 'mask')
('dm', 'particle_mass')
('dm', 'particle_ones')
('dm', 'particle_position')
('dm', 'particle_position_cylindrical_radius')
('dm', 'particle_position_cylindrical_theta')
('dm', 'particle_position_cylindrical_z')
('dm', 'particle_position_relative_x')
('dm', 'particle_position_relative_y')
('dm', 'particle_position_relative_z')
('dm', 'particle_position_spherical_phi')
('dm', 'particle_position_spherical_radius')
('dm', 'particle_position_spherical_theta')
('dm', 'particle_position_x')
('dm', 'particle_position_y')
('dm', 'particle_position_z')
('dm', 'particle_radius')
('dm', 'phi')
('dm', 'relative_particle_position')
('dm', 'relative_particle_position_x')
('dm', 'relative_particle_position_y')
('dm', 'relative_particle_position_z')
('dm', 'vr')
('dm', 'vx')
('dm', 'vy')
('dm', 'vz')
('gas', 'MMW')
('gas', 'cell_mass')
('gas', 'cell_volume')
('gas', 'cylindrical_radius')
('gas', 'cylindrical_theta')
('gas', 'cylindrical_z')
('gas', 'density')
('gas', 'dx')
('gas', 'dy')
('gas', 'dynamical_time')
('gas', 'dz')
('gas', 'fof_halo_tag')
('gas', 'id')
('gas', 'mask')
('gas', 'mass')
('gas', 'mean_molecular_weight')
('gas', 'metallicity')
('gas', 'mu')
('gas', 'number_density')
('gas', 'ones')
('gas', 'particle_mass')
('gas', 'particle_ones')
('gas', 'particle_position')
('gas', 'particle_position_cylindrical_radius')
('gas', 'particle_position_cylindrical_theta')
('gas', 'particle_position_cylindrical_z')
('gas', 'particle_position_relative_x')
('gas', 'particle_position_relative_y')
('gas', 'particle_position_relative_z')
('gas', 'particle_position_spherical_phi')
('gas', 'particle_position_spherical_radius')
('gas', 'particle_position_spherical_theta')
('gas', 'particle_position_x')
('gas', 'particle_position_y')
('gas', 'particle_position_z')
('gas', 'particle_radius')
('gas', 'path_element_x')
('gas', 'path_element_y')
('gas', 'path_element_z')
('gas', 'phi')
('gas', 'position')
('gas', 'radius')
('gas', 'relative_particle_position')
('gas', 'relative_particle_position_x')
('gas', 'relative_particle_position_y')
('gas', 'relative_particle_position_z')
('gas', 'relative_x')
('gas', 'relative_y')
('gas', 'relative_z')
('gas', 'smoothing_length')
('gas', 'spherical_phi')
('gas', 'spherical_radius')
('gas', 'spherical_theta')
('gas', 'temperature')
('gas', 'uu')
('gas', 'volume')
('gas', 'vr')
('gas', 'vx')
('gas', 'vy')
('gas', 'vz')
('gas', 'x')
('gas', 'xh')
('gas', 'y')
('gas', 'yhe')
('gas', 'z')
('gas', 'zmet')
('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', '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')
('nbody', 'fof_halo_tag')
('nbody', 'id')
('nbody', 'mask')
('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', 'phi')
('nbody', 'relative_particle_position')
('nbody', 'relative_particle_position_x')
('nbody', 'relative_particle_position_y')
('nbody', 'relative_particle_position_z')
('nbody', 'vr')
('nbody', 'vx')
('nbody', 'vy')
('nbody', 'vz')
('star', 'age')
('star', 'fof_halo_tag')
('star', 'id')
('star', 'initial_star_mass')
('star', 'mask')
('star', 'particle_mass')
('star', 'particle_ones')
('star', 'particle_position')
('star', 'particle_position_cylindrical_radius')
('star', 'particle_position_cylindrical_theta')
('star', 'particle_position_cylindrical_z')
('star', 'particle_position_relative_x')
('star', 'particle_position_relative_y')
('star', 'particle_position_relative_z')
('star', 'particle_position_spherical_phi')
('star', 'particle_position_spherical_radius')
('star', 'particle_position_spherical_theta')
('star', 'particle_position_x')
('star', 'particle_position_y')
('star', 'particle_position_z')
('star', 'particle_radius')
('star', 'phi')
('star', 'relative_particle_position')
('star', 'relative_particle_position_x')
('star', 'relative_particle_position_y')
('star', 'relative_particle_position_z')
('star', 'vr')
('star', 'vx')
('star', 'vy')
('star', 'vz')
('star', 'yhe')
('star', 'zmet')
('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')