Python iraf.imarith函数代码示例

def flatCorrection():
	print "FlatCorrection start"
	#put all others params to no , they may be set by previous actions
	#WITH ccdproc
#	iraf.ccdproc.setParam('flatcor', 'yes')
#	iraf.ccdproc.setParam('fixpix', 'no')
#	iraf.ccdproc.setParam('darkcor', 'no')
#	iraf.ccdproc.setParam('illumcor', 'no')
#	iraf.ccdproc.setParam('trim', 'no')
#	iraf.ccdproc.setParam('overscan', 'no')
#	iraf.ccdproc.setParam('zerocor', 'no')
#	iraf.ccdproc.setParam('trimsec', '')
#	iraf.ccdproc.setParam('biassec', '')
#	#online
#	iraf.ccdproc.setParam('output', '')
	for f in FILTERS:	
		flatfilename = os.path.join(OUTPUTDIR, "flat", f, "FlatNorm.fits")
		if os.path.isfile(flatfilename):
#			#WITH ccdproc
#			iraf.ccdproc.setParam('flat', flatfilename)
#			iraf.ccdproc.setParam("images", os.path.join(OUTPUTDIR, "object", f) + "/*.fits")
#			iraf.ccdproc()
			with open(os.path.join(OUTPUTDIR, "object", f, "list")) as file1:
				for line in file1:
					objfilename = line.strip()
					if objfilename!="":
						objFullfilename = os.path.join(OUTPUTDIR, "object", f,	line.strip())
						iraf.imarith(objFullfilename, '/', flatfilename, objFullfilename)
		else:
			print "Flat file %s not present" % flatfilename
	print "FlatCorrection end"

 def sub_darks(self,data):
     for exp in data.keys():
         inputlist = data[exp]['raw']
         outputlist = [s[s.rfind('/')+1:s.find('.FIT')]+\
                           '_ds.fits'\
                           for s in inputlist]
         
         if self.darks[exp]['combined'] == None: self.gen_dark(exp)
         if len(inputlist) > 10:
             while len(inputlist) > 0:
                 inputstring = ','.join(inputlist[:10])
                 outputstring = ','.join(outputlist[:10])
                 iraf.imarith(inputstring,'-',\
                                  self.darks[exp]['combined'],\
                                  outputstring)
                 data[exp]['ds'] += outputlist[:10]
                 inputlist = inputlist[10:]
                 outputlist = outputlist[10:]
             
         else:
             iraf.imarith(','.join(inputlist),\
                              '-',\
                              self.darks[exp]['combined'],\
                              ','.join([s[s.rfind('/')+1:s.find('.FIT')]+\
                                            '_ds.fits' for s in inputlist]))
             data[exp]['ds'] += outputlist

def skySub(dfView):
    """
	given the dfview, and using the skyflat calculated in makeSkyFlat,
	subtract the skyflat from each image
	the output will be one sky-subtracted image for each image in the dfView
	the names of the sky-subtracted fits images have 'sky' infront of them
	"""
    # grab the paths to images to be sky-subtracted
    images = dfView.file.values.tolist()
    # adding a 's' infront of each file name, save to a seperate list
    skySubImages = ["s-" + i[5:] for i in images]
    # grab the YYMMDD observation date from this view
    date = str(dfView.Date.values[0])

    # organize the input, skyflat, and output in iraf-friendly ways
    inputFiles = joinStrList(images)
    skyFlat = "scratch/" + date + "sky.fits"
    outputSkySub = joinStrList(skySubImages, scratch=True)
    iraf.imarith(
        inputFiles,
        "-",
        skyFlat,
        outputSkySub,
        divzero=0.0,
        hparams="",
        pixtype="",
        calctype="",
        verbose="yes",
        noact="no",
    )
    return

def createFlatFiles():
	import re
	print "CreateFlatFiles start"
	for f in FILTERS:
		if(os.listdir(os.path.join(OUTPUTDIR, "flat", f))):
			iraf.imcombine.setParam("input", os.path.join(OUTPUTDIR, "flat", f) + "/*.fits")
			flatFile = os.path.join(OUTPUTDIR, "flat", f , "Flat.fits")
			if os.path.exists(flatFile):
				print("flatFile %s alreday exists deleting"  % flatFile)
				os.remove(flatFile)
			iraf.imcombine.setParam("output", flatFile)
			#from doc:	
			#http://www.iac.es/sieinvens/siepedia/pmwiki.php?n=HOWTOs.PythonianIRAF
			#--> iraf.listpix(mode='ql')     # confirms parameter
			#--> iraf.listpix(mode='h')     # doesn't ask for [email protected]@
			iraf.imcombine(mode="h")
			#NORMALIZE
			#imstat
			res = iraf.imstat(flatFile, Stdout=1)
			print(res[0].strip()) 
			print(res[1].strip()) 
			resArray = re.split("\s+", res[1].strip())
			#max value
			#toDivValue = float(resArray[5])
			#meanValue
			toDivValue = float(resArray[2])
			flatNormFile = os.path.join(OUTPUTDIR, "flat", f , "FlatNorm.fits")
			if os.path.exists(flatNormFile):
				print("flatNormFile %s alreday exists deleting"  % flatNormFile)
				os.remove(flatNormFile)
			#divide by max value
			iraf.imarith(flatFile, '/', toDivValue, flatNormFile)
		else:
			print("NO FLAT FILES for filter %s PRESENT" %f)
	print "CreateFlatFiles end"

def ExpNormalize(images, outbase="_n"):
			
	# Build the list of output image names
	out = [os.path.splitext(i) for i in images]	# Split off the extension
	out = [i[0] + outbase + '.fits' for i in out]	# Paste the outbase at the end of the filename 
												# and put the extension back on
	# Get a list of exposure times.
	exp_times = [GetHeaderKeyword(i, 'exptime') for i in images]
	exp_times = [str(e) for e in exp_times]	
	
	# run imarith to do the normalization
	iraf.imarith.unlearn()
	iraf.imarith.op = '/'
	iraf.imarith.mode = 'h'

	for i in range(len(images)):
		iraf.imarith.operand1 = images[i]
		iraf.imarith.operand2 = exp_times[i]
		iraf.imarith.result = out[i]
	
		iraf.imarith()

		# update the exptime keyword		
		iraf.hedit.unlearn()
		iraf.hedit.verify='no'
		iraf.hedit.show='yes'
		iraf.hedit.update='yes'
		iraf.hedit.images=out[i]
		iraf.hedit.fields='exptime'
		iraf.hedit.value=1
		iraf.hedit.mode='h'

		iraf.hedit(Stdout=1)

	return out

 def unit_converter(self):
     '''
     Converts DRZ images from units of counts/sec to counts
     '''
     
     for image, exptime in zip(self.imlist, self.explist):
         iraf.imarith(operand1=image + '[1]', op='*', operand2=exptime, \
                      result=image[:9] + '_counts.fits')

    def test_hedit(self):
        if os.path.exists('image.real.fits'):
            os.remove('image.real.fits')

        iraf.imarith('dev$pix', '/', '1', 'image.real', pixtype='r')
        iraf.hedit('image.real', 'title', 'm51 real', verify=False,
                   Stdout="/dev/null")
        with fits.open('image.real.fits') as f:
            assert f[0].header['OBJECT'] == 'm51 real'

def custom1(filename): # for NACO timing mode cubes - removes horizontal banding
    #iraf.imarith(filename,'-','dark','temp')
    iraf.imarith(filename,'/','flatK','temp')
    im = pyfits.getdata('temp.fits')
    med = median(im.transpose())
    out = ((im).transpose()-med).transpose()
    (pyfits.ImageHDU(out)).writeto("temp2.fits",clobber=True)
    iraf.imdel('temp')
    iraf.imcopy('temp2[1]','temp')

 def normalise(self, out=None):
     """
     Normalise the image.
     """
     if out is None:
         name, ext = os.path.splitext(self.filename)
         out = name + "_norm" + ext
     iraf.imarith(self.filename, "/", self.MEAN, out)
     return out

 def test_imarith(self):
     iraf.imarith('dev$pix', '/', '1', 'image.real', pixtype='r')
     with fits.open('image.real.fits') as f:
         assert f[0].header['BITPIX'] == -32
         assert f[0].data.shape == (512, 512)
     iraf.imarith('dev$pix', '/', '1', 'image.dbl', pixtype='d')
     with fits.open('image.dbl.fits') as f:
         assert f[0].header['BITPIX'] == -64
         assert f[0].data.shape == (512, 512)

def run_scldark(input):
    iraf.image(_doprint=0)
    iraf.image.imutil(_doprint=0)
    scl=float(input[1])/float(input[2])
#    print input[1],scl
    iraf.imarith(input[0],"*",scl,input[3],divzero="0.0",hparams="",pixtype="",calctype="",verbose="no",noact="no")
    iraf.hedit(input[3],"EXPTIME",input[1], add="yes", addonly="yes", delete="no", verify="no", show="no", update="yes")
    iraf.hedit(input[3],"DARKTIME",input[1], add="yes", addonly="yes", delete="no", verify="no", show="no", update="yes")
    iraf.hedit(input[3],"EXPOSURE","", add="no", addonly="no", delete="yes", verify="no", show="no", update="yes")

def MatchNSubtract(TargetImg,Template,OutputImage,fitgeometry="general"):
    """ Creates OutputImage =  TargetImg - Template after scaling and matching Template to TargetImg.
    fitgeometry can be set to 'rotate' when the Template is also TIRSPEC data
    Otherwise if Template is 2MASS or other instrument set it as 'general'  """
    
    AlignedImg = os.path.splitext(TargetImg)[0]+"_"+os.path.basename(Template)
    AlignedImg = os.path.splitext(AlignedImg)[0][:115]+'.fits' # Reduce filename length for iraf geopmap
    TransformDBfile = AlignImage(TargetImg,Template,AlignedImg,fitgeometry=fitgeometry)
    
    # Now get the Good sky region coordinates
    SkyCoordsFile = os.path.splitext(TargetImg)[0]+'_BlankSky.coo'
    if not os.path.isfile(SkyCoordsFile) :
        iraf.display(TargetImg,1)
        print ('For taking coordinates of good sky. Press _x_ over blank sky areas.')
        imx=iraf.imexam(Stdout=1)
        with open(SkyCoordsFile,'w') as foo :    #Creating blank sky coords files
            for line in imx :               
                foo.write(line.split()[0] +'  '+line.split()[1]+'\n')

    # Now get the regions in the image whose brightness has to be cancelled by scaling
    FluxCoordsFile = os.path.splitext(TargetImg)[0]+'_FluxRegions.coo'
    if not os.path.isfile(FluxCoordsFile) :
        iraf.display(TargetImg,1)
        print ('Press _x_ over areas you want to minimise the residual flux after subtraction')
        imx=iraf.imexam(Stdout=1)
        with open(FluxCoordsFile,'w') as foo :    #Creating Flux ares which we have to remove in subtraction
            for line in imx :               
                foo.write(line.split()[0] +'  '+line.split()[1]+'\n')

    #Now we first has to remove background from both the images.
    TargetSkySubtractedFile = os.path.splitext(TargetImg)[0]+'_SkyS.fits'
    if not os.path.isfile(TargetSkySubtractedFile):
        skyvalue = SkySubtractImage(TargetImg,TargetSkySubtractedFile,SkyCoordsFile)
    else:
        print('Warning: Using old {0} file'.format(TargetSkySubtractedFile))

    AlignedSkySubtractedFile = os.path.splitext(AlignedImg)[0]+'_SkyS.fits'
    if not os.path.isfile(AlignedSkySubtractedFile):
        skyvalue = SkySubtractImage(AlignedImg,AlignedSkySubtractedFile,SkyCoordsFile)
    else:
        print('Warning: Using old {0} file'.format(AlignedSkySubtractedFile))

    #We shall now extract the totel Flux in each tiles from both the images
    TargetFluxinTiles = ExtractTiles(TargetSkySubtractedFile,FluxCoordsFile,Summeryfunction=np.sum,hsize=7*1.5)
    TemplateFluxinTiles = ExtractTiles(AlignedSkySubtractedFile,FluxCoordsFile,Summeryfunction=np.sum,hsize=7*1.5)
    
    def DiffSquareSum(x):
        return np.sum([(targetF - x*templateF)**2 for targetF,templateF in zip(TargetFluxinTiles,TemplateFluxinTiles)])
    
    res = scipy.optimize.minimize_scalar(DiffSquareSum)
    Scale = res.x
    print('Scaling to match the fluxes is {0}'.format(Scale))
    iraf.imarith(operand1=AlignedSkySubtractedFile,op="*",operand2=Scale,result=os.path.splitext(AlignedSkySubtractedFile)[0]+'M.fits')

    iraf.imarith(operand1=TargetSkySubtractedFile,op="-",operand2=os.path.splitext(AlignedSkySubtractedFile)[0]+'M.fits',result=OutputImage)

def imarith(operand1, op, operand2, result):
    from pyraf import iraf
    iraf.images()

    pars = iraf.imarith.getParList()
    iraf.imcombine.unlearn()

    print "%s %s %s -> %s" % (operand1, op, operand2, result)
    iraf.imarith(operand1=operand1, op=op, operand2=operand2, result=result)

    iraf.imarith.setParList(pars)

	def get_median(self):
		'''
		Calculates the median subtracted image used to compute the pixel mask
		'''
		iraf.median(input = self.image, output = 'tmp_med.fits', 
					xwindow = 40, ywindow = 40, verbose ='No')

		iraf.imarith(operand1 = self.image, operand2 = 'tmp_med.fits[0]', 
					 op = '-', result = self.output + '_sub.fits')
		
		os.remove('tmp_med.fits')

def mask_norm(mask):
    # normalizes mask.
    # test if trimmed mask exists:
    if os.path.isfile(mask):
        iraf.cd(os.path.dirname(mask))
	hdulist = pyfits.open(mask)
	imdata = hdulist[0].data
	hdulist.close()
	outname = mask.replace('_trimmed.fits', '_normed.fits')
	print imdata.max()
	iraf.imarith(os.path.basename(mask),"/",imdata.max(),os.path.basename(outname))