Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

DM-10156: Replace all uses of Calib with PhotoCalib #49

Merged
merged 2 commits into from
Apr 5, 2019
Merged

DM-10156: Replace all uses of Calib with PhotoCalib #49

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6e395ac
Replace Calib->PhotoCalib
parejkoj Mar 19, 2019
2fb35f4
Cleanup handling of colorterms
parejkoj Mar 30, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
20 changes: 10 additions & 10 deletions python/lsst/meas/mosaic/checkMosaicTask.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ def makeDiffPosFlux(self, allMat, allSource, wcsDic, calibDic, ffpDic):
ra, dec = wcs.pixelToSky(x, y).getPosition(afwGeom.degrees)
dx_m.append((ra - ra_cat) * 3600)
dy_m.append((dec - dec_cat) * 3600)
mag = 2.5*math.log10(calibDic[iexp][ichip].getFluxMag0()[0]/ss[j].getFlux())
mag = 2.5*math.log10(calibDic[iexp][ichip].getInstFluxAtZeroMagnitude()/ss[j].getFlux())
mcor = ffpDic[iexp][ichip].eval(x,y)
jcor = -2.5*math.log10(measMosaic.computeJacobian(wcs, afwGeom.Point2D(x, y)))
m0_m.append(mag_cat)
Expand Down Expand Up @@ -86,7 +86,7 @@ def makeDiffPosFlux(self, allMat, allSource, wcsDic, calibDic, ffpDic):
ra_cat += ra
dec_cat += dec

mag = 2.5*math.log10(calibDic[iexp][ichip].getFluxMag0()[0]/ss[j].getFlux())
mag = 2.5*math.log10(calibDic[iexp][ichip].getInstFluxAtZeroMagnitude()/ss[j].getFlux())
err = 2.5 / math.log(10) * ss[j].getFluxErr() / ss[j].getFlux()
mcor = ffpDic[iexp][ichip].eval(x,y)
jcor = -2.5*math.log10(measMosaic.computeJacobian(wcs, afwGeom.Point2D(x, y)))
Expand Down Expand Up @@ -131,7 +131,7 @@ def makeDiffPosFlux(self, allMat, allSource, wcsDic, calibDic, ffpDic):
ra_cat += ra
dec_cat += dec

mag = 2.5*math.log10(calibDic[iexp][ichip].getFluxMag0()[0]/ss[j].getFlux())
mag = 2.5*math.log10(calibDic[iexp][ichip].getInstFluxAtZeroMagnitude()/ss[j].getFlux())
err = 2.5 / math.log(10) * ss[j].getFluxErr() / ss[j].getFlux()
mcor = ffpDic[iexp][ichip].eval(x,y)
jcor = -2.5*math.log10(measMosaic.computeJacobian(wcs, afwGeom.Point2D(x, y)))
Expand Down Expand Up @@ -174,7 +174,7 @@ def makeFluxStat(self, allMat, allSource, calibDic, ffpDic, wcsDic):
iexp = ss[j].getExp()
ichip = ss[j].getChip()
if ss[j].getFlux() > 0.0:
mag = calibDic[iexp][ichip].getMagnitude(ss[j].getFlux())
mag = calibDic[iexp][ichip].instFluxToMagnitude(ss[j].getFlux())
err = 2.5 / math.log(10) * ss[j].getFluxErr() / ss[j].getFlux()
xs, ys = ss[j].getX(), ss[j].getY()
mcor = ffpDic[iexp][ichip].eval(xs, ys)
Expand All @@ -201,9 +201,9 @@ def makeFluxStat(self, allMat, allSource, calibDic, ffpDic, wcsDic):
for j in range(1,len(ss)):
iexp = ss[j].getExp()
ichip = ss[j].getChip()
#print iexp, ichip, calibDic[iexp][ichip].getFluxMag0()[0], ss[j].getFlux()
if calibDic[iexp][ichip].getFluxMag0()[0] > 0 and ss[j].getFlux() > 0.0:
mag = calibDic[iexp][ichip].getMagnitude(ss[j].getFlux())
#print iexp, ichip, calibDic[iexp][ichip].getInstFluxAtZeroMagnitude(), ss[j].getFlux()
if calibDic[iexp][ichip].getInstFluxAtZeroMagnitude() > 0 and ss[j].getFlux() > 0.0:
mag = calibDic[iexp][ichip].instFluxToMagnitude(ss[j].getFlux())
err = 2.5 / math.log(10) * ss[j].getFluxErr() / ss[j].getFlux()
xs, ys = ss[j].getX(), ss[j].getY()
mcor = ffpDic[iexp][ichip].eval(xs, ys)
Expand Down Expand Up @@ -411,7 +411,7 @@ def writeCatalog(self, allSource, wcsDic, calibDic, ffpDic):
ra, dec = wcs.pixelToSky(x, y).getPosition(afwGeom.degrees)
outData.ra[i] = ra
outData.dec[i] = dec
fluxMag0 = calibDic[iexp][ichip].getFluxMag0()[0]
fluxMag0 = calibDic[iexp][ichip].getInstFluxAtZeroMagnitude()
flux = src.getFlux()
if flux > 0 and fluxMag0 > 0:
mcor = ffpDic[iexp][ichip].eval(x, y)
Expand Down Expand Up @@ -462,7 +462,7 @@ def run(self, dataRefList, ct=None, debug=False, verbose=False):

md = dataRef.get('fcr_md')
ffp = measMosaic.FluxFitParams(md)
calib = afwImage.Calib(md)
photoCalib = afwImage.makePhotoCalibFromMetadata(md)

sources = dataRef.get('src',
flags=afwTable.SOURCE_IO_NO_FOOTPRINTS,
Expand Down Expand Up @@ -528,7 +528,7 @@ def run(self, dataRefList, ct=None, debug=False, verbose=False):
match[1].setChip(dataRef.dataId['ccd'])
mlVisit[dataRef.dataId['visit']].append(match)
wcsDic[dataRef.dataId['visit']][dataRef.dataId['ccd']] = wcs
calibDic[dataRef.dataId['visit']][dataRef.dataId['ccd']] = calib
calibDic[dataRef.dataId['visit']][dataRef.dataId['ccd']] = photoCalib
ffpDic[dataRef.dataId['visit']][dataRef.dataId['ccd']] = ffp
dataRefListUsed.append(dataRef)

Expand Down
43 changes: 17 additions & 26 deletions python/lsst/meas/mosaic/mosaicTask.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
import os
import math
import numpy
import astropy.units

import multiprocessing

Expand Down Expand Up @@ -392,19 +393,19 @@ def readSrc(self, dataRef):
newMatch = table.makeRecord()
newMatch.assign(match[0], mapper)
match[0] = newMatch
primaryFluxKey = refSchema.find(refSchema.join(self.cterm.primary, "flux")).key
secondaryFluxKey = refSchema.find(refSchema.join(self.cterm.secondary, "flux")).key
primaryFluxErrKey = refSchema.find(refSchema.join(self.cterm.primary, "fluxErr")).key
secondaryFluxErrKey = refSchema.find(refSchema.join(self.cterm.secondary, "fluxErr")).key
refFlux1 = numpy.array([m[0].get(primaryFluxKey) for m in matches])
refFlux2 = numpy.array([m[0].get(secondaryFluxKey) for m in matches])
refFluxErr1 = numpy.array([m[0].get(primaryFluxErrKey) for m in matches])
refFluxErr2 = numpy.array([m[0].get(secondaryFluxErrKey) for m in matches])
refMag1 = -2.5*numpy.log10(refFlux1)
refMag2 = -2.5*numpy.log10(refFlux2)
refMag = self.cterm.transformMags(refMag1, refMag2)
refFlux = numpy.power(10.0, -0.4*refMag)
refFluxErr = self.cterm.propagateFluxErrors(refFluxErr1, refFluxErr2)

# extract the matched refCat as a Catalog for the colorterm code
refCat = afwTable.SimpleCatalog(matches[0].first.schema)
refCat.reserve(len(matches))
for x in matches:
record = refCat.addNew()
record.assign(x.first)

refMag, refMagErr = self.cterm.getCorrectedMagnitudes(refCat,
afwImage.Filter(calexp_md).getName())
# NOTE: mosaic assumes fluxes are in Jy
refFlux = (refMag*astropy.units.ABmag).to_value(astropy.units.Jy)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a weird mix of astropy to compute the refFlux and the stack to compute refFluxErr.
Does afwImage.fluxFromABMag work?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm trying to get rid of all uses of that in the stack prior to deprecation. See DM-16903.

Unfortunately, astropy doesn't currently have a way to handle the error calculations directly (astropy.uncertainty isn't finished yet), otherwise I would have replaced the other one too.

refFluxErr = afwImage.fluxErrFromABMagErr(refMagErr, refMag)
matches = [self.setCatFlux(m, flux, fluxKey, fluxErr, fluxErrKey) for
m, flux, fluxErr in zip(matches, refFlux, refFluxErr) if flux == flux]
else:
Expand Down Expand Up @@ -677,9 +678,8 @@ def writeFcr(self, dataRefList):
exp = afwImage.ExposureI(0,0)
exp.getMetadata().combine(metadata)
scale = self.fexp[iexp]*self.fchip[ichip]
calib = afwImage.Calib()
calib.setFluxMag0(1.0/scale, 1.0/scale*std*M_LN10*0.4)
exp.setCalib(calib)
photoCalib = afwImage.makePhotoCalibFromCalibZeroPoint(1.0/scale, 1.0/scale*std*M_LN10*0.4)
exp.setPhotoCalib(photoCalib)
try:
dataRef.put(exp, "fcr")
except Exception as e:
Expand All @@ -700,18 +700,9 @@ def writeFcr(self, dataRefList):
except Exception as e:
print("failed to read Wcs for PhotoCalib: %s" % (e))
continue
instFluxMag0, instFluxMag0Err = calib.getFluxMag0()
# TODO: need to scale these until DM-10153 is completed and PhotoCalib has replaced Calib
referenceFlux = 1e23 * 10**(48.6 / -2.5) * 1e9
bf = measMosaic.FluxFitBoundedField(bbox, newP, wcs,
zeroPoint=instFluxMag0,
zeroPoint=photoCalib.getInstFluxAtZeroMagnitude(),
nQuarter=nQuarter)
photoCalib = afwImage.PhotoCalib(
referenceFlux / instFluxMag0,
referenceFlux * instFluxMag0Err/instFluxMag0**2,
bf,
isConstant=False
)
dataRef.put(photoCalib, "jointcal_photoCalib")

def outputDiagWcs(self):
Expand Down
4 changes: 1 addition & 3 deletions python/lsst/meas/mosaic/photometricSolution.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,9 +169,7 @@ def writeFcr(self, butler, ccdIdList, ccdSet, filterName,
exp = afwImage.ExposureI(0,0)
exp.getMetadata().combine(metadata)
scale = fexp[iexp] * fchip[ichip]
calib = afwImage.Calib()
calib.setFluxMag0(1.0/scale)
exp.setCalib(calib)
exp.setPhotoCalib(afwImage.makePhotoCalibFromCalibZeroPoint(1.0/scale))
exp.setFilter(afwImage.Filter(filterName))
try:
butler.put(exp, 'fcr', {'visit': iexp, 'ccd': ichip})
Expand Down
21 changes: 9 additions & 12 deletions python/lsst/meas/mosaic/updateExposure.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,8 @@ def applyMosaicResultsExposure(dataRef, calexp=None):
mosaic.wcs = measAstrom.rotateWcsPixelsBy90(mosaic.wcs, 4 - nQuarter, dims)
calexp.setWcs(mosaic.wcs)

calexp.getCalib().setFluxMag0(mosaic.calib.getFluxMag0())
fluxMag0 = mosaic.calib.getInstFluxAtZeroMagnitude()
calexp.setPhotoCalib(afwImage.makePhotoCalibFromCalibZeroPoint(fluxMag0, 0.0))

mi = calexp.getMaskedImage()
# rotate photometric correction to LSST coordiantes
Expand All @@ -99,7 +100,7 @@ def getFluxFitParams(dataRef):
ffpHeader = dataRef.get("fcr_hsc_md", immediate=True)
else:
ffpHeader = dataRef.get("fcr_md", immediate=True)
calib = afwImage.Calib(ffpHeader)
photoCalib = afwImage.makePhotoCalibFromMetadata(ffpHeader)
ffp = FluxFitParams(ffpHeader)

wcs = getWcs(dataRef)
Expand All @@ -112,7 +113,7 @@ def getFluxFitParams(dataRef):
import lsst.meas.astrom as measAstrom
dimensions = dataRef.get("calexp_bbox").getDimensions()
wcs = measAstrom.rotateWcsPixelsBy90(wcs, nQuarter, dimensions)
return Struct(ffp=ffp, calib=calib, wcs=wcs)
return Struct(ffp=ffp, calib=photoCalib, wcs=wcs)


def getWcs(dataRef):
Expand Down Expand Up @@ -206,7 +207,7 @@ def applyMosaicResultsCatalog(dataRef, catalog, addCorrection=True):
return Struct(catalog=catalog, wcs=wcs, ffp=ffp)


def applyCalib(catalog, calib, hscRun=None):
def applyCalib(catalog, photoCalib, hscRun=None):
"""Convert all fluxes in a catalog to magnitudes

The fluxes are converted in-place, so that the "_flux*" are now really
Expand Down Expand Up @@ -243,21 +244,17 @@ def applyCalib(catalog, calib, hscRun=None):
aliasMap.set(name, newName)
aliasMap.set(name + errName, newName + errName)

calib.setThrowOnNegativeFlux(False)

newCatalog = afwTable.SourceCatalog(mapper.getOutputSchema())
newCatalog.extend(catalog, mapper=mapper)

for name, key in newFluxKeys.items():
flux = newCatalog[key]
if name in newErrKeys:
fluxErr = newCatalog[newErrKeys[name]]
magArray = numpy.array([calib.getMagnitude(f, e) for f, e in zip(flux, fluxErr)])
mag = magArray[:, 0]
fluxErr[:] = magArray[:, 1]
result = photoCalib.instFluxToMagnitude(newCatalog, name.strip('_mag'))
flux[:] = result[:, 0]
newCatalog[newErrKeys[name]] = result[:, 1]
else:
mag = numpy.array([calib.getMagnitude(f) for f in flux])
flux[:] = mag
flux[:] = numpy.array([photoCalib.instFluxToMagnitude(f) for f in flux])

return newCatalog

Expand Down
14 changes: 7 additions & 7 deletions tests/test_fluxFitBoundedField.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -266,21 +266,21 @@ def checkPhotoCalibCatalog(self, nQuarter):
catalog2 = lsst.meas.mosaic.applyCalib(catalog2, results2.ffp.calib)
mag2 = catalog2["example_mag"]
# Check that the non-spatially varying part of the correction is the same.
fluxMag0 = results2.ffp.calib.getFluxMag0()
self.assertFloatsAlmostEqual(photoCalib.getInstFluxAtZeroMagnitude(), fluxMag0[0],
self.assertFloatsAlmostEqual(photoCalib.getCalibrationMean(),
results2.ffp.calib.getCalibrationMean(),
rtol=1E-14)
self.assertFloatsAlmostEqual(photoCalib.getCalibrationErr(), referenceFlux*fluxMag0[1]/fluxMag0[0]**2,
self.assertFloatsAlmostEqual(photoCalib.getCalibrationErr(),
results2.ffp.calib.getCalibrationErr(),
rtol=1E-14)
# Compute partially-calibrated magnitudes that don't account for the spatially-varying part.
mag0, magErr0 = results2.ffp.calib.getMagnitude(catalog.get("example_instFlux"),
catalog.get("example_instFluxErr"))
mag0 = results2.ffp.calib.instFluxToMagnitude(catalog, "example")
# Check that both approaches yield similar results overall...
rtol = 1E-10 if nQuarter == 0 else 1E-6 # rotating SIP Wcses involves a big loss of precision
self.assertFloatsAlmostEqual(mag1, mag2, rtol=rtol)
# ...and in just the spatially-varying part (but with less precision, partially because of
# round-off error).
magDiff2 = mag2 - mag0
magDiff1 = mag1 - mag0
magDiff2 = mag2 - mag0[:, 0]
magDiff1 = mag1 - mag0[:, 0]
self.assertFloatsAlmostEqual(magDiff1, magDiff2, rtol=rtol*3E3)

def checkPhotoCalibExposure(self, nQuarter):
Expand Down