import math import numpy import pyfits import os infile = 'AGESVC1_small_raw.fits' clip = 2.0 niter =5 order =2 outfile = str(infile.split('.fits')[0])+str('_poly_')+str(order)+str('.fits') os.system('cls') # Changes to the directory where the script is located abspath = os.path.abspath(__file__) dname = os.path.dirname(abspath) os.chdir(dname) FitsFile = pyfits.open(infile) image = FitsFile[0].data header =FitsFile[0].header sizez = image.shape[0] sizey = image.shape[1] sizex = image.shape[2] # Create an array to hold the channel numbers channel = numpy.array([i for i in range(0,sizez)]) for xp in range(0,sizex): print xp for yp in range(0,sizey): spectra = numpy.array(image[:,yp,xp]) nanlist=[] for i in range(0,len(spectra)): if numpy.isnan(spectra[i]): nanlist.append(i) tempspec = numpy.delete(spectra,nanlist) tempchan = numpy.delete(channel,nanlist) # Create duplicate arrays to fit sigma-clipped polynomials #tempspec = numpy.array(spectra) if len(tempspec)>order : for i in range(0,niter): length = len(tempchan) # Fit a polynomial temppoly = numpy.polyfit(tempchan, tempspec,order) tempmodl = numpy.array([0.0 for k in range(0,length)]) for f in range(0,order+1): fn = float(-f + order) tempmodl = tempmodl + temppoly[f]*pow(tempchan,fn) tempspec = tempspec - tempmodl rms = numpy.std(tempspec) av = numpy.median(tempspec) # Find the location of values to remove deviants = numpy.where(abs(tempspec - av) > clip*rms) if (len(deviants)>0 and len(deviants) < length): tempspec = numpy.delete(tempspec,deviants) tempchan = numpy.delete(tempchan,deviants) # Fit a polynomial to the original spectrum, masking all deviant values poly = numpy.polyfit(tempchan,spectra[tempchan],order) model = numpy.array( [0.0 for k in range(0,len(spectra))]) for f in range(0,order+1): fn = float(-f + order) model = model + poly[f]*pow(channel,fn) fitted = spectra - model image[:,yp,xp] = fitted FitsFile = pyfits.writeto(outfile, image, header=header)