diff --git a/source/service_docs/python_examples.rst b/source/service_docs/python_examples.rst index 65c5b7f..3722b98 100644 --- a/source/service_docs/python_examples.rst +++ b/source/service_docs/python_examples.rst @@ -10,11 +10,14 @@ in python. Here, we provide links to these examples. specified region, and save part of a snapshot to a local file ` + * :doc:`Read SOAP halo properties using swiftsimio` + + * :doc:`Identify a halo of interest and read in the particles + belonging to it ` + * :doc:`Read and plot the matter power spectrum at several redshifts ` - * :doc:`Read SOAP halo properties using swiftsimio` - * :doc:`Identify the same halo between different simulations ` * :ref:`Plot the evolution of a halo across snapshots ` diff --git a/source/snapshots/halo_particles.rst b/source/snapshots/halo_particles.rst new file mode 100644 index 0000000..40b8c0f --- /dev/null +++ b/source/snapshots/halo_particles.rst @@ -0,0 +1,136 @@ +Reading particles belonging to a halo +===================================== + +Here we show how to identify a halo of interest in a :doc:`SOAP halo +catalogue ` and then find its particles in the +corresponding :doc:`snapshot `. This could be used +to compute a halo density profile, for example. + +In the SOAP halo catalogues, each halo is assigned an index which is +unique within the snapshot (but not between snapshots). This index is +stored in the HDF5 dataset ``InputHalos/HaloCatalogueIndex``. There is +a corresponding index ``PartTypeX/HaloCatalogueIndex`` associated with +each particle in the snapshot. + +We can read in all particles belonging to a halo as follows: + + * Use swiftsimio to read the halo catalogue (see :doc:`/soap/swiftsimio`) + * Identify a halo of interest in the catalogue + * Extract the position and ``HaloCatalogueIndex`` of this halo from the SOAP catalogue + * Read particles from the snapshot in some radius around the halo position (see :ref:`swiftsimio_cutouts`) + * Discard particles which do not have the required ``HaloCatalogueIndex`` + +This can be done either by using the :doc:`hdfstream python module +` for remote access, or by downloading +the SOAP catalogue and snapshot as HDF5 files and reading them +directly with h5py. Both methods are illustrated below. + +.. tab-set:: + + .. tab-item:: Using remote file access + + .. code-block:: python + + import numpy as np + import unyt + import hdfstream + import swiftsimio as sw + + # Connect to the hdfstream service and open the root directory + root_dir = hdfstream.open("cosma", "/") + + # Open the z=0 halo catalogue from the L1_m10_DMO simulation + soap_file = root_dir["FLAMINGO/L1_m10/L1_m10_DMO/SOAP-HBT/halo_properties_0077.hdf5"] + soap = sw.load(soap_file) + + # Get halo positions, masses and indexes + halo_pos = soap.input_halos.halo_centre + halo_m200c = soap.spherical_overdensity_200_crit.total_mass + halo_index = soap.input_halos.halo_catalogue_index + + # Select the most massive halo + i = np.argmax(halo_m200c) + target_pos = halo_pos[i,:] + target_index = halo_index[i] + + # Choose a region to read in. Note that we need to choose a radius large + # enough to enclose all halo particles. + radius = 5.0*unyt.Mpc + region = [[x-radius, x+radius] for x in target_pos] + + # Open the z=0 snapshot from the L1_m10_DMO simulation and select this region + snap_file = root_dir["FLAMINGO/L1_m10/L1_m10_DMO/snapshots/flamingo_0077/flamingo_0077.hdf5"] + mask = sw.mask(snap_file) + mask.constrain_spatial(region) + snap = sw.load(snap_file, mask=mask) + + # Read position and halo index of dark matter particles in this region + dm_pos = snap.dark_matter.coordinates + dm_halo_index = snap.dark_matter.halo_catalogue_index + + # Of the particles we read, identify those which belong to the halo + in_halo = (dm_halo_index == target_index) + halo_dm_pos = dm_pos[in_halo,:] + + # Make a simple dotplot as a sanity check + import matplotlib.pyplot as plt + plt.plot(halo_dm_pos[:,0], halo_dm_pos[:,1], "k,") + plt.gca().set_aspect("equal") + plt.xlabel("x [Mpc]") + plt.ylabel("y [Mpc]") + plt.title("Most massive halo in L1_m10_DMO") + plt.show() + + .. tab-item:: Reading local HDF5 files + + .. code-block:: python + + import numpy as np + import unyt + import swiftsimio as sw + + # Open the z=0 halo catalogue from the L1_m10_DMO simulation + soap_file = "./FLAMINGO/L1_m10/L1_m10_DMO/SOAP-HBT/halo_properties_0077.hdf5" + soap = sw.load(soap_file) + + # Get halo positions, masses and indexes + halo_pos = soap.input_halos.halo_centre + halo_m200c = soap.spherical_overdensity_200_crit.total_mass + halo_index = soap.input_halos.halo_catalogue_index + + # Select the most massive halo + i = np.argmax(halo_m200c) + target_pos = halo_pos[i,:] + target_index = halo_index[i] + + # Choose a region to read in. Note that we need to choose a radius large + # enough to enclose all halo particles. + radius = 5.0*unyt.Mpc + region = [[x-radius, x+radius] for x in target_pos] + + # Open the z=0 snapshot from the L1_m10_DMO simulation and select this region + snap_file = "./FLAMINGO/L1_m10/L1_m10_DMO/snapshots/flamingo_0077/flamingo_0077.hdf5" + mask = sw.mask(snap_file) + mask.constrain_spatial(region) + snap = sw.load(snap_file, mask=mask) + + # Read position and halo index of dark matter particles in this region + dm_pos = snap.dark_matter.coordinates + dm_halo_index = snap.dark_matter.halo_catalogue_index + + # Of the particles we read, identify those which belong to the halo + in_halo = (dm_halo_index == target_index) + halo_dm_pos = dm_pos[in_halo,:] + + # Make a simple dotplot as a sanity check + import matplotlib.pyplot as plt + plt.plot(halo_dm_pos[:,0], halo_dm_pos[:,1], "k,") + plt.gca().set_aspect("equal") + plt.xlabel("x [Mpc]") + plt.ylabel("y [Mpc]") + plt.title("Most massive halo in L1_m10_DMO") + plt.show() + +This results in the following plot showing dark matter particles in the halo: + +.. image:: images/L1_m10_halo.png diff --git a/source/snapshots/images/L1_m10_halo.png b/source/snapshots/images/L1_m10_halo.png new file mode 100644 index 0000000..a5b3039 Binary files /dev/null and b/source/snapshots/images/L1_m10_halo.png differ diff --git a/source/snapshots/index.rst b/source/snapshots/index.rst index fe69e87..6178122 100644 --- a/source/snapshots/index.rst +++ b/source/snapshots/index.rst @@ -40,6 +40,7 @@ The following sections describe the layout and contents of the snapshots. Output redshifts Particle properties Reading with swiftsimio + Reading halo particles Partial snapshots For more information about the SWIFT simulation snapshot format used diff --git a/source/snapshots/swiftsimio.rst b/source/snapshots/swiftsimio.rst index bb0a38a..d8bf9d3 100644 --- a/source/snapshots/swiftsimio.rst +++ b/source/snapshots/swiftsimio.rst @@ -156,6 +156,8 @@ Mpc:: [987.07806218 995.72685218 88.39880218] [987.07910218 995.73126218 88.39935218]] Mpc (Comoving) +.. _swiftsimio_cutouts: + Cutting out a region -------------------- diff --git a/source/soap/swiftsimio.rst b/source/soap/swiftsimio.rst index 12aecdc..548a810 100644 --- a/source/soap/swiftsimio.rst +++ b/source/soap/swiftsimio.rst @@ -16,6 +16,9 @@ The latter method might be preferable if you only need a small fraction of the data, such as a subset of halo properties or halos in a small region of interest. +To read the particles associated with a halo selected from a SOAP +catalogue, see :doc:`/snapshots/halo_particles`. + Installation ------------