Python iraf.set函数代码示例

    def __init__(self, input, **kwargs):
        """
        Init.

        :param input: FITS association
        :type input: string
        :param kwargs: additional keyword arguments
        :type kwargs: dictionary
        """
        self.input = input
        self.settings = dict(asndir='asn', rawdir='raw',
                             jref='/grp/hst/cdbs/jref/',
                             mtab='/grp/hst/cdbs/mtab/',
                             sourceImage='POL0V_drz.fits',
                             sigma=5.0)
        self.settings.update(kwargs)

        #add env variables to both system and IRAF
        os.putenv('jref', self.settings['jref'])
        os.putenv('mtab', self.settings['mtab'])
        iraf.set(jref=self.settings['jref'])
        iraf.set(mtab=self.settings['mtab'])

        #IRAF has funny True and False
        self.yes = iraf.yes
        self.no = iraf.no

def get_seeing(file, scale, ref):
    coords = file.replace('fits','dao')
    log = file.replace('fits','psf')
    out = file.replace('fits','seeing')

    iraf.noao()
    iraf.noao.obsutil()
    iraf.set(stdgraph="uepsf")

    data = iraf.psfmeasure(file, coords='markall', wcs='logical', display='no', frame=1, level=0.5, size='FWHM', radius=10.0, sbuffer=1.0, swidth=3.0, iterations=1, logfile=log, imagecur=coords, Stdout=1)

    fwhm = data.pop().split().pop()

    fwhm_pix = float(fwhm)
    fwhm = fwhm_pix*scale
    seeing = dimm_seeing(fwhm_pix, ref, scale)
    print "Seeing: ", seeing

    fp = open(out, 'w')
    fp.write("%f %f\n" % (fwhm_pix, seeing))
    fp.close()
    os.system("echo \"image;text 85 500 # text={Spot FWHM = %5.2f pixels}\" | xpaset WFS regions" % fwhm_pix)
    os.system("echo \'image;text 460 500 # text={Seeing = %4.2f\"}\' | xpaset WFS regions" % seeing)
    os.system("echo \"set wfs_seeing %4.2f\" | nc hacksaw 7666" % seeing)
    # New function call to insert seeing value into MySQL database.
    fits2mysql(file, mode, fwhm, seeing)

 def initialize_instrument(self, output_dir):
     self.logger.info("Initializing instrument.")
     inst_file_name = output_dir + '/tmp/cp.dat'
     if self.debug:
         iraf.set(debug=1)
     iraf.set(use_new_imt='no')
     iraf.noao.imred(_doprint=0)
     iraf.noao.imred.ccdred(_doprint=0)
     if not os.path.exists(inst_file_name):
         with open(inst_file_name, 'w+') as inst_file:
             inst_file.write('subset          FILTER\n\n')
             inst_file.write('darktime        EXPTIME\n\n')
             inst_file.write('\'Dark Frame\'    dark\n')
             inst_file.write('\'Bias Frame\'    zero\n')
             inst_file.write('\'Light Frame\'   object\n')
             inst_file.write('\'Flat Field\'    flat\n')
     iraf.noao.imred.ccdred.setinst(
         instrument = 'cp',
         review = 'no',
         mode = 'h',
         dir = output_dir + '/tmp/',
         site = ''
     )
     iraf.noao.digiphot(_doprint=0)
     iraf.noao.digiphot.apphot(_doprint=0)
     iraf.dataio(_doprint=0)
     iraf.noao.digiphot.ptools(_doprint=0)

    def display(self,extn=None):
        """Display SCI extensions from the Astrodata object
           and mark objects and reference stars found in the
           input file tables OBJCAT and REFCAT.

           DS9 must be running.
        """
        try:
            from stsci.numdisplay import display as ndisplay
        except ImportError:
            from numdisplay import display as ndisplay

        #file = self.file
        #ad = AstroData(file)

        ad = self.adout
        nscext = ad.count_exts('SCI')

        if extn != None:
            si,ei = extn,extn+1
        else:
            si,ei = 1,nscext+1
            
        fr = 1
        for xt in range(si,ei):
            #ff = '%s[%d]' % (file,xt)
            ff = ad['SCI',xt].data
            
            iraf.set (stdimage='imt47')
            ndisplay(ff, frame=fr)
    
            xtnad = ad['OBJCAT',xt]
            if xtnad:
                x = xtnad.data.field('x')
                y = xtnad.data.field('y')
                print 'Marking %d objects with red.' % len(x)
                np.savetxt('/tmp/xy.txt', zip(x,y), fmt='%.2f', delimiter=' ')
                tvmark(frame=fr,coords='/tmp/xy.txt',color=204)
            else:
                print "No stars found in extension %d."%xt
        
            xtnad = ad['REFCAT',xt]
            if xtnad:
                x = xtnad.data.field('x')
                y = xtnad.data.field('y')
                print 'Marking %d refstars with green.' % len(x)
                np.savetxt('/tmp/xy.txt', zip(x,y), fmt='%.2f', delimiter=' ')
                tvmark(frame=fr,coords='/tmp/xy.txt',color=205)
            else:
                print "No reference stars found in extension %d."%xt

            fr += 1

    def build(self, fwhm=None, maxlin=40000, zeropoint=None):
        """

        @param fwhm:  initial guess of the FWHM, used for selecting sources that might be stars
        @param thresh: the detection threshold for stars (in sky sigma)
        @param maxlin: maximum linearity of the CCD, brighter sources are saturated.
        @param zeropoint: the cts/exptime mag zeropoint of the image
        @return:
        """
        if not zeropoint:
            zeropoint = storage.get_zeropoint(self.expnum,
                                              self.ccd,
                                              self.prefix,
                                              self.version)
        if not fwhm:
            fwhm = storage.get_fwhm(self.expnum,self.ccd,prefix=self.prefix, version=self.version)

        filename = storage.get_image(self.expnum,
                                     self.ccd,
                                     version=self.version,
                                     prefix=self.prefix)
        basename = os.path.splitext(filename)[0]

        iraf.set(writepars=0)

        iraf.kinky
        iraf.pipeline
        phot_filename = "{}.{}".format(basename, PHOT_EXT)
        psf_filename = "{}.{}".format(basename, PSF_EXT)
        zeropt_filename = "{}.{}".format(basename, ZERO_PT_USED)
        with open(zeropt_filename, 'w') as z:
            z.write("{}\n".format(zeropoint))

        apmax = AP_OUT_SCALE * fwhm
        swidth = fwhm

        # Turn off verbosity and parameter checking in DAOPHOT
        iraf.daophot.verify=iraf.no
        iraf.daophot.verbose=iraf.no
        iraf.daophot.update=iraf.no
        iraf.psf.showplots=iraf.no
        iraf.psf.interactive=iraf.no

        iraf.centerpars.calgori="centroid"
        iraf.centerpars.maxshift=3

        # set the basic 'datapars' for the images.
	    iraf.datapars.exposure=""
        # Note the itime is set to 1 as the SGwyn ZP includes the exposure time.
	    iraf.datapars.itime=1

def marksn2(img,fitstab,frame=1,fitstab2='',verbose=False):
   from pyraf import iraf
   from numpy import array   #,log10
   import lsc
   iraf.noao(_doprint=0)
   iraf.digiphot(_doprint=0)
   iraf.daophot(_doprint=0)
   iraf.images(_doprint=0)
   iraf.imcoords(_doprint=0)
   iraf.proto(_doprint=0)
   iraf.set(stdimage='imt1024')
   hdr=lsc.util.readhdr(fitstab)
   _filter=lsc.util.readkey3(hdr,'filter')
   column=lsc.lscabsphotdef.makecatalogue([fitstab])[_filter][fitstab]

   rasex=array(column['ra0'],float)
   decsex=array(column['dec0'],float)


   if fitstab2:
      hdr=lsc.util.readhdr(fitstab2)
      _filter=lsc.util.readkey3(hdr,'filter')
      _exptime=lsc.util.readkey3(hdr,'exptime')
      column=lsc.lscabsphotdef.makecatalogue([fitstab2])[_filter][fitstab2]
      rasex2=array(column['ra0'],float)
      decsex2=array(column['dec0'],float)

   iraf.set(stdimage='imt1024')
   iraf.display(img + '[0]',frame,fill=True,Stdout=1)
   vector=[]
   for i in range(0,len(rasex)):
      vector.append(str(rasex[i])+' '+str(decsex[i]))

   xy = iraf.wcsctran('STDIN',output="STDOUT",Stdin=vector,Stdout=1,image=img+'[0]',inwcs='world',units='degrees degrees',outwcs='logical',\
                         formats='%10.1f %10.1f',verbose='yes')[3:]
   iraf.tvmark(frame,'STDIN',Stdin=list(xy),mark="circle",number='yes',label='no',radii=10,nxoffse=5,nyoffse=5,color=207,txsize=2)

   if verbose:
 #     print 2.5*log10(_exptime)
      for i in range(0,len(column['ra0'])):
         print xy[i],column['ra0'][i],column['dec0'][i],column['magp3'][i],column['magp4'][i],column['smagf'][i],column['magp2'][i]

   if fitstab2:
      vector2=[]
      for i in range(0,len(rasex2)):
         vector2.append(str(rasex2[i])+' '+str(decsex2[i]))
      xy1 = iraf.wcsctran('STDIN',output="STDOUT",Stdin=vector2,Stdout=1,image=img+'[0]',inwcs='world',units='degrees degrees',outwcs='logical',\
                            formats='%10.1f %10.1f',verbose='yes')[3:]
      iraf.tvmark(frame,'STDIN',Stdin=list(xy1),mark="cross",number='yes',label='no',radii=10,nxoffse=5,nyoffse=5,color=205,txsize=2)

def buildtmc(tmcname):
    from pyraf import iraf
    from iraf import stsdas,hst_calib,synphot
    out=open('buildtmc.log','w')
    f=pyfits.open(tmcname)
    flist=f[1].data.field('filename')
    iraf.set(crrefer='./') #work locally

    for k in range(len(flist)):
        oldname=iraf.osfn(flist[k]).split('[')[0]
        newname=fincre(oldname)
        if os.path.exists(newname):
            flist[k]=fincre(flist[k])
        else:
            out.write("%s: no change necessary\n"%oldname)
    f.writeto(tmcname.replace(".fits","_new.fits"))
    out.close()

def atmofile(imgstd, imgout=''):
    # print "LOGX:: Entering `atmofile` method/function in %(__file__)s" %
    # globals()
    from pyraf import iraf
    import os
    import ntt

    iraf.noao(_doprint=0)
    iraf.onedspec(_doprint=0)
    iraf.set(direc=ntt.__path__[0] + '/')
    _cursor = 'direc$standard/ident/cursor_sky_0'
    if not imgout:
        imgout = 'atmo_' + imgstd
    os.system('rm -rf ' + imgout)
    iraf.noao.onedspec.bplot(imgstd, cursor=_cursor,
                             spec2=imgstd, new_ima=imgout, overwri='yes')
    return imgout

def doPhot():
	print 'Running doPhot'
	try:
		mapFile = open('map.dat')
	except:
		print 'map.dat not found. Exiting.'
		sys.exit()
	try:
		changeFile = open('change.dat')
	except:
		print 'change.dat not found. Exiting.'
		sys.exit()

	coordList = list()
	for line in mapFile:
		coords = map(float,line.split())
		coordList.append(coords)
	coordList = np.array(coordList)

	iraf.noao()
	iraf.digiphot()
	iraf.daophot()
	iraf.ptools()
	iraf.set(clobber='yes')
	
	photFile = open('phot.dat', 'w')
	for line in changeFile.readlines():
		elms = line.split()
		imageName = elms[0]
		changeCoords = np.array([float(elms[1]),float(elms[2])])
		newCoords = coordList + changeCoords
		print 'Image: '+imageName
		coordFile = makeCoordFile(newCoords)
		
		iraf.noao.digiphot.daophot.phot(image=imagesDir+imageName, coords=coordFile, output='.temp-phot', skyfile='', verify='no', fwhmpsf=fwhmpsf, sigma=sigma, readnoise=readnoise, epadu=epadu, exposure=exposureHeader, obstime=obstimeHeader, calgorithm=calgorithm, cbox=cbox, apertures=apertures, annulus=annulus, dannulus=dannulus)

		result = iraf.noao.digiphot.ptools.txdump(Stdout=1, textfiles='.temp-phot', fields='mag, merr, otime', expr='yes')
		writeString = result[0].split()[-1] +' '+ ' '.join([' '.join(x.split()[:2]) for x in result])
		photFile.write(writeString+"\n")
		
	photFile.close()
	raw_input('doPhot Done. Hit return key to continue.')

 def setiraf(self):
     """
     Set IRAF global parameters and load DAOPHOT package for aperture
     photometry.
     
     Parameters
     ----------
     
     Returns
     -------
     None
     """
     iraf.prcacheOff()
     iraf.set(writepars=0)
 
     # Load IRAF packages
     iraf.noao(_doprint = 0)
     iraf.noao.digiphot(_doprint = 0)
     iraf.noao.digiphot.daophot(_doprint = 0)
 
     return

def get_seeing(file, scale, ref):

#    log = file.replace('fits','psf')		# never used

# Gets called with full path to the fits file
#    print "daofind get_seeing: ", file

    iraf.noao()
    iraf.noao.obsutil()
    iraf.set(stdgraph="uepsf")

    # name of file with coordinates of star images.
    coords = file.replace('fits','dao')

    data = iraf.psfmeasure(file, coords='markall', wcs='logical', display='no', frame=1, level=0.5, size='FWHM', radius=10.0, sbuffer=1.0, swidth=3.0, iterations=1, logfile=log, imagecur=coords, Stdout=1)

    # print "daofind psfmeasure data: ", data
    # last = data.pop()
    # data.append(last)
    # print "daofind psfmeasure last line of data: ", last
    # psfmeasure returns lots of information (data for each spot), but we just
    # want the last line (the average fwhm), which looks like:
    # Average full width at half maximum (FWHM) of 6.0163

    fwhm = data.pop().split().pop()
    print "daofind psfmeasure fwhm: ", fwhm

    fwhm_pix = float(fwhm)
    fwhm = fwhm_pix*scale
    seeing = dimm_seeing(fwhm_pix, ref, scale)
    print "Seeing: ", seeing

    seeingfile = file.replace('fits','seeing')
    fp = open(seeingfile, 'w')
    fp.write("%f %f\n" % (fwhm_pix, seeing))
    fp.close()

    os.system("echo \"image;text 85 500 # text={Spot FWHM = %5.2f pixels}\" | xpaset WFS regions" % fwhm_pix)
    os.system("echo \'image;text 460 500 # text={Seeing = %4.2f\"}\' | xpaset WFS regions" % seeing)
    os.system("echo \"set wfs_seeing %4.2f\" | nc hacksaw 7666" % seeing)

def get_seeing(file, scale, ref):
    coords = file.replace("fits", "dao")
    log = file.replace("fits", "psf")
    out = file.replace("fits", "seeing")

    iraf.noao()
    iraf.noao.obsutil()
    iraf.set(stdgraph="uepsf")

    data = iraf.psfmeasure(
        file,
        coords="markall",
        wcs="logical",
        display="no",
        frame=1,
        level=0.5,
        size="FWHM",
        radius=10.0,
        sbuffer=1.0,
        swidth=3.0,
        iterations=1,
        logfile=log,
        imagecur=coords,
        Stdout=1,
    )

    fwhm = data.pop().split().pop()

    fwhm_pix = float(fwhm)
    fwhm = fwhm_pix * scale
    seeing = dimm_seeing(fwhm_pix, ref, scale)
    print seeing
    fp = open(out, "w")
    fp.write("%f %f\n" % (fwhm_pix, seeing))
    fp.close()
    os.system('echo "image;text 85 500 # text={Spot FWHM = %5.2f pixels}" | xpaset WFS regions' % fwhm_pix)
    os.system("echo 'image;text 460 500 # text={Seeing = %4.2f\"}' | xpaset WFS regions" % seeing)
    os.system('echo "set wfs_seeing %4.2f" | nc hacksaw 7666' % seeing)

def dark_sky_flat(filter):
    med_otalist = []
    print 'making dark sky flats for',filter
    for key in tqdm(odi.OTA_dictionary):
        image_list = odi.OTA_dictionary[key]+'.'+filter+'.lis'
        med_out = image_list.replace('.lis','.med.fits')
        med_otalist.append(med_out)
        iraf.unlearn(iraf.immatch.imcombine)
        iraf.immatch.imcombine.setParam('input','@'+str(image_list))
        iraf.immatch.imcombine.setParam('output',odi.skyflatpath+med_out)
        iraf.immatch.imcombine.setParam('combine','median')
        iraf.immatch.imcombine.setParam('masktype','goodvalue')
        iraf.immatch.imcombine.setParam('maskvalue',0)
        iraf.immatch.imcombine.setParam('scale','median')
        # iraf.immatch.imcombine.setParam('zero','none')
        iraf.immatch.imcombine.setParam('zero','median')
        iraf.immatch.imcombine(logfile='imcombine.log.txt', mode='h')
        if key == 1:
            data,header = odi.fits.getdata(odi.skyflatpath+med_out,header=True)
            mean, median, std = odi.sigma_clipped_stats(data, sigma=3.0)
            normalization_factor = median
    iraf.set(clobber = 'yes')
    print 'smoothing dark sky flats for',filter
    for i in tqdm(range(len(med_otalist))):
        iraf.unlearn(iraf.imutil.imarith,iraf.imfilter.median)
        iraf.imutil.imarith.setParam('operand1',odi.skyflatpath+med_otalist[i])
        iraf.imutil.imarith.setParam('op','/')
        iraf.imutil.imarith.setParam('operand2',normalization_factor)
        iraf.imutil.imarith.setParam('result',odi.skyflatpath+med_otalist[i])
        iraf.imutil.imarith(verbose='no', mode='h')
        iraf.imfilter.median.setParam('input',odi.skyflatpath+med_otalist[i])
        iraf.imfilter.median.setParam('output',odi.skyflatpath+med_otalist[i])
        iraf.imfilter.median.setParam('xwindow',3)
        iraf.imfilter.median.setParam('ywindow',3)
        iraf.imfilter.median(verbose='no', mode='h')
    iraf.set(clobber = 'no')

    return normalization_factor

def get_star_data(asteroid_id, mag, expnum, header):
    """
    From ossos psf fitted image, calculate mean of the flux of each row of the rotated PSF
    """

    # calculate mean psf
    uri = storage.get_uri(expnum.strip('p'), header[_CCD].split('d')[1])
    ossos_psf = '{}.psf.fits'.format(uri.strip('.fits'))
    local_psf = '{}{}.psf.fits'.format(expnum, header[_CCD].split('d')[1])
    local_file_path = '{}/{}'.format(_STAMPS_DIR, local_psf)
    storage.copy(ossos_psf, local_file_path)

    # pvwcs = wcs.WCS(header)
    # x, y = pvwcs.sky2xy(asteroid_id['ra'].values, asteroid_id['dec'].values)
    x = asteroid_id[_XMID_HEADER].values[0]
    y = asteroid_id[_YMID_HEADER].values[0]

    # run seepsf on the mean psf image
    iraf.set(uparm="./")
    iraf.digiphot(_doprint=0)
    iraf.apphot(_doprint=0)
    iraf.daophot(_doprint=0)
    iraf.seepsf(local_file_path, local_psf, xpsf=x, ypsf=y, magnitude=mag)

    with fits.open(local_psf) as hdulist:
        data = hdulist[0].data

    th = math.degrees(asteroid_id[_THETA_HEADER].values[0])
    data_rot = rotate(data, th)
    data_rot = np.ma.masked_where(data_rot == 0, data_rot)

    data_mean = np.ma.mean(data_rot, axis=1)

    os.unlink(local_psf)
    os.unlink(local_file_path)

    return data_mean[np.nonzero(np.ma.fix_invalid(data_mean, fill_value=0))[0]]

def mask_create(images, no_unlearn=False):
    output_files = []
    # If the output file list doesn't exsist, create it.
    if len(output_files) == 0:
        output_files = [os.path.join(os.path.dirname(i), "mask_" + os.path.basename(i)) for i in images]
    # Create a temporary file to hold the input image list
    tfile = tempfile.mkstemp()
    os.write(tfile[0], "\n".join(images) + "\n")
    os.close(tfile[0])

    # Create a temporary file to hold the input image list
    tofile = tempfile.mkstemp()
    os.write(tfile[0], "\n".join(output_files) + "\n")
    os.close(tfile[0])

    try:

        # Set up to call the task
        iraf.set(clobber="yes")
        iraf.nproto()

        if no_unlearn == False:
            print "Objmasks is unlearning!!"
            iraf.objmasks.unlearn()

        # set parameters
        iraf.objmasks.images = "@" + tfile[1]
        iraf.objmasks.objmasks = "@" + tofile[1]
        iraf.objmasks.omtype = "boolean"
        iraf.objmasks()

    except iraf.IrafError, e:
        print "objmasks failed"
        print "error #" + str(e.errno)
        print "Msg: " + e.errmsg
        print "Task: " + e.errtask
        raise e