Ctrl+K

Quick-start: first steps with pynbody

Quick-start: first steps with pynbody#

Pynbody includes some essential tools that allow you to retrieve, explore and visualize physically interesting quantities. The most basic functionality is found in a pair of classes:

  • the SimSnap class, which makes the data from a simulation snapshot available, as well as providing facilities to convert units and perform transformations;

  • the HaloCatalogue class, which provides access to halo catalogues and the halos they contain.

But to make best use of pynbody it is also necessary to use the pynbody.analysis and pynbody.plot modules. This brief walkthrough will show you some of the capabilities in these, and subsequent walkthroughs expand in more detail.

In all walkthroughs in the pynbody documentation we will use the ipython interpreter which offers a much richer interactive environment over the vanilla python interpreter. This is also the same as using a Jupyter notebook. However, you can type exactly the same commands into vanilla python; only the formatting will look slightly different.

Note

Before you start make sure pynbody is properly installed. See Pynbody Installation for more information. You will also need the standard pynbody test files, so that you can load the exact same data as used to write the tutorial. You need to download these separately here: testdata.tar.gz. You can then extract them in a directory of your choice with tar -zxvf testdata.tar.gz

Code snippets can be copied from this page and pasted into python, ipython or jupyter. Hover over the code and click the button that appears; the commands will be copied to your clipboard.

Loading snapshots and halos#

The first step of any analysis is to load the data. Afterwards, we will want to center it on the halo of interest (in this case the main halo) to analyze its contents.

In [1]: import pynbody

In [2]: import pylab

In [3]: s = pynbody.load('testdata/gasoline_ahf/g15784.lr.01024.gz')

This loads the snapshot s (make sure you use the correct path to the testdata directory). The object s is an instance of SimSnap, and handles loading data on your behalf. Note that on the initial call to load, only the header is loaded. Actual data will be loaded only when it is needed.

Now we load the halos:

In [4]: h = s.halos()

Note that h is now a HaloCatalogue object, containing the halo catalogue which pynbody could locate for the snapshot s. Generally, it is not necessary to tell pynbody what format the halos are in or their location; it can infer it in most situations.

For later convenience, we can store the first halo in a separate variable. The particular snapshot we have loaded here is a zoom cosmological simulation, and halo 0 contains the central galaxy.

In [5]: main_halo = h[0]

Note

The halo numbers by default are those used by the halo finder, which (depending on your specific finder) may not start at zero, and may even be random numbers! You can see all the available halos using h.keys().

Older versions of pynbody renumbered AHF halos to start at 1, regardless of the internal numbering used by AHF. This inconsistency has been fixed in version 2, but to get the same results as in the previous versions, you need to specifically request it. h = s.halos(halo_number='v1') provides this backwards-compatibility.

We can check quickly how many particles of each type are identified there:

In [6]: print('ngas = %e, ndark = %e, nstar = %e\n'%(len(main_halo.gas),len(main_halo.dark),len(main_halo.star)))
ngas = 7.906000e+04, ndark = 1.610620e+05, nstar = 2.621780e+05

pynbody refers to different particle types as “families”. Here, we have accessed the gas, dark and star families of the halo. There are also convenient one-letter aliases for these regularly-used families: .g, .d and .s respectively. And, as you might expect, the python len function returns the number of particles in each family.

We could similarly have applied similar code to the entire snapshot, or to any other halo:

In [7]: print('Whole snapshot ngas = %e, ndark = %e, nstar = %e\n'%(len(s.gas),len(s.dark),len(s.star)))
Whole snapshot ngas = 1.587550e+05, ndark = 1.293231e+06, nstar = 2.651700e+05


In [8]: print('Halo 5 ngas = %e, ndark = %e, nstar = %e\n'%(len(h[5].gas),len(h[5].dark),len(h[5].star)))
Halo 5 ngas = 1.880000e+02, ndark = 1.009500e+04, nstar = 0.000000e+00

See also

Making some images#

Let’s skip straight to making some images. The following code will make a simple density interpolation of the gas particles in the main halo.

In [9]: s.physical_units()

In [10]: pynbody.analysis.center(main_halo)
Out[10]: <Transformation translate, offset_velocity>

In [11]: image_values = pynbody.plot.image(main_halo.gas, width=100, cmap='Blues')
plots/snapshot_manipulation_fig1.png

This has used three of pynbody’s routines:

  1. physical_units() to convert the units of the snapshot to physical units (unless otherwise specified, this means kpc, Msol and km/s for distances masses and velocities respectively);

  2. pynbody.analysis.center() to center the halo on the central density peak of the halo;

  3. pynbody.plot.image() to make an SPH-interpolated image of main_halo gas particles.

The latter automatically estimates smoothing lengths and densities if needed, even if these are not stored in the file explicitly. The returned image_values from image() is a numpy array of the pixel values, which you can then manipulate further if you wish.

Here’s another example showing the larger-scale dark-matter distribution – note that you can conveniently specify the width as a string with a unit. The units keyword is used to specify the units of the output, and notice here that we have specified a mass per unit area, which pynbody takes as an indication that we want a projected density map (rather than a slice through z=0, which is what we obtained in the gas case above).

In [12]: pynbody.plot.image(s.d[pynbody.filt.Sphere('10 Mpc')],
   ....:                    width='10 Mpc', units = 'Msol kpc^-2',
   ....:                    cmap='Greys')
   ....:
plots/snapshot_manipulation_fig1_wide.png

See also

See the Pictures in Pynbody tutorial for more examples and help regarding images.

pynbody also has a companion package, topsy, which enables real-time rendering of snapshots on a GPU. See its separate website for more information.

Aligning the Snapshot#

In the above example, the disk seems to be aligned more or less face-on. Pynbody images are always in the x-y plane; if they are projected, then the z-axis is the line of sight. To cut or project the simulation along another direction, we need to align it. For example, to align the disk, we can use the sideon() function:

In [13]: pynbody.analysis.sideon(main_halo)
Out[13]: <Transformation translate, offset_velocity, rotate>

In [14]: pynbody.plot.image(main_halo.g, width=100, cmap='Blues');
plots/snapshot_manipulation_fig2.png

Note that the function sideon() also calls center() to center the halo, so it doesn’t matter if the halo isn’t centered when you start. It then calculates the angular momentum vector in a sphere around the center and rotates the snapshot such that the angular momentum vector is parallel to the y-axis. If, instead, you’d like the disk face-on, you can call the equivalent pynbody.analysis.faceon().

Note

High-level snapshot manipulation functions defined in pynbody.analysis transform the entire simulation, even if you only pass in a subset of particles like a halo. That is why we could pass main_halo to center() but still plot all the dark matter particles in the simulation in the example in the previous section. The particles in main_halo were used to calculate the right center, but the transformation was applied to all particles. If this is not the behaviour you want, you can pass move_all = False to these routines, and only the particles you pass in will be transformed.

By contrast, core routines (i.e. those that are not part of the pynbody.analysis module) always operate on exactly what you apply them to, so s.g.rotate_x(90) rotates only the gas while s.rotate_x(90) rotates the entire simulation.

See also

See the next tutorial’s section on centering and reference documentation for the transformation module for more information about how coordinate transformations are handled in pynbody, including how to revert back to the original orientation.

Quick-look at the data and units#

Most analyses require you to get closer to the raw data arrays, and pynbody makes these readily accessible through a dictionary-like interface. The 3D position array is always known as pos, the velocity array as vel, and the mass array as mass. The units of these arrays are accessible through the units attribute, and may be converted to something more useful using the in_units method.

In [15]: s['pos']
Out[15]: 
SimArray([[-7.33556137e+02, -3.06744657e+03,  1.02756519e+02],
          [-6.86376128e+02, -3.05824610e+03, -1.67443704e+02],
          [-6.70360617e+02, -3.09768352e+03, -1.44824057e+02],
          ...,
          [-2.31221493e+00, -5.06106268e-01, -4.50602701e+00],
          [-1.59198626e-03, -2.84809511e-02, -1.01323468e+00],
          [ 7.31314832e+00,  4.49151520e-01,  2.23054824e+00]], 'kpc')

In [16]: s['pos'].units
Out[16]: Unit("kpc")

Earlier on, we converted the snapshot to physical units. We can easily undo that and see the data in its original units:

In [17]: s.original_units()

In [18]: s['pos']
Out[18]: 
SimArray([[-1.07099119e-02, -4.47846877e-02,  1.50024409e-03],
          [-1.00210843e-02, -4.46503609e-02, -2.44467632e-03],
          [-9.78725797e-03, -4.52261468e-02, -2.11442971e-03],
          ...,
          [-3.37583137e-05, -7.38914619e-06, -6.57879469e-05],
          [-2.32429826e-08, -4.15821587e-07, -1.47932157e-05],
          [ 1.06771889e-04,  6.55760748e-06,  3.25659809e-05]], '6.85e+04 kpc a')

Equally, we can manually convert units to whatever we wish:

In [19]: s['pos'].in_units('Mpc')
Out[19]: 
SimArray([[-7.33556137e-01, -3.06744657e+00,  1.02756519e-01],
          [-6.86376128e-01, -3.05824610e+00, -1.67443704e-01],
          [-6.70360617e-01, -3.09768352e+00, -1.44824057e-01],
          ...,
          [-2.31221493e-03, -5.06106268e-04, -4.50602701e-03],
          [-1.59198626e-06, -2.84809511e-05, -1.01323468e-03],
          [ 7.31314832e-03,  4.49151520e-04,  2.23054824e-03]], 'Mpc')

In [20]: s['pos'].in_units('Mpc a h**-1')
Out[20]: 
SimArray([[-5.35580029e-01, -2.23958746e+00,  7.50240321e-02],
          [-5.01133217e-01, -2.23287006e+00, -1.22253089e-01],
          [-4.89440059e-01, -2.26166389e+00, -1.05738155e-01],
          ...,
          [-1.68818183e-03, -3.69515564e-04, -3.28991600e-03],
          [-1.16233237e-06, -2.07943576e-05, -7.39777411e-04],
          [ 5.33943619e-03,  3.27932073e-04,  1.62855579e-03]], 'Mpc a h**-1')

Note here that the a is the cosmological expansion factor, i.e. its appearance in a unit indicates that the unit is comoving. The h is the Hubble parameter in units of 100 km/s/Mpc. The in_units() method makes a copy of the array in the new units, leaving the original array unchanged. There is also a convert_units() method that changes the units of the array in-place.

Now let’s convert the entire snapshot back to kpc, Msol and km/s, and check the units of the pos array again:

In [21]: s.physical_units()

In [22]: s['pos']
Out[22]: 
SimArray([[-7.33556137e+02, -3.06744657e+03,  1.02756519e+02],
          [-6.86376128e+02, -3.05824610e+03, -1.67443704e+02],
          [-6.70360617e+02, -3.09768352e+03, -1.44824057e+02],
          ...,
          [-2.31221493e+00, -5.06106268e-01, -4.50602701e+00],
          [-1.59198626e-03, -2.84809511e-02, -1.01323468e+00],
          [ 7.31314832e+00,  4.49151520e-01,  2.23054824e+00]], 'kpc')

Of course, vel and mass arrays can be handled in exactly the same way. Pynbody also loads all the other arrays inside a snapshot, standardizing the names where possible. If no standardized name is available, the array is loaded with the name it has in the snapshot file.

See also

Making a density profile#

Another component of pynbody’s scientific analysis tools is the ability to make profiles of any quantity. The pynbody.analysis.profile module is powerful and flexible, but here we will simply make a simple density profile of the gas, dark matter, and stars in the main halo.

Remember that the halo is already centred on the origin. We can therefore make 3d density profiles as follows:

In [23]: star_profile = pynbody.analysis.Profile(main_halo.s, min=0.2, max=50,
   ....:                                         type='log', nbins=50, ndim=3)
   ....: 

In [24]: dm_profile = pynbody.analysis.Profile(main_halo.d, min=0.2, max=50,
   ....:                                       type='log', nbins=50, ndim=3)
   ....: 

In [25]: gas_profile = pynbody.analysis.Profile(main_halo.g, min=0.2, max=50,
   ....:                                        type='log', nbins=50, ndim=3)
   ....:

The min and max arguments specify the minimum and maximum radii of the profile, and the nbins argument specifies the number of bins. The type argument specifies the binning scheme, which can be ‘log’, ‘lin’ or ‘equaln’. Finally, the ndim argument specifies the dimensionality. Note the use of the s, d and g shortcuts for the star, dark matter and gas families respectively.

Let’s now plot the profiles:

In [26]: pylab.plot(star_profile['rbins'], star_profile['density'], 'r', label='Stars')
Out[26]: [<matplotlib.lines.Line2D at 0x7f196b915bd0>]

In [27]: pylab.plot(dm_profile['rbins'], dm_profile['density'], 'k', label='Dark Matter')
Out[27]: [<matplotlib.lines.Line2D at 0x7f196ba51f10>]

In [28]: pylab.plot(gas_profile['rbins'], gas_profile['density'], 'b', label='Gas')
Out[28]: [<matplotlib.lines.Line2D at 0x7f196ba53810>]

In [29]: pylab.loglog()
Out[29]: []

In [30]: pylab.xlabel('r [kpc]')
Out[30]: Text(0.5, 0, 'r [kpc]')

In [31]: pylab.ylabel(r'$\rho$ [M$_\odot$/kpc$^3$]')
Out[31]: Text(0, 0.5, '$\\rho$ [M$_\\odot$/kpc$^3$]')

In [32]: pylab.legend()
Out[32]: <matplotlib.legend.Legend at 0x7f196b90c6d0>
plots/snapshot_manipulation_denpro.png

Where next?#

This tutorial has shown you how to load a snapshot, access the data, make some simple images and a density profile.

Contents