Skip to content

Commit

Permalink
photometry/astrometry can be run independently, with passing tests.
Browse files Browse the repository at this point in the history 
Split photo/astro code with Config parameters, and a test for each option in the
lsstSim test class. The tests check for empty Wcs/Calib, if they shouldn't have
been written.

Pulled apart _testJointcalTask so I could use the run/plot bits in the "only
one" tests, since those have quite different assert conditions.

Also support the split in JointcalStatistics: the nested ifs are a bit
excessive, but they get the job done.

Fixed inconsistent capitalization of jointcal in tests.
  • Loading branch information
@parejkoj
parejkoj committed Jan 24, 2017
1 parent 5dcb5c4 commit 3d1afc3
Show file tree
Hide file tree
Showing 7 changed files with 261 additions and 111 deletions.
72 changes: 45 additions & 27 deletions python/lsst/jointcal/jointcal.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,9 @@

__all__ = ["JointcalConfig", "JointcalTask"]

Photometry = collections.namedtuple('Photometry', ('fit', 'model'))
Astrometry = collections.namedtuple('Astrometry', ('fit', 'model', 'sky_to_tan_projection'))


class JointcalRunner(pipeBase.TaskRunner):
"""Subclass of TaskRunner for jointcalTask (copied from the HSC MosaicRunner)
Expand Down Expand Up @@ -76,6 +79,16 @@ def __call__(self, args):
class JointcalConfig(pexConfig.Config):
"""Config for jointcalTask"""

doAstrometry = pexConfig.Field(
doc="Fit astrometry and write the fitted result.",
dtype=bool,
default=True
)
doPhotometry = pexConfig.Field(
doc="Fit photometry and write the fitted result.",
dtype=bool,
default=True
)
coaddName = pexConfig.Field(
doc="Type of coadd, typically deep or goodSeeing",
dtype=str,
Expand Down Expand Up @@ -139,7 +152,7 @@ def _makeArgumentParser(cls):
parser = pipeBase.ArgumentParser(name=cls._DefaultName)
parser.add_argument("--profile_jointcal", default=False, action="store_true",
help="Profile steps of jointcal separately.")
parser.add_id_argument("--id", "calexp", help="data ID, e.g. --selectId visit=6789 ccd=0..9",
parser.add_id_argument("--id", "calexp", help="data ID, e.g. --id visit=6789 ccd=0..9",
ContainerClass=PerTractCcdDataIdContainer)
return parser

Expand Down Expand Up @@ -257,8 +270,6 @@ def run(self, dataRefs, profile_jointcal=False):
print("Using", filt, "band for reference flux")
refCat = loader.loadSkyCircle(center, afwGeom.Angle(radius, afwGeom.radians), filt).refCat

# associations.CollectRefStars(False) # To use USNO-A catalog

associations.CollectLSSTRefStars(refCat, filt)
associations.SelectFittedStars()
associations.DeprojectFittedStars() # required for AstromFit
Expand All @@ -271,18 +282,24 @@ def run(self, dataRefs, profile_jointcal=False):
if associations.fittedStarListSize() == 0:
raise RuntimeError('No stars in the fittedStarList!')

load_cat_prof_file = 'jointcal_fit_astrometry.prof' if profile_jointcal else ''
with pipeBase.cmdLineTask.profile(load_cat_prof_file):
astrometry = self._fit_astrometry(associations)
load_cat_prof_file = 'jointcal_fit_photometry.prof' if profile_jointcal else ''
with pipeBase.cmdLineTask.profile(load_cat_prof_file):
photometry = self._fit_photometry(associations)

# TODO: not clear that this is really needed any longer?
# TODO: makeResTuple should at least be renamed, if we do want to keep that big data-dump around.
# Fill reference and measurement n-tuples for each tract
tupleName = "res_" + str(dataRefs[0].dataId["tract"]) + ".list"
astrometry.fit.makeResTuple(tupleName)
if self.config.doAstrometry:
load_cat_prof_file = 'jointcal_fit_astrometry.prof' if profile_jointcal else ''
with pipeBase.cmdLineTask.profile(load_cat_prof_file):
astrometry = self._fit_astrometry(associations)
# TODO: not clear that this is really needed any longer?
# TODO: makeResTuple should at least be renamed, if we do want to keep that big data-dump around.
# Fill reference and measurement n-tuples for each tract
tupleName = "res_" + str(dataRefs[0].dataId["tract"]) + ".list"
astrometry.fit.makeResTuple(tupleName)
else:
astrometry = Astrometry(None, None, None)

if self.config.doPhotometry:
load_cat_prof_file = 'jointcal_fit_photometry.prof' if profile_jointcal else ''
with pipeBase.cmdLineTask.profile(load_cat_prof_file):
photometry = self._fit_photometry(associations)
else:
photometry = Photometry(None, None)

load_cat_prof_file = 'jointcal_write_results.prof' if profile_jointcal else ''
with pipeBase.cmdLineTask.profile(load_cat_prof_file):
Expand Down Expand Up @@ -322,7 +339,6 @@ def _fit_photometry(self, associations):
chi2 = fit.computeChi2()
print(chi2)

Photometry = collections.namedtuple('Photometry', ('fit', 'model'))
return Photometry(fit, model)

def _fit_astrometry(self, associations):
Expand Down Expand Up @@ -384,7 +400,6 @@ def _fit_astrometry(self, associations):
break
print("unxepected return code from minimize")

Astrometry = collections.namedtuple('Astrometry', ('fit', 'model', 'sky_to_tan_projection'))
return Astrometry(fit, model, sky_to_tan_projection)

def _write_results(self, associations, astrom_model, photom_model, visit_ccd_to_dataRef):
Expand All @@ -405,19 +420,22 @@ def _write_results(self, associations, astrom_model, photom_model, visit_ccd_to_

ccdImageList = associations.getCcdImageList()
for ccdImage in ccdImageList:
tanSip = astrom_model.ProduceSipWcs(ccdImage)
frame = ccdImage.ImageFrame()
tanWcs = afwImage.TanWcs.cast(jointcalLib.GtransfoToTanWcs(tanSip, frame, False))

# TODO: there must be a better way to identify this ccdImage?
# TODO: there must be a better way to identify this ccdImage than a visit,ccd pair?
ccd = ccdImage.getCcdId()
visit = ccdImage.getVisit()
dataRef = visit_ccd_to_dataRef[(visit, ccd)]
print("Updating WCS for visit: %d, ccd: %d" % (visit, ccd))
exp = afwImage.ExposureI(0, 0, tanWcs)
# start with the original calib saved to the ccdImage
fluxMag0, fluxMag0Sigma = ccdImage.getCalib().getFluxMag0()
exp.getCalib().setFluxMag0(fluxMag0*photom_model.photomFactor(ccdImage), fluxMag0Sigma)
exp = afwImage.ExposureI(0, 0)
if self.config.doAstrometry:
print("Updating WCS for visit: %d, ccd: %d" % (visit, ccd))
tanSip = astrom_model.ProduceSipWcs(ccdImage)
frame = ccdImage.ImageFrame()
tanWcs = afwImage.TanWcs.cast(jointcalLib.GtransfoToTanWcs(tanSip, frame, False))
exp.setWcs(tanWcs)
if self.config.doPhotometry:
print("Updating Calib for visit: %d, ccd: %d" % (visit, ccd))
# start with the original calib saved to the ccdImage
fluxMag0, fluxMag0Sigma = ccdImage.getCalib().getFluxMag0()
exp.getCalib().setFluxMag0(fluxMag0*photom_model.photomFactor(ccdImage), fluxMag0Sigma)
try:
dataRef.put(exp, 'wcs')
except pexExceptions.Exception as e:
Expand Down
146 changes: 89 additions & 57 deletions python/lsst/jointcal/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,9 @@ class JointcalStatistics(object):
diagnostic plots.
"""

def __init__(self, match_radius=0.1*arcseconds, flux_limit=100.0, verbose=False):
def __init__(self, match_radius=0.1*arcseconds, flux_limit=100.0,
do_photometry=True, do_astrometry=True,
verbose=False):
"""
Parameters
----------
Expand All @@ -44,11 +46,17 @@ def __init__(self, match_radius=0.1*arcseconds, flux_limit=100.0, verbose=False)
flux_limit : float
Signal/Noise (flux/fluxSigma) for sources to be included in the RMS cross-match.
100 is a balance between good centroids and enough sources.
verbose : bool
do_photometry : bool, optional
Perform calculations/make plots for photometric metrics.
do_astrometry : bool, optional
Perform calculations/make plots for astrometric metrics.
verbose : bool, optional
Print extra things
"""
self.match_radius = match_radius
self.flux_limit = flux_limit
self.do_photometry = do_photometry
self.do_astrometry = do_astrometry
self.verbose = verbose
self.log = lsst.log.Log.getLogger('JointcalStatistics')

Expand Down Expand Up @@ -109,27 +117,34 @@ def compute(catalogs, calibs):
new_cats = [ref.get('src') for ref in data_refs]
new_wcss = [ref.get('wcs') for ref in data_refs]
new_calibs = [wcs.getCalib() for wcs in new_wcss]
for wcs, cat in zip(new_wcss, new_cats):
# update in-place the object coordinates based on the new wcs
lsst.afw.table.utils.updateSourceCoords(wcs.getWcs(), cat)
if self.do_astrometry:
for wcs, cat in zip(new_wcss, new_cats):
# update in-place the object coordinates based on the new wcs
lsst.afw.table.utils.updateSourceCoords(wcs.getWcs(), cat)

self.new_dist, self.new_flux, self.new_ref_flux, self.new_source = compute(new_cats, new_calibs)

self._photometric_rms()

if self.verbose:
print('"photometric factor" for each data ref:')
for ref, old, new in zip(data_refs, old_calibs, new_calibs):
print(tuple(ref.dataId.values()), new.getFluxMag0()[0]/old.getFluxMag0()[0])
if self.do_photometry:
self._photometric_rms()
if self.verbose:
print('"photometric factor" for each data ref:')
for ref, old, new in zip(data_refs, old_calibs, new_calibs):
print(tuple(ref.dataId.values()), new.getFluxMag0()[0]/old.getFluxMag0()[0])
else:
self.new_PA1 = None

def rms_total(data):
"""Compute the total rms across all sources."""
total = sum(sum(dd**2) for dd in data.values())
n = sum(len(dd) for dd in data.values())
return np.sqrt(total/n)

self.old_dist_total = MatchDict(*(tuple(map(rms_total, self.old_dist))*u.radian).to(u.arcsecond))
self.new_dist_total = MatchDict(*(tuple(map(rms_total, self.new_dist))*u.radian).to(u.arcsecond))
if self.do_astrometry:
self.old_dist_total = MatchDict(*(tuple(map(rms_total, self.old_dist))*u.radian).to(u.arcsecond))
self.new_dist_total = MatchDict(*(tuple(map(rms_total, self.new_dist))*u.radian).to(u.arcsecond))
else:
self.old_dist_total = MatchDict(None, None)
self.new_dist_total = MatchDict(None, None)

Rms_result = collections.namedtuple("Rms_result", ["dist_relative", "dist_absolute", "pa1"])
return Rms_result(self.new_dist_total.relative, self.new_dist_total.absolute, self.new_PA1)
Expand Down Expand Up @@ -169,32 +184,35 @@ def make_plots(self, data_refs, old_wcs_list,
if interactive:
plt.ion()

plot_flux_distributions(plt, self.old_mag, self.new_mag,
self.old_weighted_rms, self.new_weighted_rms,
self.faint, self.bright, self.old_PA1, self.new_PA1,
name=name, outdir=outdir)
self.log.info("N data_refs: %d", len(data_refs))

if self.do_photometry:
plot_flux_distributions(plt, self.old_mag, self.new_mag,
self.old_weighted_rms, self.new_weighted_rms,
self.faint, self.bright, self.old_PA1, self.new_PA1,
name=name, outdir=outdir)

def rms_per_source(data):
"""Each element of data must already be the "delta" of whatever measurement."""
return (np.sqrt([np.mean(dd**2) for dd in data.values()])*u.radian).to(u.arcsecond)

old_dist_rms = MatchDict(*(tuple(map(rms_per_source, self.old_dist))))
new_dist_rms = MatchDict(*(tuple(map(rms_per_source, self.new_dist))))

self.log.info("N data_refs: %d", len(data_refs))
self.log.info("relative RMS (old, new): {:.2e} {:.2e}".format(self.old_dist_total.relative,
self.new_dist_total.relative))
self.log.info("absolute RMS (old, new): {:.2e} {:.2e}".format(self.old_dist_total.absolute,
self.new_dist_total.absolute))
plot_rms_histogram(plt, old_dist_rms.relative, old_dist_rms.absolute,
new_dist_rms.relative, new_dist_rms.absolute,
self.old_dist_total.relative, self.old_dist_total.absolute,
self.new_dist_total.relative, self.new_dist_total.absolute,
name=name, outdir=outdir)

plot_all_wcs_deltas(plt, data_refs, self.visits_per_dataRef, old_wcs_list,
per_ccd_plot=per_ccd_plot,
name=name, outdir=outdir)
if self.do_astrometry:
old_dist_rms = MatchDict(*(tuple(map(rms_per_source, self.old_dist))))
new_dist_rms = MatchDict(*(tuple(map(rms_per_source, self.new_dist))))

self.log.info("relative RMS (old, new): {:.2e} {:.2e}".format(self.old_dist_total.relative,
self.new_dist_total.relative))
self.log.info("absolute RMS (old, new): {:.2e} {:.2e}".format(self.old_dist_total.absolute,
self.new_dist_total.absolute))
plot_rms_histogram(plt, old_dist_rms.relative, old_dist_rms.absolute,
new_dist_rms.relative, new_dist_rms.absolute,
self.old_dist_total.relative, self.old_dist_total.absolute,
self.new_dist_total.relative, self.new_dist_total.absolute,
name=name, outdir=outdir)

plot_all_wcs_deltas(plt, data_refs, self.visits_per_dataRef, old_wcs_list,
per_ccd_plot=per_ccd_plot,
name=name, outdir=outdir)

if interactive:
plt.show()
Expand Down Expand Up @@ -281,42 +299,56 @@ def _make_match_dict(self, reference, visit_catalogs, calibs, refcalib=None):
else:
ref_flux_key = '{}_flux'

def get_fluxes(match):
"""Return (flux, ref_flux) or None if either is invalid."""
# NOTE: Protect against negative fluxes: ignore this match if we find one.
flux = match[1]['slot_CalibFlux_flux']
if flux < 0:
return None
else:
# convert to magnitudes and then Janskys, for a useable flux.
flux = fluxFromABMag(calib.getMagnitude(flux))

# NOTE: Have to protect against negative reference fluxes too.
if 'slot' in ref_flux_key:
ref_flux = match[0][ref_flux_key]
if ref_flux < 0:
return None
else:
ref_flux = fluxFromABMag(refcalib.getMagnitude(ref_flux))
else:
# a.net fluxes are already in Janskys.
ref_flux = match[0][ref_flux_key.format(filter)]
if ref_flux < 0:
return None

Flux = collections.namedtuple('Flux', ('flux', 'ref_flux'))
return Flux(flux, ref_flux)

for cat, calib, filter in zip(visit_catalogs, calibs, self.filters):
good = (cat.get('base_PsfFlux_flux')/cat.get('base_PsfFlux_fluxSigma')) > self.flux_limit
# things the classifier called sources are not extended.
good &= (cat.get('base_ClassificationExtendedness_value') == 0)
matches = lsst.afw.table.matchRaDec(reference, cat[good], self.match_radius)
for m in matches:
# NOTE: Protect against negative fluxes: ignore this match if we find one.
flux = m[1]['slot_CalibFlux_flux']
if flux < 0:
continue
else:
# convert to magnitudes and then Janskys, for a useable flux.
flux = fluxFromABMag(calib.getMagnitude(flux))

# NOTE: Have to protect against negative reference fluxes too.
if 'slot' in ref_flux_key:
ref_flux = m[0][ref_flux_key]
if ref_flux < 0:
if self.do_photometry:
flux = get_fluxes(m)
if fluxes is None:
continue
else:
ref_flux = fluxFromABMag(refcalib.getMagnitude(ref_flux))
else:
# a.net fluxes are already in Janskys.
ref_flux = m[0][ref_flux_key.format(filter)]
if ref_flux < 0:
continue
fluxes[m[0].getId()].append(flux.flux)
# we can just use assignment here, since the value is always the same.
ref_fluxes[m[0].getId()] = flux.ref_flux

if self.do_astrometry:
# Just use the computed separation distance directly.
distances[m[0].getId()].append(m[2])

# Just use the computed separation distance directly.
distances[m[0].getId()].append(m[2])
fluxes[m[0].getId()].append(flux)
# we can just use assignment here, since the value is always the same.
ref_fluxes[m[0].getId()] = ref_flux
sources[m[0].getId()].append(m[1])
# Convert to numpy array for easier math
for source in distances:
distances[source] = np.array(distances[source])
for source in fluxes:
fluxes[source] = np.array(fluxes[source])

return distances, fluxes, ref_fluxes, sources
Expand Down

0 comments on commit 3d1afc3

Please sign in to comment.