compareModel.py

Go to the documentation of this file.

"""
Compare the input fake galaxy model to the output one
"""
import argparse
 
import numpy as np
import astropy.table
 
import lsst.pipe.base
import lsst.afw.table
import lsst.afw.geom.ellipses
import lsst.daf.persistence as dafPersist
 
from .matchFakes import getFakeSources
 
 
def getParams(record, galType='sersic'):
    """
    return the semi-major axis, axis ratio and position angle (in degrees)
    """
    fluxType = {'dev': 'cmodel.dev',
                'exp': 'cmodel.exp',
                'sersic': 'cmodel',
                'cmodel': 'cmodel'}[galType]
 
    e = lsst.afw.geom.ellipses.Axes(record.get(fluxType + '.ellipse'))
    q = (e.getB() / e.getA())
    reff = e.getA()
    pa = (e.getTheta() * 180.0 / np.pi)
    return reff, q, pa
 
 
def getMag(record, fluxType='cmodel'):
    """
    return the magnitude and error
    """
    flux, fluxerr = record.get(fluxType), record.get(fluxType+'.err')
    mag, magerr = -2.5 * np.log10(flux), 2.5/np.log(10.0)*fluxerr/flux
    mag += record.get('zeropoint')
    return mag, magerr
 
 
def writeNumpyTable(fakeTable):
    """Writes output to numpy
    """
    npTable = np.recarray(len(fakeTable),
                          dtype={'names': ['id', 'fakeid', 'visit', 'ccd',
                                           'cmodelMag', 'expMag',
                                           'devMag',
                                           'expReff', 'devReff',
                                           'expBA', 'devBA',
                                           'expPosAng', 'devPosAng',
                                           'cmodelMagErr', 'expMagErr',
                                           'devMagErr', 'nchild', 'parent'],
                                 'formats': [int, int, int, int,
                                             float, float,
                                             float, float, float, float,
                                             float, float, float, float,
                                             float, float, int, int]})
    for indFake, fake in enumerate(fakeTable):
        npTable[indFake]['id'] = fake.get('id')
        npTable[indFake]['fakeid'] = fake.get('fakeId')
        npTable[indFake]['ccd'] = fake.get('ccd')
        npTable[indFake]['visit'] = fake.get('visit')
        nameMatch = {'sersic': 'cmodel',
                     'exp': 'cmodel.exp',
                     'cmodel': 'cmodel',
                     'dev': 'cmodel.dev'}
        for name in ('cmodel', 'exp', 'dev'):
            m1, m2 = getMag(fake, nameMatch[name]+'.flux')
            npTable[indFake][name+'Mag'] = m1
            npTable[indFake][name+'MagErr'] = m2
        for name in ('exp', 'dev'):
            params = getParams(fake, name)
            npTable[indFake][name+'Reff'] = params[0]
            npTable[indFake][name+'BA'] = params[1]
            npTable[indFake][name+'PosAng'] = params[2]
        npTable[indFake]['nchild'] = fake.get('deblend.nchild')
        npTable[indFake]['parent'] = fake.get('parent')
    return npTable
 
 
def main(root, visit, ccds, galType='sersic', output='outputs/'):
    """
    main function for controlling fake source comparison
    """
    butler = dafPersist.Butler(root)
    fakeTable = None
    for ccd in ccds:
        dataId = {'visit': visit, 'ccd': ccd}
        try:
            temp = getFakeSources(butler, dataId, tol=1.0,
                                  visit=True, ccd=True)
        except Exception:
            continue
        if fakeTable is None:
            fakeTable = temp.copy(True)
        else:
            fakeTable.extend(temp, True)
 
    npTable = writeNumpyTable(fakeTable)
    rerunName = root.split('/')[-2]
    tableFileName = output + '/' + rerunName + '_galMags.txt'
    astropy.table.Table(npTable).write(tableFileName, format='ascii')
 
 
if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('root', help="Root directory of data repository")
    parser.add_argument("visit", type=int, help="Visit")
    parser.add_argument("--ccd", nargs='+', type=int, help="CCD(s)")
    parser.add_argument('-g', '--galtype', type=str, dest='galType',
                        choices=['exp', 'dev', 'sersic'])
    parser.add_argument('-o', '--outputpath', dest='outpath',
                        help='path for output')
    args = parser.parse_args()
 
    main(args.root, args.visit, args.ccd, galType=args.galType,
         output=args.outpath)