setPrimaryFlags.py

Go to the documentation of this file.

#!/usr/bin/env python
#
# LSST Data Management System
# Copyright 2008-2016 LSST/AURA
#
# This product includes software developed by the
# LSST Project (http://www.lsst.org/).
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.    See the
# GNU General Public License for more details.
#
# You should have received a copy of the LSST License Statement and
# the GNU General Public License along with this program.  If not,
# see <http://www.lsstcorp.org/LegalNotices/>.
#
import numpy as np
from lsst.pex.config import Config, Field, ListField
from lsst.pipe.base import Task
from lsst.geom import Box2D
 
 
def getPatchInner(sources, patchInfo):
    """Set a flag for each source if it is in the innerBBox of a patch.
 
    Parameters
    ----------
    sources : `lsst.afw.table.SourceCatalog`
        A sourceCatalog with pre-calculated centroids.
    patchInfo : `lsst.skymap.PatchInfo`
        Information about a `SkyMap` `Patch`.
 
    Returns
    --------
    isPatchInner : array-like of `bool`
        `True` for each source that has a centroid
        in the inner region of a patch.
    """
    # Extract the centroid position for all the sources
    x = sources["slot_Centroid_x"]
    y = sources["slot_Centroid_y"]
    centroidFlag = sources["slot_Centroid_flag"]
 
    # set inner flags for each source and set primary flags for
    # sources with no children (or all sources if deblend info not available)
    innerFloatBBox = Box2D(patchInfo.getInnerBBox())
    inInner = innerFloatBBox.contains(x, y)
 
    # When the centroider fails, we can still fall back to the peak,
    # but we don't trust that quite as much -
    # so we use a slightly smaller box for the patch comparison.
    shrunkInnerFloatBBox = Box2D(innerFloatBBox)
    shrunkInnerFloatBBox.grow(-1)
    inShrunkInner = shrunkInnerFloatBBox.contains(x, y)
 
    # Flag sources contained in the inner region of a patch
    isPatchInner = (centroidFlag & inShrunkInner) | (~centroidFlag & inInner)
    return isPatchInner
 
 
def getTractInner(sources, tractInfo, skyMap):
    """Set a flag for each source that the skyMap includes in tractInfo.
 
    Parameters
    ----------
    sources : `lsst.afw.table.SourceCatalog`
        A sourceCatalog with pre-calculated centroids.
    tractInfo : `lsst.skymap.TractInfo`
        Tract object
    skyMap : `lsst.skymap.BaseSkyMap`
        Sky tessellation object
 
    Returns
    -------
    isTractInner : array-like of `bool`
        True if the skyMap.findTract method returns
        the same tract as tractInfo.
    """
    tractId = tractInfo.getId()
    isTractInner = np.array([skyMap.findTract(s.getCoord()).getId() == tractId for s in sources])
    return isTractInner
 
 
def getPseudoSources(sources, pseudoFilterList, schema, log):
    """Get a flag that marks pseudo sources.
 
    Some categories of sources, for example sky objects,
    are not really detected sources and should not be considered primary
    sources.
 
    Parameters
    ----------
    sources : `lsst.afw.table.SourceCatalog`
        The catalog of sources for which to identify "pseudo"
        (e.g. sky) objects.
    pseudoFilterList : `list` of `str`
        Names of filters which should never be primary
 
    Returns
    -------
    isPseudo : array-like of `bool`
        True for each source that is a pseudo source.
        Note: to remove pseudo sources use `~isPseudo`.
    """
    # Filter out sources that should never be primary
    isPseudo = np.zeros(len(sources), dtype=bool)
    for filt in pseudoFilterList:
        try:
            pseudoFilterKey = schema.find("merge_peak_%s" % filt).getKey()
            isPseudo |= sources[pseudoFilterKey]
        except KeyError:
            log.warning("merge_peak is not set for pseudo-filter %s", filt)
    return isPseudo
 
 
def getDeblendPrimaryFlags(sources):
    """Get flags generated by the deblender
 
    scarlet is different than meas_deblender in that it is not
    (necessarily) flux conserving. For consistency in scarlet,
    all of the parents with only a single child (isolated sources)
    need to be deblended. This creates a question: which type
    of isolated source should we make measurements on, the
    undeblended "parent" or the deblended child?
    For that reason we distinguish between a DeblendedSource,
    which is a source that has no children and uses the
    isolated parents, and a DeblendedModelSource, which uses
    the scarlet models for both isolated and blended sources.
    In the case of meas_deblender, DeblendedModelSource is
    `None` because it is not contained in the output catalog.
 
    Parameters
    ----------
    sources : `lsst.afw.table.SourceCatalog`
        A sourceCatalog that has already been deblended using
        either meas_extensions_scarlet or meas_deblender.
 
    Returns
    -------
    fromBlend : array-like of `bool`
        True for each source modeled by the deblender from a `Peak`
        in a parent footprint that contained at least one other `Peak`.
        While these models can be approximated as isolated,
        and measurements are made on them as if that's the case,
        we know deblending to introduce biases in the shape and centroid
        of objects and it is important to know that the sources that these
        models are based on are all bleneded in the true image.
    isIsolated : array-like of `bool`
        True for isolated sources, regardless of whether or not they
        were modeled by the deblender.
    isDeblendedSource : array-like of `bool`
        True for each source that is a "DeblendedSource" as defined above.
    isDeblendedModelSource : array-like of `bool`
        True for each source that is a "DeblendedSourceModel"
        as defined above.
    """
    nChildKey = "deblend_nChild"
    nChild = sources[nChildKey]
    parent = sources["parent"]
 
    if "deblend_scarletFlux" in sources.schema:
        # The number of peaks in the sources footprint.
        # This (may be) different than nChild,
        # the number of deblended sources in the catalog,
        # because some peaks might have been culled during deblending.
        nPeaks = sources["deblend_nPeaks"]
        parentNChild = sources["deblend_parentNChild"]
        # It is possible for a catalog to contain a hierarchy of sources,
        # so we mark the leaves (end nodes of the hierarchy tree with no
        # children).
        isLeaf = nPeaks == 1
        fromBlend = parentNChild > 1
        isIsolated = isLeaf & ((parent == 0) | parentNChild == 1)
        isDeblendedSource = (fromBlend & isLeaf) | (isIsolated & (parent == 0))
        isDeblendedModelSource = (parent != 0) & isLeaf
    else:
        # Set the flags for meas_deblender
        fromBlend = parent != 0
        isIsolated = (nChild == 0) & (parent == 0)
        isDeblendedSource = nChild == 0
        isDeblendedModelSource = None
    return fromBlend, isIsolated, isDeblendedSource, isDeblendedModelSource
 
 
class SetPrimaryFlagsConfig(Config):
    nChildKeyName = Field(dtype=str, default="deprecated",
                          doc="Deprecated. This parameter is not used.")
    pseudoFilterList = ListField(dtype=str, default=['sky'],
                                 doc="Names of filters which should never be primary")
 
 
class SetPrimaryFlagsTask(Task):
    """Add isPrimaryKey to a given schema.
 
    Parameters
    ----------
    schema : `lsst.afw.table.Schema`
        The input schema.
    isSingleFrame : `bool`
        Flag specifying if task is operating with single frame imaging.
    includeDeblend : `bool`
        Include deblend information in isPrimary and
        add isDeblendedSource field?
    kwargs :
        Keyword arguments passed to the task.
    """
 
    ConfigClass = SetPrimaryFlagsConfig
 
    def __init__(self, schema, isSingleFrame=False, **kwargs):
        Task.__init__(self, **kwargs)
        self.schemaschema = schema
        self.isSingleFrameisSingleFrame = isSingleFrame
        self.includeDeblendincludeDeblend = False
        if not self.isSingleFrameisSingleFrame:
            primaryDoc = ("true if source has no children and is in the inner region of a coadd patch "
                          "and is in the inner region of a coadd tract "
                          "and is not \"detected\" in a pseudo-filter (see config.pseudoFilterList)")
            self.isPatchInnerKeyisPatchInnerKey = self.schemaschema.addField(
                "detect_isPatchInner", type="Flag",
                doc="true if source is in the inner region of a coadd patch",
            )
            self.isTractInnerKeyisTractInnerKey = self.schemaschema.addField(
                "detect_isTractInner", type="Flag",
                doc="true if source is in the inner region of a coadd tract",
            )
        else:
            primaryDoc = "true if source has no children and is not a sky source"
        self.isPrimaryKeyisPrimaryKey = self.schemaschema.addField(
            "detect_isPrimary", type="Flag",
            doc=primaryDoc,
        )
 
        if "deblend_nChild" in schema.getNames():
            self.includeDeblendincludeDeblend = True
            self.isDeblendedSourceKeyisDeblendedSourceKey = self.schemaschema.addField(
                "detect_isDeblendedSource", type="Flag",
                doc=primaryDoc + " and is either an unblended isolated source or a "
                                 "deblended child from a parent with 'deblend_nChild' > 1")
            self.fromBlendKeyfromBlendKey = self.schemaschema.addField(
                "detect_fromBlend", type="Flag",
                doc="This source is deblended from a parent with more than one child."
            )
            self.isIsolatedKeyisIsolatedKey = self.schemaschema.addField(
                "detect_isIsolated", type="Flag",
                doc="This source is not a part of a blend."
            )
            if "deblend_scarletFlux" in schema.getNames():
                self.isDeblendedModelKeyisDeblendedModelKey = self.schemaschema.addField(
                    "detect_isDeblendedModelSource", type="Flag",
                    doc=primaryDoc + " and is a deblended child")
            else:
                self.isDeblendedModelKeyisDeblendedModelKey = None
 
    def run(self, sources, skyMap=None, tractInfo=None, patchInfo=None):
        """Set isPrimary and related flags on sources.
 
        For coadded imaging, the `isPrimary` flag returns True when an object
        has no children, is in the inner region of a coadd patch, is in the
        inner region of a coadd trach, and is not detected in a pseudo-filter
        (e.g., a sky_object).
        For single frame imaging, the isPrimary flag returns True when a
        source has no children and is not a sky source.
 
        Parameters
        ----------
        sources : `lsst.afw.table.SourceCatalog`
            A sourceTable. Reads in centroid fields and an nChild field.
            Writes is-patch-inner, is-tract-inner, and is-primary flags.
        skyMap : `lsst.skymap.BaseSkyMap`
            Sky tessellation object
        tractInfo : `lsst.skymap.TractInfo`
            Tract object
        patchInfo : `lsst.skymap.PatchInfo`
            Patch object
        """
        # Mark whether sources are contained within the inner regions of the
        # given tract/patch and are not "pseudo" (e.g. sky) sources.
        if not self.isSingleFrameisSingleFrame:
            isPatchInner = getPatchInner(sources, patchInfo)
            isTractInner = getTractInner(sources, tractInfo, skyMap)
            isPseudo = getPseudoSources(sources, self.config.pseudoFilterList, self.schemaschema, self.log)
            isPrimary = isTractInner & isPatchInner & ~isPseudo
 
            sources[self.isPatchInnerKeyisPatchInnerKey] = isPatchInner
            sources[self.isTractInnerKeyisTractInnerKey] = isTractInner
        else:
            # Mark all of the sky sources in SingleFrame images
            # (if they were added)
            if "sky_source" in sources.schema:
                isSky = sources["sky_source"]
            else:
                isSky = np.zeros(len(sources), dtype=bool)
            isPrimary = ~isSky
 
        if self.includeDeblendincludeDeblend:
            result = getDeblendPrimaryFlags(sources)
            fromBlend, isIsolated, isDeblendedSource, isDeblendedModelSource = result
            sources[self.fromBlendKeyfromBlendKey] = fromBlend
            sources[self.isIsolatedKeyisIsolatedKey] = isIsolated
            sources[self.isDeblendedSourceKeyisDeblendedSourceKey] = isDeblendedSource
            if self.isDeblendedModelKeyisDeblendedModelKey is not None:
                sources[self.isDeblendedModelKeyisDeblendedModelKey] = isDeblendedModelSource
            isPrimary = isPrimary & isDeblendedSource
 
        sources[self.isPrimaryKeyisPrimaryKey] = isPrimary