Astronomical Dendrograms in Python¶

Requirements¶

This package has the following depdenencies:

• Python 2.6 or later (Python 3.x is supported)
• Numpy 1.4.1 or later
• Astropy 0.2.0 or later, optional (needed for reading/writing FITS files and for analysis code)
• h5py 0.2.0 or later, optional (needed for reading/writing HDF5 files)

Installation¶

To install the latest stable release, you can type:

pip install astrodendro

or you can download the latest tar file from PyPI and install it using:

python setup.py install

Developer version¶

If you want to install the latest developer version of the dendrogram code, you can do so from the git repository:

git clone https://github.com/dendrograms/astrodendro.git
cd astrodendro
python setup.py install

You may need to add the --user option to the last line if you do not have root access.

Basic example¶

The following diagram shows a one-dimensional dataset (with flux versus position) in the solid black line, with the corresponding dendrogram for that dataset overplotted in green:

In the rest of this document, we will refer to the individual points in this dataset as pixels.

The way the algorithm works is to construct the tree starting from the brightest pixels in the dataset, and progressively adding fainter and fainter pixels. We will illustrate this by showing the current value being considered, with the following blue dashed line:

Let’s now start moving this line down, starting from the peak pixel in the dataset. We create our first structure from this pixel. We then move to the pixel with the next largest value, and each time, we decide whether to join the pixel to an existing structure, or create a new structure. We only start a new structure if the value of the pixel is greater than its immediate neighbors, and therefore is a local maximum. The first structure being constructed is shown below:

We have now found a local maximum, so rather than add this pixel to the first structure, we create a new structure. As we move further down, both structures keep growing, until we reach a pixel that is not a local maximum, and is adjacent to both existing structures:

At this point, we merge the structures into a branch, which is shown by a green horizontal line. As we move down further, the branch continues to grow, and we very quickly find two more local maxima which cause new structures to be created:

These structures eventually merge, and we end up with a single tree:

Accounting for noise¶

Setting a minimum value (min_value)¶

Most real-life datasets are likely to contain some level of noise, and below a certain value of the noise, there is no point in expanding the tree since it will not be measuring anything physical; new branches will be ‘noise spikes’. By default, the minimum value is set to negative infinity, which means all pixels are added to the tree. However, you will very likely want to change this so that only significant features above the noise are included.

Let’s go back to the original data. We have left the outline of the complete tree for reference. We now set a minimum value, which we show below with the purple line. This is controlled by the min_value option in compute().

The effect on the tree is simply to get rid of (or prune) any structure peaking below this minimum. In this case, the peak on the right is no longer part of the tree since it is below the minimum specified value.

Setting a minimum significance for structures (min_delta)¶

If our data are noisy, we also want to avoid including local maxima that - while above the minimum absolute value specified above - are only identified because of noise, so we can also define a minimum height required for a structure to be retained. This is the min_delta parameter in compute(). We show the value corresponding to the current value being considered plus this minimum height:

In this case, min_delta is set to 0.01. As we now move down in flux as before, the structure first appears red. This indicates that the structure is not yet part of the tree:

Once the height of the structure exceeds the minimum specified, the structure can now be considered part of the tree:

In this case, all structures that are above the minimum value are also all large enough to be significant, so the tree is the same as before:

We can now repeat this experiment, but this time, with a larger minimum height for structures to be retained (min_delta=0.025). Once we reach the point where the second peak would have been merged, we can see that it is not high enough above the merging point to be considered an independent structure:

and the pixels are then simply added to the first structure, rather than creating a branch:

We can now see that the final tree looks a little different to the original one, because the second largest peak was deemed insignificant:

Additional options¶

In addition to the minimum height of a structure, it is also possible to specify the minimum number of pixels that a structure should contain in order to remain an independent structure (min_npix), and in the future, it will be possible to specify arbitrary criteria, such as the proximity to a given point or set of coordinates.

Computing a Dendrogram¶

Dendrograms can be computed from an n-dimensional array using:

>>> from astrodendro import Dendrogram
>>> d = Dendrogram.compute(array)

where array is a Numpy array and d is then an instance of the Dendrogram class, which can be used to access the computed dendrogram (see Exploring the Dendrogram below). Where the array comes from is not important - for example it can be read in from a FITS file, from an HDF5 file, or it can be generated in memory. If you are interested in making a dendrogram from data in a FITS file, you can do:

>>> from astropy.io import fits
>>> array = fits.getdata('observations.fits')
>>> from astrodendro import Dendrogram
>>> d = Dendrogram.compute(array)

• min_value: the minimum value to consider in the dataset - any value lower than this will not be considered in the dendrogram. If you are working with observations, it is likely that you will want to set this to the “detection level,” for example 3- or 5-sigma, so that only significant values are included in the dendrogram. By default, all values are used.
• min_delta: how significant a leaf has to be in order to be considered an independent entity. The significance is measured from the difference between its peak flux and the value at which it is being merged into the tree. If you are working with observational data, then you could set this to, e.g., 1-sigma, which means that any leaf that is locally less than 1-sigma tall is combined with its neighboring leaf or branch and is no longer considered a separate entity.
• min_npix: the minimum number of pixels/values needed for a leaf to be considered an independent entity. When the dendrogram is being computed, and when a leaf is about to be joined onto a branch or another leaf, if the leaf has fewer than this number of pixels, then it is combined with the branch or leaf it is being merged with and is no longer considered a separate entity. By default, this parameter is set to zero, so there is no minimum number of pixels required for leaves to remain independent entities.

These options are illustrated graphically in An Illustrated Description of the Core Algorithm.

As an example, we can use a publicly available extinction map of the Perseus star-formation region from the The COordinated Molecular Probe Line Extinction Thermal Emission (COMPLETE) Survey of Star Forming Regions (PerA_Extn2MASS_F_Gal.fits, originally obtained from http://hdl.handle.net/10904/10080). The units of the map are magnitudes of extinction, and we want to make a dendrogram of all structures above a minimum value of 2 magnitudes, and we only consider leaves with at least 10 pixels and which have a peak to base difference larger than one magnitude of extinction:

>>> from astrodendro import Dendrogram
>>> from astropy.io import fits
>>> image = fits.getdata('PerA_Extn2MASS_F_Gal.fits')
>>> d = Dendrogram.compute(image, min_value=2.0, min_delta=1., min_npix=10)

By default, the computation will be silent, but for large dendrograms, it can be useful to have an idea of how long the computation will take:

>>> d = Dendrogram.compute(image, min_value=2.0, min_delta=1., min_npix=10,
                           verbose=True)
Generating dendrogram using 6,386 of 67,921 pixels (9% of data)
[=========================>               ] 64%

The ‘9% of data’ indicates that only 9% of the data are over the min_value threshold.

Exploring the Dendrogram¶

Once the dendrogram has been computed, you will want to explore/visualize it. You can access the full tree from the computed dendrogram. Assuming that you have computed a dendrogram with:

>>> d = Dendrogram.compute(array, ...)

you can now access the full tree from the d variable.

The first place to start is the trunk of the tree (the trunk attribute), which is a list of all the structures at the lowest level. Unless you left min_value to the default setting, which would mean that all values in the dataset are used, it’s likely that not all structures are connected. So the trunk is a collection of items at the lowest level, each of which could be a leaf or a branch (itself having leaves and branches). In the case of the Perseus extinction map, we get:

>>> d.trunk
[<Structure type=leaf idx=101>,
 <Structure type=branch idx=2152>,
 <Structure type=leaf idx=733>,
 <Structure type=branch idx=303>]

In the above case, the trunk contains two leaves and two branches. Since trunk is just a list, you can access items in it with e.g.:

>>> d.trunk[1]
<Structure type=branch idx=2152>

Branches have a children attribute that returns a list of all sub-structures, which can include branches and leaves. Thus, we can return the sub-structures of the above branch with:

>>> d.trunk[1].children
[<Structure type=branch idx=1680>,
 <Structure type=branch idx=5771>]

which shows that the child branch is composed of two more branches. We can therefore access the sub-structures of these branch with e.g.:

>>> d.trunk[1].children[0].children
[<Structure type=leaf idx=1748>,
 <Structure type=leaf idx=1842>]

which shows this branch splitting into two leaves.

We can access the properties of leaves as follows:

>>> leaf = d.trunk[1].children[0].children[0]
>>> leaf.indices
(array([143, 142, 142, 142, 139, 141, 141, 141, 143, 140, 140]),
 array([116, 114, 115, 116, 115, 114, 115, 116, 115, 115, 114]))
>>> leaf.values
array([ 2.7043395 ,  2.57071948,  3.4551146 ,  3.29953575,  2.53844047,
        2.59633183,  3.11309052,  2.70936489,  2.81024122,  2.76864815,
        2.52840114], dtype=float32)

A full list of attributes and methods for leaves and branches (i.e. structures) is available from the Structure page, while a list of attributes and methods for the dendrogram itself is available from the Dendrogram page.

Saving and loading the dendrogram¶

A Dendrogram object can be exported to an HDF5 file (requires h5py) or FITS file (requires astropy). To export the dendrogram to a file, use:

>>> d.save_to('my_dendrogram.hdf5')

or:

>>> d.save_to('my_dendrogram.fits')

and to load and existing dendrogram:

>>> d = Dendrogram.load_from('my_other_dendrogram.hdf5')

or:

>>> d = Dendrogram.load_from('my_other_dendrogram.fits')

If you wish, you can use this to separate the computation and analysis of the dendrogram into two scripts, to ensure that the dendrogram is only computed once. For example, you could have a script compute.py that contains:

from astropy.io import fits
from astrodendro import Dendrogram

array = fits.getdata('observations.fits')
d = Dendrogram.compute(array)
d.save_to('dendrogram.fits')

and a second file containing:

from astrodendro import Dendrogram
d = Dendrogram.load_from('dendrogram.fits')

# any analysis code here

Interactive Visualization¶

One you have computed your dendrogram, the easiest way to view it interactively is to use the viewer() method:

d = Dendrogram.compute(...)
v = d.viewer()
v.show()

This will launch an interactive window showing the original data, and the dendrogram itself. Note that the viewer is only available for 2 or 3-d datasets. The main window will look like this:

Within the viewer, you can:

Highlight structures: either click on structures in the dendrogram to highlight them, which will also show them in the image view on the left, or click on pixels in the image and have the corresponding structure be highlighted in the dendrogram plot. Clicking on a branch in the dendrogram plot or in the image will highlight that branch and all sub-structures.

Multiple structures can be highlighted in different colors using the three mouse buttons: Mouse button 1 (Left-click or “regular” click), button 2 (Middle-click or “alt+click”), and button 3 (Right-click/”ctrl+click”). Each selection is independent of the other two; any of the three can be selected either by clicking on the image or the dendrogram.

Change the image stretch: use the vmin and vmax sliders above the image to change the lower and upper level of the image stretch.

Change slice in a 3-d cube: if you select a structure in the dendrogram for a 3-d cube, the cube will automatically change to the slice containing the peak pixel of the structure (including sub-structures). However, you can also change slice manually by using the slice slider.

View the selected structure ID: in a computed dendrogram, every structure has a unique integer ID (the .idx attribute) that can be used to recognize the identify the structure when computing catalogs or making plots manually (see below).

Display astronomical coordinates: If your data has an associated WCS object (for example, if you loaded your data from a FITS file with astronomical coordinate information), the interactive viewer will display the coordinates using wcsaxes:

from astropy.io.fits import getdata
from astropy import wcs

data, header = getdata('astrodendro/docs/PerA_Extn2MASS_F_Gal.fits', header=True)
wcs = wcs.WCS(header)
d = astrodendro.Dendrogram.compute(data, wcs=wcs)
v = d.viewer()
v.show()

Note that this functionality requires that the wcsaxes package is installed. Installation instructions can be found here: http://wcsaxes.readthedocs.org/en/latest/

Linked scatter plots: If you have built a catalog (see Computing Dendrogram Statistics), you can also display a scatterplot of two catalog columns, linked to the viewer. The available catalog columns can be accessed as `catalog.colnames`. Selections in the main viewer update the colors of the points in this plot:

from astrodendro.scatter import Scatter
... code to create a dendrogram (d) and catalog ...
dv = d.viewer()
ds = Scatter(d, dv.hub, catalog, 'radius', 'v_rms')
dv.show()

The catalog properties of dendrogram structures will be plotted here. You can select structures directly from the scatter plot by clicking and dragging a lasso, and the selected structures will be highlighted in other plots:

To set logarithmic scaling on either the x axis, the y axis, or both, the following convenience methods are defined:

ds.set_semilogx()
ds.set_semilogy()
ds.set_loglog()

# To unset logarithmic scaling, pass `log=False` to the above methods, i.e.
ds.set_loglog(False)

Making plots for publications¶

While the viewer is useful for exploring the dendrogram, it does not allow one to produce publication-quality plots. For this, you can use the non-interactive plotting interface. To do this, you can first use the plotter() method to provide a plotting tool:

d = Dendrogram.compute(...)
p = d.plotter()

and then use this to make the plot you need. The following complete example shows how to make a plot of the dendrogram of the extinction map of the Perseus region (introduced in Computing and exploring dendrograms) using plot_tree(), highlighting two of the main branches:

import matplotlib.pyplot as plt
from astropy.io import fits
from astrodendro import Dendrogram

image = fits.getdata('PerA_Extn2MASS_F_Gal.fits')
d = Dendrogram.compute(image, min_value=2.0, min_delta=1., min_npix=10)
p = d.plotter()

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

# Plot the whole tree
p.plot_tree(ax, color='black')

# Highlight two branches
p.plot_tree(ax, structure=8, color='red', lw=2, alpha=0.5)
p.plot_tree(ax, structure=24, color='orange', lw=2, alpha=0.5)

# Add axis labels
ax.set_xlabel("Structure")
ax.set_ylabel("Flux")

You can find out the structure ID you need either from the interactive viewer presented above, or programmatically by accessing the idx attribute of a Structure.

A plot_contour() method is also provided to outline the contours of structures. Calling plot_contour() without any arguments results in a contour corresponding to the value of min_value used being shown.

import matplotlib.pyplot as plt
from astropy.io import fits
from astrodendro import Dendrogram

image = fits.getdata('PerA_Extn2MASS_F_Gal.fits')
d = Dendrogram.compute(image, min_value=2.0, min_delta=1., min_npix=10)
p = d.plotter()

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.imshow(image, origin='lower', interpolation='nearest', cmap=plt.cm.Blues, vmax=4.0)

# Show contour for ``min_value``
p.plot_contour(ax, color='black')

# Highlight two branches
p.plot_contour(ax, structure=8, lw=3, colors='red')
p.plot_contour(ax, structure=24, lw=3, colors='orange')

Plotting contours of structures in third-party packages¶

In some cases you may want to plot the contours in third party packages such as APLpy or DS9. For these cases, the best approach is to output FITS files with a mask of the structures to plot (one mask file per contour color you want to show).

Let’s first take the plot above and make a contour plot in APLpy outlining all the leaves. We can use the get_mask() method to retrieve the footprint of a given structure:

import aplpy
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astrodendro import Dendrogram

hdu = fits.open('PerA_Extn2MASS_F_Gal.fits')[0]
d = Dendrogram.compute(hdu.data, min_value=2.0, min_delta=1., min_npix=10)

# Create empty mask. For each leaf we do an 'or' operation with the mask so
# that any pixel corresponding to a leaf is set to True.
mask = np.zeros(hdu.data.shape, dtype=bool)
for leaf in d.leaves:
    mask = mask | leaf.get_mask()

# Now we create a FITS HDU object to contain this, with the correct header
mask_hdu = fits.PrimaryHDU(mask.astype('short'), hdu.header)

# We then use APLpy to make the final plot
fig = aplpy.FITSFigure(hdu, figsize=(8, 6))
fig.show_colorscale(cmap='Blues', vmax=4.0)
fig.show_contour(mask_hdu, colors='red', linewidths=0.5)
fig.tick_labels.set_xformat('dd')
fig.tick_labels.set_yformat('dd')

Now let’s take the example from Making plots for publications and try and reproduce the same plot. As described there, one way to find interesting structures in the dendrogram is to use the Interactive Visualization tool. This tool will give the ID of a structure as an integer (idx).

Because we are starting from this ID rather than a Structure object, we need to first get the structure, which can be done with:

structure = d[idx]

where d is a Dendrogram instance. We also want to create a different mask for each contour so as to have complete control over the colors:

import aplpy
from astropy.io import fits
from astrodendro import Dendrogram

hdu = fits.open('PerA_Extn2MASS_F_Gal.fits')[0]
d = Dendrogram.compute(hdu.data, min_value=2.0, min_delta=1., min_npix=10)

# Find the structures
structure_08 = d[8]
structure_24 = d[24]

# Extract the masks
mask_08 = structure_08.get_mask()
mask_24 = structure_24.get_mask()

# Create FITS HDU objects to contain the masks
mask_hdu_08 = fits.PrimaryHDU(mask_08.astype('short'), hdu.header)
mask_hdu_24 = fits.PrimaryHDU(mask_24.astype('short'), hdu.header)

# Use APLpy to make the final plot
fig = aplpy.FITSFigure(hdu, figsize=(8, 6))
fig.show_colorscale(cmap='Blues', vmax=4.0)
fig.show_contour(hdu, levels=[2.0], colors='black', linewidths=0.5)
fig.show_contour(mask_hdu_08, colors='red', linewidths=0.5)
fig.show_contour(mask_hdu_24, colors='orange', linewidths=0.5)
fig.tick_labels.set_xformat('dd')
fig.tick_labels.set_yformat('dd')

Deriving statistics for individual structures¶

In order to derive statistics for a given structure, you will need to use the PPStatistic or the PPVStatistic classes from the astrodendro.analysis module, e.g.:

>>> from astrodendro.analysis import PPStatistic
>>> stat = PPStatistic(structure)

where structure is a Structure instance from a dendrogram. The resulting object then has methods to compute various statistics. Using the example data from Computing and exploring dendrograms:

>>> from astrodendro import Dendrogram
>>> from astropy.io import fits
>>> image = fits.getdata('PerA_Extn2MASS_F_Gal.fits')
>>> d = Dendrogram.compute(image, min_value=2.0, min_delta=1., min_npix=10)

we can get statistics for the first structure in the trunk, which is a leaf:

>>> from astrodendro.analysis import PPStatistic
>>> d.trunk[0]
<Structure type=leaf idx=101>
>>> stat = PPStatistic(d.trunk[0])
>>> stat.major_sigma  # length of major axis on the sky
<Quantity 1.882980574564531 pix>
>>> stat.minor_sigma  # length of minor axis on the sky
<Quantity 1.4639300383020182 pix>
>>> stat.position_angle  # position angle on the sky
<Quantity 134.61988014787443 deg>

>>> q = 1.882980574564531 * u.pix
>>> q.unit
Unit("pix")
>>> q.value
1.882980574564531

Specifying meta-data when computing statistics¶

In some cases, meta-data can or should be specified. To demonstrate this, we will use a different data set which is a small section (L1551_scuba_850mu.fits) of a SCUBA 850 micron map from the SCUBA Legacy Catalog. This map has a pixel scale of 6 arcseconds per pixel and a circular beam with a full-width at half maximum (FWHM) of 22.9 arcseconds. First, we compute the dendrogram as usual:

>>> from astropy.io import fits
>>> from astrodendro import Dendrogram
>>> image = fits.getdata('L1551_scuba_850mu.fits')
>>> d = Dendrogram.compute(image, min_value=0.1, min_delta=0.02)

then we set up a Python dictionary containing the required meta-data:

>>> from astropy import units as u
>>> metadata = {}
>>> metadata['data_unit'] = u.Jy / u.beam
>>> metadata['spatial_scale'] =  6 * u.arcsec
>>> metadata['beam_major'] =  22.9 * u.arcsec # FWHM
>>> metadata['beam_minor'] =  22.9 * u.arcsec # FWHM

Note that the meta-data required depends on the units of the data and whether you are working in position-position or position-position-velocity (see Required metadata).

Finally, as before, we use the PPStatistic class to extract properties for the first structure:

>>> from astrodendro.analysis import PPStatistic
>>> stat = PPStatistic(d.trunk[0], metadata=metadata)
>>> stat.major_sigma
<Quantity 20.34630778380526 arcsec>
>>> stat.minor_sigma
<Quantity 8.15504176035544 arcsec>
>>> stat.position_angle
<Quantity 85.14309012311242 deg>
>>> stat.flux
<Quantity 0.24119688679751278 Jy>

Note that the major and minor sigma on the sky of the structures are now in arcseconds since the spatial scale was specified, and the flux (density) has been converted from Jy/beam to Jy. Note also that the flux does not include any kind of background subtraction, and is just a plain sum of the values in the structure, converted to Jy

Making a catalog¶

In order to produce a catalog of properties for all structures, it is also possible to make use of the pp_catalog() and ppv_catalog() functions. We demonstrate this using the same SCUBA data as used above:

>>> from astropy.io import fits
>>> from astrodendro import Dendrogram, pp_catalog
>>> image = fits.getdata('L1551_scuba_850mu.fits')
>>> d = Dendrogram.compute(image, min_value=0.1, min_delta=0.02)

>>> from astropy import units as u
>>> metadata = {}
>>> metadata['data_unit'] = u.Jy / u.beam
>>> metadata['spatial_scale'] =  6 * u.arcsec
>>> metadata['beam_major'] =  22.9 * u.arcsec
>>> metadata['beam_minor'] =  22.9 * u.arcsec

>>> cat = pp_catalog(d, metadata)
>>> cat.pprint(show_unit=True, max_lines=10)
_idx       flux       major_sigma   minor_sigma  ...     radius        x_cen         y_cen
            Jy           arcsec        arcsec    ...     arcsec         pix           pix
---- --------------- ------------- ------------- ... ------------- ------------- -------------
   7  0.241196886798 20.3463077838 8.15504176036 ... 12.8811874315 168.053017504 3.98809714744
  51  0.132470059814 14.2778133293 4.81100492125 ...  8.2879810685  163.25495657 9.13394216473
  60 0.0799106574322 9.66298008473 3.47364264736 ... 5.79359471511 169.278409915 15.1884110291
 ...             ...           ...           ... ...           ...           ...           ...
1203  0.183438198239 22.7202518034 4.04690367115 ... 9.58888264776 15.3760934458 100.136384362
1384   2.06217635837 38.1060171889  19.766115194 ... 27.4446338168 136.429313911 107.190835447
1504   1.90767291972 8.64476839751 8.09070477357 ... 8.36314946298  68.818705665 120.246719845

The catalog function returns an Astropy Table object.

Note that pp_catalog() and ppv_catalog() generate warnings if required meta-data is missing and sensible defaults can be assumed. If no sensible defaults can be assumed (e.g. for data_unit) then an exception is raised.

Unlike clumpfind-style algorithms, there is not a one-to-one mapping between identifiers and pixels in the map: each pixel may belong to multiple nested branches in the catalog.

Available statistics¶

For a full list of available statistics for each type of statistic class, see PPStatistic and PPVStatistic. For more information on the quantities listed in these pages, consult the paper on Bias Free Measurements of Molecular Cloud Properties or the original dendrogram paper. In the terminology of the dendrogram paper, the quantities in PPStatistic and PPVStatistic adopt the “bijection” paradigm.

Required metadata¶

As described above, the metadata needed by the statistic routines depends on what statistics are required and on the units of the data. With the exception of wcs, all meta-data should be specified as Astropy Quantity objects (e.g., 3 * u.arcsec):

• data_unit is required in order to compute the flux, so it is needed for both the pp_catalog() and ppv_catalog() functions, as well as for the flux attribute of the PPStatistic and PPVStatistic classes. Note: if data_unit is given as K, it is interpreted as units of main beam brightness temperature, following the conventions in the astropy Brightness Temperature / Flux Density equivalency .
• spatial_scale is required if the data are in units of surface brightness (e.g. MJy/sr, Jy/beam, or K) so as to be able to convert the surface brightness to the flux in each pixel. Even if the data are not in units of surface brightness, the spatial_scale can optionally be specified, causing any derived size (e.g. major_sigma) to be in the correct units instead of in pixels.
• velocity_scale can optionally be specified for PPV data, causing v_rms to be in the correct units instead of in pixels.
• beam_major and beam_minor are required if the data units depend on the beam (e.g. Jy/beam or K).
• vaxis can optionally be specified when using 3-dimensional data to indicate which dimension corresponds to the velocity. By default, this is 0, which corresponds to the third axis in a FITS file (because the dimensions are reversed in numpy).
• wavelength is required if the data are in monochromatic flux densities per unit wavelength because the fluxes need to be converted to monochromatic flux densities per unit frequency. It is also required if the data are in brightness temperature units of K.
• wcs can optionally be specified and should be a WCS instance. If specified, it allows x_cen, y_cen, and v_cen to be in the correct world coordinates rather than in pixel coordinates.

Example¶

The following example shows how to combine the plotting functionality in Plotting Dendrograms and the analysis tools shown above, to overlay ellipses approximating the structures on top of the structures themselves:

from astropy.io import fits

from astrodendro import Dendrogram
from astrodendro.analysis import PPStatistic

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

hdu = fits.open('PerA_Extn2MASS_F_Gal.fits')[0]

d = Dendrogram.compute(hdu.data, min_value=2.0, min_delta=1., min_npix=10)
p = d.plotter()

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

ax.imshow(hdu.data, origin='lower', interpolation='nearest',
          cmap=plt.cm.Blues, vmax=6.0)

for leaf in d.leaves:

    p.plot_contour(ax, structure=leaf, lw=3, colors='red')

    s = PPStatistic(leaf)
    ellipse = s.to_mpl_ellipse(edgecolor='orange', facecolor='none')

    ax.add_patch(ellipse)

ax.set_xlim(75., 170.)
ax.set_ylim(120., 260.)

As shown above, the PPStatistic and PPVStatistic classes have a to_mpl_ellipse() method to convert the first and second moments of the structures into schematic ellipses.

Specifying a custom structure merging strategy¶

By default, the decision about whether a leaf remains independent (i.e., whether it remains a leaf or its pixels get incorporated into another branch) when merged is made based on the min_delta and min_npix parameters, but in some cases, you may want to use more specialized criteria. For example, you may want only leaves overlapping with a certain position, or you may want leaves with a certain spatial or velocity extent, or a minimum peak value, to be considered independent structures.

In order to accomodate this, the compute() method can optionally take an is_independent argument which should be a function with the following call signature:

def is_independent(structure, index=None, value=None):
    ...

where structure is the Structure object that is being considered, and index and value are the pixel index and value of the pixel that is linking the structure to the rest of the tree. These last two values are only set when calling the is_independent function during the tree computation, but the is_independent function is also used at the end of the computation to prune leaves that are not attached to the tree, and in this case index and value are not set.

The following example compares the dendrogram obtained with and without a custom is_independent function:

import matplotlib.pyplot as plt
from astropy.io import fits
from astrodendro import Dendrogram

image = fits.getdata('PerA_Extn2MASS_F_Gal.fits')

fig = plt.figure(figsize=(15,5))

# Default merging strategy

d1 = Dendrogram.compute(image, min_value=2.0)
p1 = d1.plotter()

ax1 = fig.add_subplot(1, 3, 1)
p1.plot_tree(ax1, color='black')
ax1.hlines(3.5, *ax1.get_xlim(), color='b', linestyle='--')
ax1.set_xlabel("Structure")
ax1.set_ylabel("Flux")
ax1.set_title("Default merging")

# Require minimum peak value
# this is equivalent to
# custom_independent = astrodendro.pruning.min_peak(3.5)
def custom_independent(structure, index=None, value=None):
    peak_index, peak_value = structure.get_peak()
    return peak_value > 3.5

d2 = Dendrogram.compute(image, min_value=2.0,
                        is_independent=custom_independent)
p2 = d2.plotter()

ax2 = fig.add_subplot(1, 3, 2)
p2.plot_tree(ax2, color='black')
ax2.hlines(3.5, *ax2.get_xlim(), color='b', linestyle='--')
ax2.set_xlabel("Structure")
ax2.set_ylabel("Flux")
ax2.set_title("Custom merging")

# For comparison, this is what changing the min_value does:
d3 = Dendrogram.compute(image, min_value=3.5)
p3 = d3.plotter()

ax3 = fig.add_subplot(1, 3, 3)
p3.plot_tree(ax3, color='black')
ax3.hlines(3.5, *ax3.get_xlim(), color='b', linestyle='--')
ax3.set_xlabel("Structure")
ax3.set_ylabel("Flux")
ax3.set_title("min_value=3.5 merging")
ax3.set_ylim(*ax2.get_ylim())

Several pre-implemented functions suitable for use as is_independent tests are provided in astrodendro.pruning. In addition, the astrodendro.pruning.all_true() function can be used to combine several criteria. For example, the following code builds a dendrogram where each leaf contains a pixel whose value >=20, and whose pixels sum to >= 100:

from astrodendro.pruning import all_true, min_peak, min_sum

custom_independent = all_true((min_peak(20), min_sum(100)))
Dendrogram.compute(image, is_independent=custom_independent)

Handling custom adjacency logic¶

By default, neighbours to a given pixel are considered to be the adjacent pixels in the array. However, not all data are like this. For example, all-sky cartesian maps are periodic along the X axis.

You can specify custom neighbour logic by providing a neighbours function to Dendrogram.compute(). For example, the periodic_neighbours() utility will wrap neighbours across array edges. To correctly compute dendrograms for all-sky Cartesian maps:

periodic_axis = 1  # data wraps along longitude axis
Dendrogram.compute(data, neighbours=periodic_neighbours(periodic_axis))

Computing a dendrogram¶

Rather than computing a dendrogram using:

d = Dendrogram(data)
d.compute(...)

you should now use:

d = Dendrogram.compute(data)

In addition, the following options for compute have been renamed:

• minimum_flux is now min_value (since we expect dendrograms to be used not only for images, but also e.g. density fields).
• minimum_delta is now min_delta
• minimum_npix is now min_npix

Dendrogram methods and attributes¶

The following dendrogram methods have changed:

• get_leaves() has now been replaced by a leaves attribute (it is no longer a method.)
• the to_hdf5() and from_hdf5() methods have been replaced by save_to()

Leaf and Branch classes¶

The Leaf and Branch classes no longer exist, and have been replaced by a single Structure class that instead has is_leaf and is_branch attributes. Thus, if you were checking if something was a leaf by doing e.g.:

if type(s) == Leaf:
    # code here

or:

if isinstance(s, Leaf):
     # code here

then you will instead need to use:

if s.is_leaf:
     # code here

Leaf and branch attributes¶

The following leaf and branch attributes have changed:

• f has been replaced by a method called values() that can take a subtree= option that indicates whether pixels in sub-structures should be included.
• coords has been replaced by a method called indices() that can take a subtree= option that indicates whether pixels in sub-structures should be included.
• height now has a different definition - it is vmax for a leaf, or the smallest vmin of the children for a branch - this is used when plotting the dendrogram, to know at what height to plot the structure.

Interactive visualization¶

Visualizing the results of the dendrogram is now much easier, and does not require the additional astrocube package. To launch the interactive viewer (which requires only Matplotlib), once the dendrogram has been computed, you can do:

>>> d.viewer()

and the interactive viewer will launch. It will however no longer have the option to re-compute the dendrogram from the window, and will also no longer have an IPython terminal. For the latter, we recommend you consider using the Glue package.