levelprops

levelprops.py

from astrodendro.analysis import ppv_catalog
import numpy as np


class TruncatedStructure(object):

    def __init__(self, parent, value):
        self._parent = parent
        self.value = value
        self.idx = parent.idx

    def _mask(self, subtree):
        result = self._parent.values(subtree)
        return result > self.value

    def values(self, subtree=True):
        result = self._parent.values(subtree)
        return result[self._mask(subtree)]

    def indices(self, subtree=True):
        result = self._parent.indices(subtree)
        mask = self._mask(subtree)
        result = tuple(r[mask] for r in result)
        return result


def _is_trunk(structure):
    return structure.parent is None


def levelprops(dendrogram, levels, metadata):
    result = None
    shp = (len(levels), len(dendrogram))
    for structure in dendrogram:
        for i, val in enumerate(levels):

            if (not _is_trunk(structure) and val < structure.vmin):
                continue

            if val > structure.height:
                continue

            s = TruncatedStructure(structure, val)
            rec = ppv_catalog([s], metadata)

            # prepare output array
            if result is None:
                dtype = np.array(rec).dtype
                result = np.zeros(shp, dtype=dtype)
                result[:] = np.nan
                for j in range(shp[1]):
                    result[:, j]['_idx'] = j

            # copy this record up through all (leafward) substructures
            rec = np.array(rec[0])
            for sid in [structure] + structure.descendants:
                result[i, sid.idx] = rec

    return result

test.py

from astrodendro import Dendrogram
from levelprops import levelprops
from astropy.io import fits
from astropy import units as u
import numpy as np


def dendro():
    data = fits.getdata('L1448.13co.un.fits')
    data = np.nan_to_num(data)
    data = np.maximum(data, 0)

    data = data[::2, ::2, ::2]
    dg = Dendrogram.compute(data, min_value=4, min_npix=15, min_delta=0.5)

    metadata = {}
    metadata['data_unit'] = u.Jy
    metadata['beam_major'] = 22.9 * u.arcsec
    metadata['beam_minor'] = 22.9 * u.arcsec
    metadata['velocity_scale'] = 1 * u.km / u.s

    print 'Dendrogram properties'
    print '#structes:', len(dg)

    levels = np.linspace(4, 8, 5)
    lp = levelprops(dg, levels, metadata)

    return dg, lp, levels


def test_indexing():
    dg, lp, levels = dendro()

    for i, s in enumerate(dg):

        actual = lp[:, s.idx]['flux']

        v = s.values()
        expected = np.array([(v * (v > l)).sum()
                             for l in levels])

        expected[expected == 0] = np.nan

        expected[s.height < levels] = np.nan
        actual[s.height < levels] = np.nan

        expected[s.vmin > levels] = np.nan
        actual[s.vmin > levels] = np.nan

        np.testing.assert_allclose(expected, actual, rtol=1e-3)


def test_monotonic():
    dg, lp, levels = dendro()
    for s in dg:
        l = lp[:, s.idx]['flux']
        l = l[np.isfinite(l)]
        np.testing.assert_array_equal(np.sort(l), l[::-1])