Skip to content

Commit

Permalink
pass mask to aperture_photometry method; apphot the mask and store it…
Browse files Browse the repository at this point in the history 
… in the tractor file
  • Loading branch information
@dstndstn
dstndstn committed Apr 29, 2020
1 parent 3adb03c commit 0e2e7f5
Show file tree
Hide file tree
Showing 3 changed files with 49 additions and 24 deletions.
41 changes: 31 additions & 10 deletions py/legacypipe/coadds.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -527,16 +527,21 @@ class Duck(object):

if apertures is not None:
# Aperture photometry
# photutils.aperture_photometry: mask=True means IGNORE
mask = (cow == 0)
with np.errstate(divide='ignore'):
imsigma = 1.0/np.sqrt(coiv)
imsigma[coiv == 0] = 0
imsigma = 1.0/np.sqrt(cow)
imsigma[mask] = 0.

for irad,rad in enumerate(apertures):
apargs.append((irad, band, rad, coimg, imsigma, True, apxy))
apargs.append((irad, band, rad, cowimg, imsigma, mask,
True, apxy))
if mods is not None:
apargs.append((irad, band, rad, coresid, None, False, apxy))
apargs.append((irad, band, rad, coresid, None, None,
False, apxy))
if blobmods is not None:
apargs.append((irad, band, rad, coblobresid, None, False, apxy))
apargs.append((irad, band, rad, coblobresid, None, None,
False, apxy))

if callback is not None:
callback(band, *callback_args, **kwargs)
Expand Down Expand Up @@ -566,39 +571,47 @@ class Duck(object):
for iband,band in enumerate(bands):
apimg = []
apimgerr = []
apmask = []
if mods is not None:
apres = []
if blobmods is not None:
apblobres = []
for irad,rad in enumerate(apertures):
(airad, aband, isimg, ap_img, ap_err) = next(apresults)
(airad, aband, isimg, ap_img, ap_err, ap_mask) = next(apresults)
assert(airad == irad)
assert(aband == band)
assert(isimg)
apimg.append(ap_img)
apimgerr.append(ap_err)
apmask.append(ap_mask)

if mods is not None:
(airad, aband, isimg, ap_img, ap_err) = next(apresults)
(airad, aband, isimg, ap_img, ap_err, ap_mask) = next(apresults)
assert(airad == irad)
assert(aband == band)
assert(not isimg)
apres.append(ap_img)
assert(ap_err is None)
assert(ap_mask is None)

if blobmods is not None:
(airad, aband, isimg, ap_img, ap_err) = next(apresults)
(airad, aband, isimg, ap_img, ap_err, ap_mask) = next(apresults)
assert(airad == irad)
assert(aband == band)
assert(not isimg)
apblobres.append(ap_img)
assert(ap_err is None)
assert(ap_mask is None)

ap = np.vstack(apimg).T
ap[np.logical_not(np.isfinite(ap))] = 0.
C.AP.set('apflux_img_%s' % band, ap)
ap = 1./(np.vstack(apimgerr).T)**2
ap[np.logical_not(np.isfinite(ap))] = 0.
C.AP.set('apflux_img_ivar_%s' % band, ap)
ap = np.vstack(apmask).T
ap[np.logical_not(np.isfinite(ap))] = 0.
C.AP.set('apflux_masked_%s' % band, ap)
if mods is not None:
ap = np.vstack(apres).T
ap[np.logical_not(np.isfinite(ap))] = 0.
Expand Down Expand Up @@ -842,16 +855,24 @@ def _resample_one(args):
return itim,Yo,Xo,iv,im,mo,bmo,dq

def _apphot_one(args):
(irad, band, rad, img, sigma, isimage, apxy) = args
(irad, band, rad, img, sigma, mask, isimage, apxy) = args
import photutils
result = [irad, band, isimage]
aper = photutils.CircularAperture(apxy, rad)
p = photutils.aperture_photometry(img, aper, error=sigma)
p = photutils.aperture_photometry(img, aper, error=sigma, mask=mask)
result.append(p.field('aperture_sum'))
if sigma is not None:
result.append(p.field('aperture_sum_err'))
else:
result.append(None)

# If a mask is passed, also photometer it!
if mask is not None:
p = photutils.aperture_photometry(mask, aper)
result.append(p.field('aperture_sum'))
else:
result.append(None)

return result

def write_coadd_images(band,
Expand Down
4 changes: 2 additions & 2 deletions py/legacypipe/format_catalog.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ def format_catalog(T, hdr, primhdr, allbands, outfn, release,
'fiberflux', 'fibertotflux']
if has_ap:
keys.extend(['apflux', 'apflux_resid', 'apflux_blobresid',
'apflux_ivar'])
'apflux_ivar', 'apflux_masked'])
_expand_flux_columns(T, bands, allbands, keys)

from tractor.sfd import SFDMap
Expand Down Expand Up @@ -181,7 +181,7 @@ def add_wiselike(c, bands=None):
add_fluxlike('fibertotflux')
if has_ap:
for c in ['apflux', 'apflux_resid', 'apflux_blobresid',
'apflux_ivar']:
'apflux_ivar', 'apflux_masked']:
add_fluxlike(c)

cols.extend(trans_cols_opt)
Expand Down
28 changes: 16 additions & 12 deletions py/legacypipe/runbrick.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1690,25 +1690,29 @@ def stage_coadds(survey=None, bands=None, version_header=None, targetwcs=None,
assert(C.AP is not None)
# How many apertures?
A = len(apertures_arcsec)
T.apflux = np.zeros((len(T), len(bands), A), np.float32)
T.apflux_ivar = np.zeros((len(T), len(bands), A), np.float32)
T.apflux_resid = np.zeros((len(T), len(bands), A), np.float32)
T.apflux = np.zeros((len(T), len(bands), A), np.float32)
T.apflux_ivar = np.zeros((len(T), len(bands), A), np.float32)
T.apflux_masked = np.zeros((len(T), len(bands), A), np.float32)
T.apflux_resid = np.zeros((len(T), len(bands), A), np.float32)
T.apflux_blobresid = np.zeros((len(T), len(bands), A), np.float32)
if Nno:
T_donotfit.apflux = np.zeros((Nno, len(bands), A), np.float32)
T_donotfit.apflux_ivar = np.zeros((Nno, len(bands), A), np.float32)
T_donotfit.apflux_resid = np.zeros((Nno, len(bands), A), np.float32)
T_donotfit.apflux = np.zeros((Nno, len(bands), A), np.float32)
T_donotfit.apflux_ivar = np.zeros((Nno, len(bands), A), np.float32)
T_donotfit.apflux_masked = np.zeros((Nno, len(bands), A), np.float32)
T_donotfit.apflux_resid = np.zeros((Nno, len(bands), A), np.float32)
T_donotfit.apflux_blobresid = np.zeros((Nno, len(bands), A), np.float32)
AP = C.AP
for iband,band in enumerate(bands):
T.apflux [:,iband,:] = AP.get('apflux_img_%s' % band)[:Nyes,:]
T.apflux_ivar [:,iband,:] = AP.get('apflux_img_ivar_%s' % band)[:Nyes,:]
T.apflux_resid[:,iband,:] = AP.get('apflux_resid_%s' % band)[:Nyes,:]
T.apflux [:,iband,:] = AP.get('apflux_img_%s' % band)[:Nyes,:]
T.apflux_ivar [:,iband,:] = AP.get('apflux_img_ivar_%s' % band)[:Nyes,:]
T.apflux_masked[:,iband,:] = AP.get('apflux_masked_%s' % band)[:Nyes,:]
T.apflux_resid [:,iband,:] = AP.get('apflux_resid_%s' % band)[:Nyes,:]
T.apflux_blobresid[:,iband,:] = AP.get('apflux_blobresid_%s' % band)[:Nyes,:]
if Nno:
T_donotfit.apflux [:,iband,:] = AP.get('apflux_img_%s' % band)[Nyes:,:]
T_donotfit.apflux_ivar [:,iband,:] = AP.get('apflux_img_ivar_%s' % band)[Nyes:,:]
T_donotfit.apflux_resid[:,iband,:] = AP.get('apflux_resid_%s' % band)[Nyes:,:]
T_donotfit.apflux [:,iband,:] = AP.get('apflux_img_%s' % band)[Nyes:,:]
T_donotfit.apflux_ivar [:,iband,:] = AP.get('apflux_img_ivar_%s' % band)[Nyes:,:]
T_donotfit.apflux_masked[:,iband,:] = AP.get('apflux_masked_%s' % band)[Nyes:,:]
T_donotfit.apflux_resid [:,iband,:] = AP.get('apflux_resid_%s' % band)[Nyes:,:]
T_donotfit.apflux_blobresid[:,iband,:] = AP.get('apflux_blobresid_%s' % band)[Nyes:,:]
del AP

Expand Down

0 comments on commit 0e2e7f5

Please sign in to comment.