PPP catalog

astrodendro/__init__.py

@@ -2,7 +2,7 @@
 
 from .dendrogram import Dendrogram, periodic_neighbours
 from .structure import Structure
-from .analysis import ppv_catalog, pp_catalog
+from .analysis import ppp_catalog, ppv_catalog, pp_catalog
 from .plot import DendrogramPlotter
 from .viewer import BasicDendrogramViewer
 

astrodendro/analysis.py

@@ -17,7 +17,7 @@
 from .flux import UnitMetadataWarning
 from .progressbar import AnimatedProgressBar
 
-__all__ = ['ppv_catalog', 'pp_catalog']
+__all__ = ['ppp_catalog', 'ppv_catalog', 'pp_catalog']
 
 
 def memoize(func):
@@ -567,52 +567,210 @@ def area_exact(self):
         return self.stat.count() * dx ** 2
 
 
-class PPPStatistic(object):
+class VolumeBase(object):
 
-    def __init__(self, rhostat, vstat, metadata=None):
+    __metaclass__ = abc.ABCMeta
+
+    spatial_scale = MetadataQuantity('spatial_scale', 'Pixel width/height')
+    data_unit = MetadataQuantity('data_unit', 'Units of the pixel values', strict=True)
+
+    @abc.abstractmethod
+    def _sky_axes(self):
+        raise NotImplementedError()
+
+    def _world_pos(self):
+        xyz = self.stat.mom1()[::-1]
+        return xyz[::-1] * u.pixel
+
+
+    @abc.abstractproperty
+    def mass(self):
+        raise NotImplementedError
+
+    @abc.abstractproperty
+    def x_cen(self):
+        raise NotImplementedError()
+
+    @abc.abstractproperty
+    def y_cen(self):
+        raise NotImplementedError()
+
+    @abc.abstractproperty
+    def z_cen(self):
+        raise NotImplementedError()
+
+    @abc.abstractproperty
+    def azimuth(self):
+        raise NotImplementedError()
+
+    @abc.abstractproperty
+    def elevation(self):
+        raise NotImplementedError()
+
+    @property
+    def a_sigma(self):
         """
-        Derive properties from PPP density and velocity fields.
+        Ellispoidal semi-principal axis 'a' in the position-position-position 
+        (PPP) volume, computed from the intensity weighted second moment 
+        in direction of greatest elongation in the PPP volume.
+        """
+        dx = self.spatial_scale or u.pixel
+        a, b, c = self._sky_axes()
+        # We need to multiply the second moment by two to get the major axis
+        # rather than the half-major axis.
+        return dx * np.sqrt(self.stat.mom2_along(a))
+
+    @property
+    def b_sigma(self):
+        """
+        Ellispoidal semi-principal axis 'b' in the position-position-position 
+        (PPP) volume, computed from the intensity weighted second moment 
+        in direction of second largest elongation in the PPP volume.
+        """
+        dx = self.spatial_scale or u.pixel
+        a, b, c = self._sky_axes()
+        # We need to multiply the second moment by two to get the minor axis
+        # rather than the half-minor axis.
+        return dx * np.sqrt(self.stat.mom2_along(b))
+
+
+    @property
+    def c_sigma(self):
+        """
+        Ellispoidal semi-principal axis 'c' in the position-position-position 
+        (PPP) volume, computed from the intensity weighted second moment 
+        in direction of smallest elongation in the PPP volume.
+        """
+        dx = self.spatial_scale or u.pixel
+        a, b, c = self._sky_axes()
+        # We need to multiply the second moment by two to get the minor axis
+        # rather than the half-minor axis.
+        return dx * np.sqrt(self.stat.mom2_along(c))
+
+    @property
+    def volume_ellipsoid(self):
+        """
+        The volume of the ellipsoid defined by the second moments, where the
+        principal axes used are the HWHM (half-width at
+        half-maximum) derived from the moments.
+        """
+        return 4./3 * np.pi * self.a_sigma * self.b_sigma * self.c_sigma * (2.3548 * 0.5) ** 3
 
-        This is not currently implemented
+
+class PPPStatistic(VolumeBase):
+
+    def __init__(self, stat, metadata=None):
+        """
+        Derive properties from PPP density.
 
         Parameters
         ----------
-        rhostat : ScalarStatistic instance
-        vstat : VectorStatistic instance
+        stat : ScalarStatistic instance
         """
-        raise NotImplementedError()
+        if isinstance(stat, Structure):
+            self.stat = ScalarStatistic(stat.values(subtree=True),
+                                        stat.indices(subtree=True))
+        else:
+            self.stat = stat
+        if len(self.stat.indices) != 3:
+            raise ValueError("PPPStatistic can only be used on 3-d datasets")
+        self.metadata = metadata or {}
+
+    def _sky_axes(self):
+        ax = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
+        a, b, c = self.stat.projected_paxes(ax)
+        a = list(a)
+        b = list(b)
+        c = list(c)
+        return a, b, c
+
+
+    @property
+    def peak_density(self):
+        """
+        Maximum density in original units.
+        """
+        return np.nanmax(self.stat.values) * self.data_unit
 
     @property
     def mass(self):
-        pass
+        from .mass import compute_mass
+        return compute_mass(self.stat.mom0() * self.data_unit,
+                            u.Msun,
+                            spatial_scale=self.spatial_scale)
+
+    @property
+    def volume_exact(self):
+        """
+        The exact volume of the structure in PPP.
+        """
+        dx = self.spatial_scale or u.pixel
+        return self.stat.count() * dx**3
 
     @property
-    def volume(self):
-        pass
+    def radius(self):
+        """
+        Equivalent radius of the sphere occupying the same volume
+        of the structure.
+        """
+        return (3. * self.volume_exact / (4. * np.pi))**(1./3)
 
     @property
     def surface_area(self):
-        pass
+        """
+        Equivalent area of the circle using the equivalent
+        radius estimated in PPP.
+        """
+        return np.pi * self.radius ** 2
 
     @property
-    def virial(self):
-        pass
+    def surface_density(self):
+        """
+        Surface_density over the equivalent area.
+        """
+        return self.mass / self.surface_area
 
     @property
-    def v_rms(self):
-        pass
+    def azimuth(self):
+        """
+        The position angle of a_axis in degrees counter-clockwise
+        from the +y axis.
+        """
+        a, b, c = self._sky_axes()
+        return np.degrees(np.arctan2(a[1], a[2])) * u.degree
 
     @property
-    def vz_rms(self):
-        pass
+    def elevation(self):
+        """
+        The position angle of a_axis in degrees clockwise
+        from the +z axis.
+        """
+        a, b, c = self._sky_axes()
+        return np.degrees(np.arccos(a[0]/np.linalg.norm(a))) * u.degree
 
     @property
-    def pressure_vz(self):
-        pass
+    def x_cen(self):
+        """
+        The mean position of the structure in the x direction (in pixel
+        coordinates).
+        """
+        return self._world_pos()[2]
 
     @property
-    def pressure(self):
-        pass
+    def y_cen(self):
+        """
+        The mean position of the structure in the y direction (in pixel
+        coordinates).
+        """
+        return self._world_pos()[1]
+
+    @property
+    def z_cen(self):
+        """
+        The mean position of the structure in the z direction (in pixel
+        coordinates).
+        """
+        return self._world_pos()[0]
 
 
 def _make_catalog(structures, fields, metadata, statistic, verbose=False):
@@ -696,6 +854,37 @@ def _make_catalog(structures, fields, metadata, statistic, verbose=False):
     return result
 
 
+def ppp_catalog(structures, metadata, fields=None, verbose=True):
+    """
+    Iterate over a collection of position-position-position (PPP) structures,
+    extracting several quantities from each, and building a catalog.
+
+    Parameters
+    ----------
+    structures : iterable of Structures
+         The structures to catalog (e.g., a dendrogram)
+    metadata : dict
+        The metadata used to compute the catalog
+    fields : list of strings, optional
+        The quantities to extract. If not provided,
+        defaults to all PPV statistics
+    verbose : bool, optional
+        If True (the default), will generate warnings
+        about missing metadata
+
+    Returns
+    -------
+    table : a :class:`~astropy.table.table.Table` instance
+        The resulting catalog
+    """
+    fields = fields or ['a_sigma', 'b_sigma', 'c_sigma', 'radius', 'volume_ellipsoid', 'volume_exact',
+                        'surface_area', 'surface_density', 'azimuth', 'elevation', 
+                        'x_cen', 'y_cen', 'z_cen', 'mass', 'peak_density']
+    with warnings.catch_warnings():
+        warnings.simplefilter("once" if verbose else 'ignore', category=MissingMetadataWarning)
+        return _make_catalog(structures, fields, metadata, PPPStatistic)
+
+
 def ppv_catalog(structures, metadata, fields=None, verbose=True):
     """
     Iterate over a collection of position-position-velocity (PPV) structures,

astrodendro/mass.py

@@ -0,0 +1,94 @@
+import numpy as np
+
+from astropy import units as u
+from astropy.constants import si
+
+
+def quantity_sum(quantities):
+    """
+    In Astropy 0.3, np.sum will do the right thing for quantities, but in the mean time we need a workaround.
+    """
+    return np.sum(quantities.value) * quantities.unit
+
+
+def compute_mass(input_quantities, output_unit, spatial_scale=None):
+    """
+    Given a set of density values in Msun/pc**3 or g/cm**3 units, 
+    find the total mass in a specific set of units.
+
+    Parameters
+    ----------
+    input_quantities : `~astropy.units.quantity.Quantity` instance
+        A `~astropy.units.quantity.Quantity` instance containing an array of
+        density values to be summed.
+    output_unit : `~astropy.units.core.Unit` instance
+        The final unit to give the total mass in (should be equivalent to Msun)
+    spatial_scale : `~astropy.units.quantity.Quantity` instance
+        The pixel scale of the data
+    """
+
+    # Start off by finding the total mass in Msun
+
+    if input_quantities.unit.is_equivalent(u.Msun):
+
+        # Simply sum up the values and convert to output unit
+        total_flux = quantity_sum(input_quantities).to(u.Msun)
+
+    elif input_quantities.unit.is_equivalent(u.g / u.cm ** 3):
+
+        # Convert volume pixel to cm**3
+        if spatial_scale is None:
+            print("WARNING: spatial_scale not recognized, assumed as cm")
+            pixel_volume = u.cm ** 3          
+
+        elif spatial_scale.unit.is_equivalent(u.kpc):
+            pixel_volume = ((spatial_scale.to(u.cm)) ** 3)
+
+        elif spatial_scale.unit.is_equivalent(u.pc):        
+            pixel_volume = ((spatial_scale.to(u.cm)) ** 3)
+
+        elif spatial_scale.unit.is_equivalent(u.cm):
+            pixel_volume = (spatial_scale ** 3)
+
+        else:
+            print("WARNING: spatial_scale not recognized, assumed cm")
+            pixel_volume = u.cm ** 3            
+
+        # Find total mass in Msun
+        total_mass = quantity_sum(input_quantities) * pixel_volume# / (si.M_sun * 1e3)
+
+
+    elif input_quantities.unit.is_equivalent(u.Msun / u.pc ** 3):
+
+        # Convert volume pixel to pc**3
+        if spatial_scale is None:
+            print("WARNING: spatial_scale not recognized, assumed pc")
+            pixel_volume = u.pc ** 3
+
+        elif spatial_scale.unit.is_equivalent(u.kpc):
+            pixel_volume = ((spatial_scale.to(u.pc)) ** 3)
+
+        elif spatial_scale.unit.is_equivalent(u.pc):        
+            pixel_volume = (spatial_scale ** 3)
+
+        elif spatial_scale.unit.is_equivalent(u.cm):
+            pixel_volume = ((spatial_scale.to(u.pc)) ** 3)
+
+        else:
+            print("WARNING: spatial_scale not recognized, assumed pc")
+            pixel_volume = u.pc ** 3            
+
+        # Find total mass in Msun
+        total_mass = quantity_sum(input_quantities) * pixel_volume
+
+    else:
+
+        print("WARNING: data unit not recognized. Providing direct sum of values.")
+        output_unit = input_quantities.unit
+        total_mass = quantity_sum(input_quantities)
+        return total_mass.to(output_unit)
+
+    if not output_unit.is_equivalent(u.Msun):
+        raise ValueError("output_unit has to be equivalent to Msun")
+    else:
+        return total_mass.to(output_unit)