Python iraf.phot函数代码示例

def p60photsub(ofile, inlis, refimage, refstars, ot, filter):
   
    refstars.set_mag(filter)
    coofile=open('ot.coo','w')
    ot.wcs2pix(refimage)
    coofile.write('%10.3f%10.3f\n' % (ot[0].xval, ot[0].yval))
    coofile.close()

    outfile=open(ofile, 'a+')
  
    p60photutils.p60zeropt(inlis, refstars, filter)

    subimage=iraffiles('@%s' % inlis)
    for image in subimage:

        root,ext=image.split('.')
        iraf.phot('%s.sub' % root, coo='ot.coo', output='%s.grb' % root,
                  interac=no)

        stars=Starlist('%s.mag' % root)
        refstars.wcs2pix(image)
        zp,zpu=stars.zeropt(refstars,method='mean',rejout=0)

        otstars=Starlist('%s.grb' % root)    
        
        [mjd,utshut]=get_head(image,['OBSMJD','UTSHUT'])

        if (len(stars)>0):
            outfile.write('%s\t%s\t%.3f\t%.3f\t%.3f\t%.3f\n' % (mjd, utshut, otstars[0].mag+zp, sqrt(pow(zpu,2)+pow(otstars[0].magu,2)), otstars[0].magu, zpu))

    outfile.close()

 def daophot(self, framenum, coords, outfile, apertures, verbose = "no"):
     """
     Aperture photometry of stars in the coords file using IRAF PHOT routine.
     
     Parameters
     ----------
     framenum : int
         Frame number in the image cube to perform aperture photometry on.
         
     coords : string
         Text file with coordinate of the stars
         
     outfile : string
         Text file to which photometry results are written
     
     Returns
     -------
     None
     """
     # load iraf packages
     self.setiraf()
     
     iraf.delete(outfile)
     iraf.phot(image = self.sci_file + "[,," + str(framenum) + "]", coords = coords, output = outfile, fwhmpsf = self.fwhmpsf, sigma = self.sigma, readnoise = self.readnoise, epadu = self.epadu, exposure = self.exposure, calgorithm = self.calgorithm, cbox = self.cbox, maxshift = self.maxshift, salgorithm = self.salgorithm, annulus = self.annulus, dannulus = self.dannulus, apertures = apertures, zmag = self.zmag, interactive = "no", verify = "no", verbose = verbose)        
     
     return

    def _get_photometry(self):
        """ Get the photometry for the target.

        If the target is a standard star, aperture photometry will be performed. For the moment nothing is done with
        the others, but in due time (TODO) photometry.py will be included here. """

        basename = "standards"
        fd, coords_file = tempfile.mkstemp(prefix=basename, suffix=".coords")
        os.write(fd, "{0} {1} \n".format(self.RA, self.DEC))
        os.close(fd)

        if self.objtype == "standard":
            iraf.noao(_doprint=0)
            iraf.digiphot(_doprint=0)
            iraf.apphot(_doprint=0)
            seeing = self.header.hdr[self.header.seeingk]
            photfile_name = self.header.im_name + ".mag.1"
            utilities.if_exists_remove(photfile_name)
            kwargs =  dict(output=photfile_name, coords=coords_file,
                      wcsin='world', fwhm=seeing, gain=self.header.gaink, exposure=self.header.exptimek,
                      airmass=self.header.airmassk, annulus=6*seeing, dannulus=3*seeing,
                      apert=2*seeing, verbose="no", verify="no", interac="no")
            iraf.phot(self.header.im_name, **kwargs)
            [counts] = iraf.module.txdump(photfile_name, 'FLUX', 'yes', Stdout=subprocess.PIPE)
            utilities.if_exists_remove(coords_file)
            return float(counts)

def getSkyMeanSDinAnnulus(ann,delta=5):
	iraf.noao()
	iraf.digiphot()
	iraf.apphot()
	iraf.photpars.setParam('apertures','1')
	iraf.phot.setParam('interactive','no')
	iraf.phot.setParam('image',fitsDir+fitsFile)
	iraf.phot.setParam('coords',fitsDir+fitsFile+".coo")
	outmag=".maglim"
	try:	
		os.remove(fitsDir+fitsFile+outmag)
	except:
		print "File does not exist BEFORE running phot, so no need to delete."
	iraf.phot.setParam('output',fitsDir+fitsFile+outmag)
	iraf.phot.setParam('interac','no')	
	iraf.fitskypars.setParam('annulus',str(ann))	
	iraf.fitskypars.setParam('dannulus',str(delta))
	## NO SIGMA CLIPPING! JUST TO BE SAFE: (6/2013)
	iraf.fitskypars.setParam('sloclip',"0")
	iraf.fitskypars.setParam('shiclip',"0")

	iraf.phot(fitsDir+fitsFile,Stdin=cr)	
	aa, nn, xx, ss = ao.readPhotMagFile(fitsDir,fitsFile,outmag)
	try:	
		os.remove(fitsDir+fitsFile+outmag)
	except:
		print "File not found to delete AFTER running phot; that's odd."
	return xx

def daofind(filelist):
	iraf.datapars.setParam("exposure", "EXPTIME")
	iraf.datapars.setParam("airmass", "AIRMASS")
	iraf.datapars.setParam("filter", "INSFILTE")
	iraf.datapars.setParam("ccdread", "RDNOISE")
	iraf.datapars.setParam("gain", "GAIN")

	iraf.centerpars.setParam("calgorithm", "centroid")

	iraf.photpars.setParam("apertures", "17.56")
	
	iraf.fitskypars.setParam("annulus", "20.56")
	iraf.fitskypars.setParam("dannulus", "5")
	
	

	with open(filelist) as file1:
		for line in file1:
			filename = line.strip()
			print("Filename %s"  % filename)
			y,z = np.loadtxt(filename + "-daoedit", usecols=[3,4], unpack=True)
			sigma = np.mean(y) 
			fwhm = np.mean(z)
			iraf.datapars.setParam("fwhmpsf", "%s" % str(fwhm)  )
			iraf.datapars.setParam("sigma", "%s" % str(sigma)  )
			iraf.daofind.setParam("mode", "h")
			iraf.phot.setParam("mode", "h")
			iraf.datapars.setParam("mode", "h")
			iraf.centerpars.setParam("mode", "h")
			iraf.fitskypars.setParam("mode", "h")
			iraf.photpars.setParam("mode", "h")
	
			iraf.daofind(filename, filename + "-daofind", mode='h')
			
			iraf.phot(filename, coords=(filename + "-daofind"), output=(filename+".mag"), mode='h')

def psffit2(img, fwhm, psfstars, hdr, _datamax, psffun='gauss',  fixaperture=False):
    ''' 
    giving an image, a psffile,  calculate the magnitudes of strs in the file _psf.coo
    '''
    import lsc
    _ron = lsc.util.readkey3(hdr, 'ron')
    _gain = lsc.util.readkey3(hdr, 'gain')
    if not _ron:
        _ron = 1
        print 'warning ron not defined'
    if not _gain:
        _gain = 1
        print 'warning ron not defined'

    iraf.digiphot(_doprint=0)
    iraf.daophot(_doprint=0)
    zmag = 0.
    varord = 0  # -1 analitic 0 - numeric

    if fixaperture:
        print 'use fix aperture 5 8 10'
        hdr = lsc.util.readhdr(img+'.fits')
        _pixelscale = lsc.util.readkey3(hdr, 'PIXSCALE')
        a1, a2, a3, a4, = float(5. / _pixelscale), float(5. / _pixelscale), float(8. / _pixelscale), float(
                    10. / _pixelscale)
    else:
        a1, a2, a3, a4, = int(fwhm + 0.5), int(fwhm * 2 + 0.5), int(fwhm * 3 + 0.5), int(fwhm * 4 + 0.5)

    iraf.fitskypars.annulus = a4
    iraf.fitskypars.salgori = 'mean'  #mode,mean,gaussian
    iraf.photpars.apertures = '%d,%d,%d' % (a2, a3, a4)
    iraf.datapars.datamin = -100
    iraf.datapars.datamax = _datamax
    iraf.datapars.readnoise = _ron
    iraf.datapars.epadu = _gain
    iraf.datapars.exposure = 'EXPTIME'  
    iraf.datapars.airmass = ''
    iraf.datapars.filter = ''
    iraf.photpars.zmag = zmag
    iraf.centerpars.calgori = 'centroid'
    iraf.centerpars.cbox = a2
    iraf.daopars.recenter = 'yes'
    iraf.delete('_psf2.ma*', verify=False)

    iraf.phot(img+'[0]', '_psf2.coo', '_psf2.mag', interac=False, verify=False, verbose=False)

    iraf.daopars.psfrad = a4
    iraf.daopars.functio = psffun
    iraf.daopars.fitrad = a1
    iraf.daopars.fitsky = 'yes'
    iraf.daopars.sannulus = a4
    iraf.daopars.recenter = 'yes'
    iraf.daopars.varorder = varord
    iraf.delete("_als2,_psf.grp,_psf.nrj", verify=False)
    iraf.group(img + '[0]', '_psf2.mag', img + '.psf', '_psf.grp', verify=False, verbose=False)
    iraf.nstar(img + '[0]', '_psf.grp', img + '.psf', '_als2', '_psf.nrj', verify=False, verbose=False)
    photmag = iraf.txdump("_psf2.mag", 'xcenter,ycenter,id,mag,merr', expr='yes', Stdout=1)
    fitmag = iraf.txdump("_als2", 'xcenter,ycenter,id,mag,merr', expr='yes', Stdout=1)
    return photmag, fitmag

def get_fake_centroid(filename,x,y,instrument,filt):
    """
    Locate the center of a fake psf

    INPUTS: The fake-SN psf image in filename, the expected x,y position
    of the center of the psf, the instrument and filter being modeled.
    RETURNS: xcentroid, ycentroid, fluxcorr
    """
    from pyraf import iraf

    iraf.digiphot(_doprint=0)
    iraf.apphot(_doprint=0)
    iraf.unlearn(iraf.apphot.phot)
    iraf.unlearn(iraf.datapars)
    iraf.unlearn(iraf.centerpars)
    #Use the centroid algorithm right now as it seems more robust to geometric distortion.
    iraf.centerpars.calgorithm = 'centroid'
    iraf.centerpars.cbox = 5.0

    iraf.unlearn(iraf.fitskypars)
    iraf.unlearn(iraf.photpars)
    photparams = {
        'interac':False,
        'radplot':False,
        }
    iraf.datapars.readnoise = 0.0
    iraf.datapars.itime = 1.0
    iraf.datapars.epadu = 1.0

    # iraf.digiphot.apphot.fitskypars :
    iraf.unlearn(iraf.fitskypars)
    iraf.fitskypars.salgorithm = 'constant'
    iraf.fitskypars.skyvalue = 0.0

    # iraf.digiphot.apphot.photpars :
    iraf.unlearn(iraf.photpars)
    iraf.photpars.weighting = 'constant'
    iraf.photpars.apertures = 20   # TODO : set this more intelligently !
    iraf.photpars.zmag = 25
    iraf.photpars.mkapert = False

    #Write the coordinate file starting as position x and y
    coxyfile = 'centroid.xycoo'
    coxy = open(coxyfile, 'w')
    print >> coxy, "%10.2f  %10.2f" % (x,y)
    coxy.close()
    if os.path.exists('centroid.mag'): os.remove('centroid.mag')
    iraf.phot(image=filename, skyfile='', coords=coxyfile, output='centroid.mag',
              verify=False, verbose=True, Stdout=1, **photparams)
    f = open('centroid.mag', 'r')
    maglines = f.readlines()
    f.close()
    xcentroid = float(maglines[76].split()[0])
    ycentroid = float(maglines[76].split()[1])

    return xcentroid,ycentroid

def photometry(color, images):
	#grab the aligned images
	#images=sorted(glob.glob('reduced/'+color+'align/*fits'))

	for i in images:
		iraf.phot(i, coords='reduced/'+color+'align/coords.lis',
			fwhmpsf=3,datamin=-500.,datamax=10000.,readnoise=20.,epadu=7.2,
			  exposur="EXPTIME",airmass="SECZ",obstime="JD",filter="IRFLTID",annulus=15,
		  	  dannulu=7, apertur="9",verify='no',output=i+".MC")
	return

def psffit2(img, fwhm, psfstars, hdr, _datamax=45000, psffun='gauss', fixaperture=False):
    import agnkey

    iraf.digiphot(_doprint=0)
    iraf.daophot(_doprint=0)
    zmag = 0.
    varord = 0  # -1 analitic 0 - numeric
    if fixaperture:
        print 'use fix aperture 5 8 10'
        hdr = agnkey.util.readhdr(img+'.fits')
        _pixelscale = agnkey.util.readkey3(hdr, 'PIXSCALE')
        a1, a2, a3, a4, = float(5. / _pixelscale), float(5. / _pixelscale), float(8. / _pixelscale), float(
                    10. / _pixelscale)
    else:
        a1, a2, a3, a4, = int(fwhm + 0.5), int(fwhm * 2 + 0.5), int(fwhm * 3 + 0.5), int(fwhm * 4 + 0.5)

    _center='no'
    iraf.fitskypars.annulus = a4
    iraf.fitskypars.dannulus = a4
    iraf.noao.digiphot.daophot.daopars.sannulus = int(a4)
    iraf.noao.digiphot.daophot.daopars.wsannul = int(a4)
    iraf.fitskypars.salgori = 'mean'  #mode,mean,gaussian
    iraf.photpars.apertures = '%d,%d,%d' % (a2, a3, a4)
    #    iraf.photpars.apertures = '%d,%d,%d'%(a2,a3,a4)
    iraf.datapars.datamin = -100
    iraf.datapars.datamax = _datamax
    iraf.datapars.readnoise = agnkey.util.readkey3(hdr, 'ron')
    iraf.datapars.epadu = agnkey.util.readkey3(hdr, 'gain')
    iraf.datapars.exposure = 'exptime'  #agnkey.util.readkey3(hdr,'exptime')
    iraf.datapars.airmass = 'airmass'
    iraf.datapars.filter = 'filter2'
    iraf.centerpars.calgori = 'gauss'
    iraf.centerpars.cbox = 1
    iraf.daopars.recenter = _center
    iraf.photpars.zmag = zmag

    iraf.delete('_psf2.ma*', verify=False)

    iraf.phot(img+'[0]', '_psf2.coo', '_psf2.mag', interac=False, verify=False, verbose=False)

    iraf.daopars.psfrad = a4
    iraf.daopars.functio = psffun
    iraf.daopars.fitrad = a1
    iraf.daopars.fitsky = 'yes'
    iraf.daopars.sannulus = int(a4)
    iraf.daopars.wsannul = int(a4)
    iraf.daopars.recenter = _center
    iraf.daopars.varorder = varord
    iraf.delete("_als,_psf.grp,_psf.nrj", verify=False)
    iraf.group(img+'[0]', '_psf2.mag', img + '.psf', '_psf.grp', verify=False, verbose=False)
    iraf.nstar(img+'[0]', '_psf.grp', img + '.psf', '_als', '_psf.nrj', verify=False, verbose=False)
    photmag = iraf.txdump("_psf2.mag", 'xcenter,ycenter,id,mag,merr', expr='yes', Stdout=1)
    fitmag = iraf.txdump("_als", 'xcenter,ycenter,id,mag,merr', expr='yes', Stdout=1)
    return photmag, fitmag

def recentre(image,refcoordfile):
    """Returns improved shift by centroiding
on the reference star using phot. This can be VERY
sensitive to the parameters in centerpars."""
    xin,yin = load(refcoordfile,unpack=True)
    try:
        iraf.phot(image,refcoordfile,'temp.mag',inter="no",calgorithm='centroid',
              mode='h',verify='no',update='no',verbose='no')
        xout,yout=iraf.pdump('temp.mag','xcen,ycen','yes',Stdout=1)[0].split()
    except:
        print "Recentring failed on", image
        return 0.,0.
    xout,yout = float(xout),float(yout)
    return xout-xin,yout-yin

def apply_phot(args):
    # Temporary file
    with tempfile.NamedTemporaryFile(suffix=".mag.1") as fd:
        phot_modstars = fd.name

    with tempfile.NamedTemporaryFile(suffix=".mag.1") as fd2:
        phot_subtstars = fd2.name

    hdr = args.image.header  # for short
    seeing, sigma = utils.get_from_header(args.image.im_name, hdr.seeingk, hdr.sigmak)
    iraf.phot(
        args.image.im_name,
        output=phot_modstars,
        coords=args.model_stars,
        wcsin=args.coords,
        fwhmpsf=seeing,
        sigma=sigma,
        ccdread=hdr.ccdronk,
        gain=hdr.gaink,
        exposure=hdr.exptimek,
        airmass=hdr.airmassk,
        annulus=6 * seeing,
        dannulus=3 * seeing,
        apert=3 * seeing,
        verbose="no",
        verify="no",
        interac="no",
    )

    iraf.phot(
        args.image.im_name,
        output=phot_subtstars,
        coords=args.subt_stars,
        wcsin=args.coords,
        fwhmpsf=seeing,
        sigma=sigma,
        ccdread=hdr.ccdronk,
        gain=hdr.gaink,
        exposure=hdr.exptimek,
        airmass=hdr.airmassk,
        annulus=6 * seeing,
        dannulus=3 * seeing,
        apert=3 * seeing,
        verbose="no",
        verify="no",
        interac="no",
    )
    return phot_modstars, phot_subtstars

def p48coaddphot(inlis, refstars, ot, filter, outfile):

    outf = open(outfile, "w")
    
    for image in inlis:

        root = image.split(".")[0]
        
        # Get the seeing in the image
        iqpkg.iqobjs(image, SIGMA, get_head(image, "SATURATE"), skyval="0.0")
        seepix = get_head(image, "SEEPIX")

        # Calculate zeropoint
        refstars.wcs2pix(image)
        refstars.set_mag(filter)
        xyfile = open("%s.xy" % root, "w")
        for star in refstars:
            xyfile.write("%10.3f%10.3f\n" % (star.xval, star.yval))
        xyfile.close()
        iraf.phot(image, coords="%s.xy" % root, output="%s.mag" % root,
                  aperture=1.2*float(seepix), interac=no)
        stars = Starlist("%s.mag" % root)
        zp, zpu = stars.zeropt(refstars,method="mean",rejout=0)

        # Measure source
        ot.wcs2pix(image)
        coofile = open("%s.coo" % root, "w")
        coofile.write("%10.3f%10.3f\n" % (ot[0].xval, ot[0].yval))
        coofile.close()
        iraf.phot(image, coords="%s.coo" % root, output="%s.ot" % root,
                  aperture=1.2*float(seepix), calgorithm="none", interac=no)
        stars = Starlist("%s.ot" % root)

        # Just return 99 for non-detection
        if len(stars)==1:
            smag = stars[0].mag
            smagu = stars[0].magu
            mag = smag + zp
        else:
            mag = 99.0
            smagu = 99.0
        mjdobs = DateTimeFrom(get_head(image, OBSDATEKEY)).mjd
        exptime = get_head(image, EXPTIMEKEY)
        outf.write("%15.3f%10.1f%10.3f%10.3f%10.3f\n" % (mjdobs, exptime, mag, smagu, zpu))

    outf.close()

def simplepsfphot(image,coords,psf,refstar,centre=True,vary=False):
    """PSF photometry, with a given PSF file in psf used for every image"""
    iraf.dele('temp.mag*')
    iraf.dele('temp.als')
    iraf.dele('temp.sub.fits')
    if centre:
        xsh,ysh = recentre(image,refstar) 
        print "Fine Centring: ", xsh,ysh
    else: xsh,ysh = 0,0
    if vary:
        setaperture(image,refstar)
    shift_file_coords(coords,xsh,ysh,'tempcoords2',sort='als')
    iraf.phot(image,'tempcoords2','temp.mag2',inter="no",calgorithm='none',
                  mode='h',verify='no',update='no',verbose='no')
    iraf.allstar(image,'temp.mag2',psf,'temp.als','temp.mag.arj','temp.sub.fits',
                 mode='h',verify='no',update='no',verbose='no')
    out = iraf.pdump('temp.als','id,mag,merr,msky','yes',Stdout=1)
    return out   

def createMagFile(inputDir, inputFile, useCoordFile, outputFile, fullyAutomatic=False, aperturesString='5,10,15,20', bestSkyAperture='100'):	
	cr=['\n']
	iraf.cd(inputDir)
	iraf.noao()
	iraf.digiphot()
	iraf.apphot()
	iraf.phot.setParam('interactive','no')
	iraf.centerpars.setParam('cbox','5')
	iraf.photpars.setParam('apertures',aperturesString)
	iraf.fitskypars.setParam('annulus',bestSkyAperture)
	iraf.fitskypars.setParam('dannulus','10')
	iraf.phot.setParam('image',inputFile)
	iraf.phot.setParam('coords',useCoordFile)
	iraf.phot.setParam('output',outputFile)
	if fullyAutomatic: ## If you use this option, it will not prompt you for any of the aperture/sky settings
		iraf.phot(mode='h',Stdin=cr)
	else: # Otherwise, you can change them at will, and get to press enter a bunch
		iraf.phot(mode='h')

def apphot(image,coords,refstar=None,centre=False,vary=False):
    """Apperture photometry with centering based on a reference star.
NB: centre refers to shifting the coordinates by centroiding on the
reference star; recentering on the final phot depends on
centerpars.calgorithm ."""
    iraf.dele('temp.mag*')
    if centre:
        xsh,ysh = recentre(image,refstar)
        print "Fine centring: ", xsh,ysh
    else: #no recentreing by reference star (but could still have calgorithm!=none)
        xsh,ysh = 0,0
    if vary:
        setaperture(image,refstar)
    shift_file_coords(coords,xsh,ysh,'tempcoords')
    iraf.phot(image,'tempcoords','temp.mag2',inter="no",
                  mode='h',verify='no',update='no',verbose='no')
    out = iraf.pdump('temp.mag2','id,flux,msky,stdev','yes',Stdout=1)
    return out