#!/usr/bin/env python # Want to strip a CMF file down to only entries with defined M and dM import sys, os import re import pyfits from numpy import * infile = sys.argv[1] try: outfile= sys.argv[2] except: # Calculate outfile from infile. Insert 'good' before trailing .cmf. # If input file doesn't end in .cmf behavior is undefined (basefile,ext)=os.path.splitext(infile) outfile=basefile+'.good'+ext primaryheader = pyfits.getheader(infile,0) # Good detections should have none of the following flags set # AND 0111110100111000 == 0 PM_SOURCE_MODE_DEFAULT = 0x00000000 # < Initial value: resets all bits PM_SOURCE_MODE_PSFMODEL = 0x00000001 # < Source fitted with a psf model (linear or non-linear) PM_SOURCE_MODE_EXTMODEL = 0x00000002 # < Source fitted with an extended-source model PM_SOURCE_MODE_FITTED = 0x00000004 # < Source fitted with non-linear model (PSF or EXT; good or bad) PM_SOURCE_MODE_FAIL = 0x00000008 # < Fit (non-linear) failed (non-converge, off-edge, run to zero) PM_SOURCE_MODE_POOR = 0x00000010 # < Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?) PM_SOURCE_MODE_PAIR = 0x00000020 # < Source fitted with a double psf PM_SOURCE_MODE_PSFSTAR = 0x00000040 # < Source used to define PSF model PM_SOURCE_MODE_SATSTAR = 0x00000080 # < Source model peak is above saturation PM_SOURCE_MODE_BLEND = 0x00000100 # < Source is a blend with other sourcers PM_SOURCE_MODE_EXTERNAL = 0x00000200 # < Source based on supplied input position PM_SOURCE_MODE_BADPSF = 0x00000400 # < Failed to get good estimate of object's PSF PM_SOURCE_MODE_DEFECT = 0x00000800 # < Source is thought to be a defect PM_SOURCE_MODE_SATURATED = 0x00001000 # < Source is thought to be saturated pixels (bleed trail) PM_SOURCE_MODE_CR_LIMIT = 0x00002000 # < Source has crNsigma above limit PM_SOURCE_MODE_EXT_LIMIT = 0x00004000 # < Source has extNsigma above limit PM_SOURCE_MODE_MOMENTS_FAILURE = 0x00008000 # < could not measure the moments PM_SOURCE_MODE_SKY_FAILURE = 0x00010000 # < could not measure the local sky PM_SOURCE_MODE_SKYVAR_FAILURE = 0x00020000 # < could not measure the local sky variance PM_SOURCE_MODE_BELOW_MOMENTS_SN = 0x00040000 # < moments not measured due to low S/N PM_SOURCE_MODE_BIG_RADIUS = 0x00100000 # < poor moments for small radius, try large radius PM_SOURCE_MODE_AP_MAGS = 0x00200000 # < source has an aperture magnitude PM_SOURCE_MODE_BLEND_FIT = 0x00400000 # < source was fitted as a blend PM_SOURCE_MODE_EXTENDED_FIT = 0x00800000 # < full extended fit was used PM_SOURCE_MODE_EXTENDED_STATS = 0x01000000 # < extended aperture stats calculated PM_SOURCE_MODE_LINEAR_FIT = 0x02000000 # < source fitted with the linear fit PM_SOURCE_MODE_NONLINEAR_FIT = 0x04000000 # < source fitted with the non-linear fit PM_SOURCE_MODE_RADIAL_FLUX = 0x08000000 # < radial flux measurements calculated PM_SOURCE_MODE_SIZE_SKIPPED = 0x10000000 # < size could not be determined badbitmask = PM_SOURCE_MODE_FAIL + PM_SOURCE_MODE_POOR + PM_SOURCE_MODE_SATSTAR + PM_SOURCE_MODE_BLEND + PM_SOURCE_MODE_EXTERNAL + PM_SOURCE_MODE_BADPSF + PM_SOURCE_MODE_DEFECT + PM_SOURCE_MODE_SATURATED + PM_SOURCE_MODE_CR_LIMIT + PM_SOURCE_MODE_EXT_LIMIT + PM_SOURCE_MODE_MOMENTS_FAILURE + PM_SOURCE_MODE_SKY_FAILURE + PM_SOURCE_MODE_SKYVAR_FAILURE + PM_SOURCE_MODE_BELOW_MOMENTS_SN + PM_SOURCE_MODE_SIZE_SKIPPED goodbitmask = PM_SOURCE_MODE_PSFMODEL + PM_SOURCE_MODE_FITTED # Open up the file into an array of HDUs hdulist = pyfits.open(infile) primaryheader = hdulist[0].header # Assume all .cmf files will have extensions and thus a Primary that consists of just the header hdu0 = pyfits.PrimaryHDU(data=None,header=primaryheader) newhdulist = pyfits.HDUList(hdu0) for hdu in hdulist[1:]: tbldata = hdu.data if tbldata is None: newhdulist.append(hdu) continue # Otherwise we've got some real data and will filter it x = tbldata.field('X_PSF') y = tbldata.field('Y_PSF') flags = tbldata.field('FLAGS') mag = tbldata.field('PSF_INST_MAG') magerr = tbldata.field('PSF_INST_MAG_SIG') qf = tbldata.field('PSF_QF') qf_thresh = 0.8 eps = 1e-7 # See if we are analyzing a difference .cmf; If so we won't require # goodbitmask because those aren't really set for difference images # So, sure, anyone mixing difference and normal image .cmfs in one file # is just being annoying, but this information is not in the # Primary Header Unit so we have to check here diffcmf = hdu.header.has_key('PPSUB.KERNEL') # print "DIFFCMF: ", diffcmf # There are different cuts and selections for normal and difference images if diffcmf: w = asarray([ i for i in range(len(tbldata)) if ( isfinite(mag[i]) & isfinite(magerr[i]) & (int(flags[i]) & badbitmask == 0) & (qf[i] > qf_thresh) & (x[i] > eps) & (y[i] > eps)) ] ) else: w = asarray([ i for i in range(len(tbldata)) if ( isfinite(mag[i]) & isfinite(magerr[i]) & (int(flags[i]) & badbitmask == 0) & (qf[i] > qf_thresh) & (x[i] > eps) & (y[i] > eps) & (int(flags[i]) & goodbitmask != 0)) ] ) # if len(w) > 500: sys.exit() # No good subtraction should have more than a couple of dozen real detections in sparse fields. #(w,) = where(isfinite(tbldata.field('PSF_INST_MAG')) & isfinite(tbldata.field('PSF_INST_MAG_SIG'))) #(w,) = where(isfinite(tbldata.field('PSF_INST_MAG')) & isfinite(tbldata.field('PSF_INST_MAG_SIG'))) # (w,) = where((tbldata.field('PSF_INST_MAG') > -50) & (tbldata.field('PSF_INST_MAG_SIG') > 0)) # w = tbldata.field('PSF_INST_MAG') > -50 # print "w: ", w # if len(w) > 0: # print [int(f) & goodbitmask for f in flags[w]] if len(w) == 0: print "No good detections found in '%s'" % (infile) # print "Not creating file '%s'" % (outfile) # sys.exit() if len(w) > 0: newtbldata = tbldata[w] else: newtbldata = None #print w print "%d / %d" % (len(w), len(tbldata)) newhdu = pyfits.BinTableHDU(data=newtbldata,header=hdu.header) newhdulist.append(newhdu) # Close the _old_ file hdulist.close() print "Writing new file: %s" % outfile # Write the _new_ file newhdulist.writeto(outfile, clobber=True)