Reading and manipulating snapshots

Example of how to read and manipulate CAMELS snapshots

import numpy as np
import h5py
import hdf5plugin

Get the name of the snapshot

f_snap = '/home/jovyan/Data/Sims/IllustrisTNG/CV/CV_0/snap_033.hdf5'

Open the file

data = h5py.File(f_snap, 'r')

Read the header of the snapshot

BoxSize      = data['Header'].attrs[u'BoxSize']/1e3    #size of the snapshot in comoving Mpc/h
redshift     = data['Header'].attrs[u'Redshift']       #reshift of the snapshot
scale_factor = data['Header'].attrs[u'Time']           #scale factor
h            = data['Header'].attrs[u'HubbleParam']    #value of the hubble parameter in 100 km/s/(Mpc/h)
Masses       = data['Header'].attrs[u'MassTable']*1e10 #masses of the particles in Msun/h
Np           = data['Header'].attrs[u'NumPart_Total']  #total number of particles for specie
Omega_m      = data['Header'].attrs[u'Omega0']         #Omega_matter
Omega_L      = data['Header'].attrs[u'OmegaLambda']    #Omega_baryon
Omega_b      = data['Header'].attrs[u'OmegaBaryon']    #Omega_Lambda

print('Box size:                 %.2f Mpc/h'%BoxSize)
print('snapshot redshift:        %.2f'%redshift)
print('Number of gas  particles: %d'%Np[0])
print('Number of star particles: %d'%Np[4])
print('Omega_m:                  %.3f'%Omega_m)
print('Omega_b:                  %.3f'%Omega_b)
print('Omega_L:                  %.3f'%Omega_L)

Box size:                 25.00 Mpc/h
snapshot redshift:        0.00
Number of gas  particles: 15695635
Number of star particles: 636474
Omega_m:                  0.300
Omega_b:                  0.049
Omega_L:                  0.700

Read the positions of the gas, dark matter, stars, and black-holes particles

pos_gas   = data['PartType0/Coordinates'][:]/1e3 #Mpc/h
pos_dm    = data['PartType1/Coordinates'][:]/1e3 #Mpc/h
pos_stars = data['PartType4/Coordinates'][:]/1e3 #Mpc/h
pos_bh    = data['PartType5/Coordinates'][:]/1e3 #Mpc/h

print('%.2f < pos_gas_X < %.2f'%(np.min(pos_gas[:,0]), np.max(pos_gas[:,0])))
print('%.2f < pos_gas_Y < %.2f'%(np.min(pos_gas[:,1]), np.max(pos_gas[:,1])))
print('%.2f < pos_gas_Z < %.2f'%(np.min(pos_gas[:,2]), np.max(pos_gas[:,2])))

0.00 < pos_gas_X < 25.00
0.00 < pos_gas_Y < 25.00
0.00 < pos_gas_Z < 25.00

print('Number of star particles: %d'%pos_stars.shape[0])

Number of star particles: 636474