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; see Obtaining test data.
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:
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.
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
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 null, translate, offset_velocity> In [11]: image_values = pynbody.plot.image(main_halo.gas, width=100, cmap='Blues')
This has used three of pynbody’s routines:
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);
pynbody.analysis.center() to center the halo on the central density peak of the halo;
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')
....:
See also
See the pictures 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.
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 sideon> In [14]: pynbody.plot.image(main_halo.g, width=100, cmap='Blues');
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.
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]],
shape=(1717156, 3), '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]],
shape=(1717156, 3), '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]],
shape=(1717156, 3), '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]],
shape=(1717156, 3), '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]],
shape=(1717156, 3), '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.
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 0x7206d7d6cbd0>]
In [27]: pylab.plot(dm_profile['rbins'], dm_profile['density'], 'k', label='Dark Matter')
Out[27]: [<matplotlib.lines.Line2D at 0x7206d8239090>]
In [28]: pylab.plot(gas_profile['rbins'], gas_profile['density'], 'b', label='Gas')
Out[28]: [<matplotlib.lines.Line2D at 0x7206d8238ed0>]
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 0x7206d814e1d0>
Histograms, 2D histograms, and derived arrays#
Pynbody doesn’t make specific provision for plotting 1d histograms since this is well-supported by matplotlib and numpy. However, it is worth mentioning that there are often derived arrays that can be used to make calculating the right quantity easier.
For example, for many file formats pynbody is able to calculate the metallicity relative to solar, simply by requesting the 'feh' array.
In [33]: pylab.hist(main_halo.s['feh'], bins=100, histtype='step', color='r', label='Stars', range=(-3.0, 0.2), density=True)
....: pylab.hist(main_halo.g['feh'], bins=100, histtype='step', color='b', label='Gas', range=(-3.0, 0.2), density=True)
....: pylab.xlabel("[Fe/H]")
....: pylab.legend()
....:
Out[33]: <matplotlib.legend.Legend at 0x7206d83ea550>
Pynbody also provides a convenient hist2d() plotting routine, which can be used to make 2d histograms in a variety of ways.
In [34]: pynbody.plot.hist2d(main_halo.s['feh'], main_halo.s['ofe'], logscale=True, cmap='Blues', nbins=128,
....: x_range=(-3, 0.2), y_range=(-1, 1))
....: pylab.xlabel("[Fe/H]")
....: pylab.ylabel("[O/Fe]")
....:
Out[34]: Text(0, 0.5, '[O/Fe]')
Note that the arrays 'feh' and 'ofe' are not stored in the snapshot file, but are derived on demand. To check whether an array is derived you can use the is_derived_array() method:
In [35]: main_halo.s.is_derived_array('feh')
Out[35]: True
To find out how an array was derived, you can get the deriving function using find_deriving_function():
In [36]: help(main_halo.s.find_deriving_function('feh'))
Help on function feh in module pynbody.snapshot.tipsy:
feh(self)
Derived [TipsySnap]: Iron abundance [Fe/H] derived from tipsy array FeMassFrac, with solar
values from Asplund et al 09
Here we used python’s built-in help function to display information about the deriving function. You can see here that the derivation is made specifically for tipsy files.
There is also a special function, rho_T(), for making histograms of density-temperature relations, which handles supported multiphase ISM models:
In [37]: pynbody.plot.gas.rho_T(main_halo.g)
Out[37]:
(array([[1., 1., 1., ..., 1., 3., 2.],
[1., 1., 1., ..., 5., 1., 1.],
[1., 1., 1., ..., 1., 1., 1.],
...,
[1., 1., 1., ..., 1., 1., 1.],
[1., 1., 1., ..., 1., 1., 1.],
[1., 1., 1., ..., 1., 1., 1.]], shape=(100, 100)),
array([1.81563409, 1.8892105 , 1.96278691, 2.03636332, 2.10993972,
2.18351613, 2.25709254, 2.33066895, 2.40424536, 2.47782177,
2.55139818, 2.62497458, 2.69855099, 2.7721274 , 2.84570381,
2.91928022, 2.99285663, 3.06643304, 3.14000944, 3.21358585,
3.28716226, 3.36073867, 3.43431508, 3.50789149, 3.5814679 ,
3.6550443 , 3.72862071, 3.80219712, 3.87577353, 3.94934994,
4.02292635, 4.09650276, 4.17007916, 4.24365557, 4.31723198,
4.39080839, 4.4643848 , 4.53796121, 4.61153762, 4.68511402,
4.75869043, 4.83226684, 4.90584325, 4.97941966, 5.05299607,
5.12657248, 5.20014888, 5.27372529, 5.3473017 , 5.42087811,
5.49445452, 5.56803093, 5.64160734, 5.71518374, 5.78876015,
5.86233656, 5.93591297, 6.00948938, 6.08306579, 6.1566422 ,
6.2302186 , 6.30379501, 6.37737142, 6.45094783, 6.52452424,
6.59810065, 6.67167706, 6.74525346, 6.81882987, 6.89240628,
6.96598269, 7.0395591 , 7.11313551, 7.18671192, 7.26028832,
7.33386473, 7.40744114, 7.48101755, 7.55459396, 7.62817037,
7.70174678, 7.77532318, 7.84889959, 7.922476 , 7.99605241,
8.06962882, 8.14320523, 8.21678163, 8.29035804, 8.36393445,
8.43751086, 8.51108727, 8.58466368, 8.65824009, 8.73181649,
8.8053929 , 8.87896931, 8.95254572, 9.02612213, 9.09969854]),
array([2.3112598 , 2.36072272, 2.41018564, 2.45964856, 2.50911148,
2.5585744 , 2.60803732, 2.65750024, 2.70696316, 2.75642608,
2.805889 , 2.85535192, 2.90481484, 2.95427776, 3.00374068,
3.0532036 , 3.10266652, 3.15212944, 3.20159236, 3.25105529,
3.30051821, 3.34998113, 3.39944405, 3.44890697, 3.49836989,
3.54783281, 3.59729573, 3.64675865, 3.69622157, 3.74568449,
3.79514741, 3.84461033, 3.89407325, 3.94353617, 3.99299909,
4.04246201, 4.09192493, 4.14138785, 4.19085077, 4.24031369,
4.28977661, 4.33923953, 4.38870245, 4.43816537, 4.48762829,
4.53709122, 4.58655414, 4.63601706, 4.68547998, 4.7349429 ,
4.78440582, 4.83386874, 4.88333166, 4.93279458, 4.9822575 ,
5.03172042, 5.08118334, 5.13064626, 5.18010918, 5.2295721 ,
5.27903502, 5.32849794, 5.37796086, 5.42742378, 5.4768867 ,
5.52634962, 5.57581254, 5.62527546, 5.67473838, 5.7242013 ,
5.77366422, 5.82312715, 5.87259007, 5.92205299, 5.97151591,
6.02097883, 6.07044175, 6.11990467, 6.16936759, 6.21883051,
6.26829343, 6.31775635, 6.36721927, 6.41668219, 6.46614511,
6.51560803, 6.56507095, 6.61453387, 6.66399679, 6.71345971,
6.76292263, 6.81238555, 6.86184847, 6.91131139, 6.96077431,
7.01023723, 7.05970015, 7.10916308, 7.158626 , 7.20808892]))
Where next?#
This tutorial has shown you how to load a snapshot, access the data, make some simple images and a density profile.
In the next tutorial, we will go into more depth on how to manipulate the data inside a snapshot.
For more about images, see the Images cookbook.
For more about profiles, such as density profiles or rotation curves, see the Profiles walk-through.
For more about the low-level data access facilities, see the A deeper walk through pynbody’s data access facilities walk-through.
For more about halos, see the halos cookbook.
Or go back to the table of contents for all Pynbody Tutorials.
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4