#!/usr/bin/python2 # # This program does the least squares fitting # to zero point data for the determination of # the photometric flat field of a mosaic camera. # #import glob #import copy #from sys import argv from irafdao import * from zplib import * #from Numeric import * import MLab #import Gnuplot #from Tkinter import * # get widget classes #from tkMessageBox import askokcancel, askquestion # get canned std dialog varlistNames = ['xmax', 'xmin', 'ymax', 'ymin', 'basisType', 'basisOrder', 'maxmerr', 'xsigclip', \ 'xplotMin', 'yplotMin', 'xplotMax','yplotMax'] def printUsageMessage(): print"""\n\n USAGE: zpAnalyze.py where options can be any, all, or none of: -help # print general help message -basisType [Chebyshev] # type of basis for fit (only Chebyshev) -basisOrder [7] # order of basis -xmax [8184.0} # maximum pixel range in x -xmin [0.0] # minimum pixel range in x -ymax [8196.0] # maximum pixel range in y -ymin [0.0] # minimum pixel range in y -maxmerr [1.0] # maximum magnitud error or difference -xsigclip [3.0] # maximum allowed deviation from fit in sigmas -xplotMin [-1.0] # minimum value of x coord of star for appearing in plots -xplotMax [1.0] # maximum value of x coord of star for appearing in plots -yplotMin [-1.0] # minimum value of y coord of star for appearing in plots -yplotMax [1.0] # maximum value of y coord of star for appearing in plots The numbers in the square brackets are the default values. For further help call the program with the option -help.\n\n """ return def printHelpMessage(): print """\n\n GENERAL HELP: This routine is used to do a 2D fit to magnitude difference vs offsets in x and y data. A zero point variation as a function of position over a CCD Mosaic is determined. It has been written explicitly to be used with the Wide Field Imager at the ESO/MPG 2.2-m telescope at La Silla. For the program to work appropriately the user needs to provide data with offsets in both x-y directions. The program has been shown to work well if the offsets are of 400 pixels amplitude. It does not work properly if they are on the order of 100 pixels, probably due to small scale flat field errors, as in the latter case the ff should be accurate to within 2 millimags. DEPENDENCIES: irafdao.py # local python library zplib.py # the library with necessary routines Numeric # the well know python lybrary Gnuplot # python's interface to gnuplot ALGORITHM: The routine first read the options and changes the default values, if there is any option. It then proceeds to read the data files, one by one. It will keep all the different lists separately throughout the work, as it is necessary to keep track of the individual mean values of each list. This mean value will be subtracted from each list at each iteration step. Then, 1. it goes into the initial clipping loop in which the user can interact with the program: after reviewing plots of the data the user can modify the values of maxmerr, xsigmaclip, xplotMin, xplotMax, yplotMin, and yplotMax. Although the change of xsigmaclip does not have an effect on the graphs been displayed at this stage the user can nonetheless change its value, which will be used in step 2. 2. The program uses the data clipped according to the values above to calculate the design matrix, and its singular value decomposition, that is, it does the least squares fit of the data. Then, 3. After the fit is done the result is plotted as a surface and contour plots, together with the residuals of the fit as a function of x, the magnitude differences as a function of x, and the x-y positions of the stars used for the fit. 4. Then we go to the second data clipping stage, in which the user can set the values of the different parameters interactively. And the program returns to step 2 after asking whether the user is satified with the cliping parameters, and if the user wants one more iterattion. After the final iteration the program outputs the surface plot as a ps file, and the polynomial coeficients in the file zpChebyshev.coefs. The chebyshev polynomila can then be plotted with the command: plot2DPol.py zpChebyshev.coefs Photometric data can then be corrected with the command zpCorrectPhot -pol zpChebyshev.coefs -xcol n1 -ycol n2 This command reads its standard input, it skips lines which start with #, and append a column at the end of all the other non-empty lines with the correction calculated using the polynomila in zpChebyshev.coefs, and the x y values from columns n1 and n2 of the input line. It is important to: 1) use the lowest polynomial order which erases any traces of systematic pattern in the fit residual's plot; 2) insure that the whole mosaic area is sampled with an adequate number of stars; 3) iterate until no more stars are left outside the sample and the surface does not appear to change anymore. EXAMPLES: 1. zpAnalyze.py -basisOrder 9 Vdx.zp Vdy.zp V2dx.zp V2dy.zp This is used to fit a 9th order 2D Chebyshev polynomial to the data in the files Vdx.zp, Vdy.zp, V2dx.zp, and V2dy.zp. 2. zpAnalyze.py -basisOrder 7 V*.zp It does a fit to all the files that are globed under V*.zp. \n\n """ return if __name__ == "__main__": # here we initialize the graphic surfaces # # this plots the residuals as multiples of merr gres = Gnuplot.Gnuplot(debug=0) gres('set data style points') # this plots the positions of the stars used for # the solution. Help keep the chebs tamed gxy = Gnuplot.Gnuplot(debug=0) gxy('set data style points') # this plots the dm21 gdm = Gnuplot.Gnuplot(debug=0) gdm('set data style points') # this plots the zp surface, that is, the actual # variation of zp with position g2d = Gnuplot.Gnuplot(debug=0) g2d('set parametric') g2d('set data style lines') g2d('set hidden') g2d('set cntrparam levels 10') g2d('set contour both') # here we load the default values for # the variables # # input data normalization options xmax = 8182.0 ymax = 8196.0 xmin = 0.0 ymin = 0.0 # least squares option basisType = "Chebyshev" basisOrder = 9 maxmerr = 1.00 xsigclip = 3.0 # residual plotting options xplotMin = -1.0 xplotMax = +1.0 yplotMin = -1.0 yplotMax = +1.0 varlistNames = ['xmax', 'xmin', 'ymax', 'ymin', 'basisType', 'basisOrder', 'maxmerr', 'xsigclip', \ 'xplotMin', 'yplotMin', 'xplotMax','yplotMax'] # if there is a -help option we just print the general help message # notice that this is the only no-argument option allowed. # # if we are creating a fake dataset we do that and exit # for token in argv: if token == '-fakeData': fakeData() sys.exit(0) elif token == '-help': printHelpMessage() sys.exit(0) # parsing options lineargs = getopts(argv) # if there is an optfile we read the variable values from it # if lineargs.has_key('-optfile'): #print lineargs['-optfile'] optDict = getoptsFromFile(lineargs['-optfile']) for i in range(len(varlistNames)): if optDict.has_key(varlistNames[i]): exec(varlistNames[i] + " = " + optDict[varlistNames[i]]) #print varlist[i] #print optDict # now we fill the rest of the variables in the # command line # for i in range(len(varlistNames)): thekey = '-' + varlistNames[i] if lineargs.has_key(thekey): exec(varlistNames[i] + " = " + lineargs[thekey]) # now we strip all the commands and their arguments and # whatever is left will be considered to be the list # of input files # parsedFileList = getFileList(argv) if len(parsedFileList) == 0: printUsageMessage() sys.exit(0) # here we parse the input file list # fullFileList = [] for fileName in parsedFileList: fullFileList.append(glob.glob(fileName)) fullFileList = [s[0] for s in fullFileList] finalList = [] for fileName in fullFileList: if fileName[-2:] == 'zp': finalList = finalList + [fileName] nFiles = len(finalList) #for aFile in finalList: # print aFile # here we create the top level Gui # to be finished at La Silla # root = Tk() # here we load the data onto the arrays # keeping track of where the data in each # list starts and end. Later we will need # to know this when we set the average of # each list to zero. startList = zeros((nFiles,), Int32) endList = zeros((nFiles,), Int32) xc1 = zeros((0,),Float64) yc1 = zeros((0,),Float64) xc2 = zeros((0,),Float64) yc2 = zeros((0,),Float64) dm21 = zeros((0,),Float64) merr = zeros((0,),Float64) mavg = zeros((0,),Float64) id = zeros((0,),Int32) counter = 0 for fileName in finalList: (inxc1,inyc1,inxc2,inyc2,indm21,inmerr,inmavg,inid) =\ readPhotDiffFrom(fileName) xc1 = concatenate((xc1,inxc1)) yc1 = concatenate((yc1,inyc1)) xc2 = concatenate((xc2,inxc2)) yc2 = concatenate((yc2,inyc2)) dm21 = concatenate((dm21,indm21)) merr = concatenate((merr,inmerr)) mavg = concatenate((mavg,inmavg)) id = concatenate((id,inid)) # we keep track of were each data list was # stored. This will come in handy later when # we do the subtraction of the average value # from each list, otherwise a general brightness # change between frames will dominate the fit. # nItems = len(inxc1) if counter == 0: startList[counter] = 0 else: startList[counter] = endList[counter-1] + 1 endList[counter] = startList[counter] + nItems - 1 print "%6d lines read from %s: i1 = %5d i2 = %5d" % \ (len(inid),fileName,startList[counter], endList[counter],) counter = counter + 1 nInFiles = counter print "%6d lines from %3d files: i1 = %5d i2 = %5d\n" % \ (len(xc1),nInFiles,startList[0],endList[nInFiles-1]) # here we fix possible merr=0.000 problems for i in range(len(merr)): if merr[i] == 0.000: merr[i] = 0.001 # here we normalize the coordinates to be within # -1 and +1, the range for the approximating # Chebyshev polynomials. # xc1 = (2.0*(xc1 - xmin) / (xmax-xmin)) - 1.0 yc1 = (2.0*(yc1 - ymin) / (ymax-ymin)) - 1.0 xc2 = (2.0*(xc2 - xmin) / (xmax-xmin)) - 1.0 yc2 = (2.0*(yc2 - ymin) / (ymax-ymin)) - 1.0 # and here we go into the main fit-clip loop # # we first force the clipped dm21 to have zero mean # for counter in xrange(nInFiles): i1 = startList[counter] i2 = endList[counter] avgdm = sum(dm21[i1:i2+1])/len(dm21[i1:i2+1]) dm21[i1:i2+1] = dm21[i1:i2+1] - avgdm print " ------> = %e subtracted from list %5d." % (i1,i2,avgdm,counter,) nMaxIter = 30 nMaxBadPoints = 0 done = 0 nTimes = 0 doneClipping = 0 while not done: nTimes = nTimes + 1 print "Proceeding with iteration %2d\n\n" % (nTimes,) while not doneClipping: # note that we fall into this loop only once (auxxc1, auxyc1, auxxc2, auxyc2, auxdm21, auxmerr) = copy.deepcopy((xc1, yc1, xc2, yc2, dm21, merr)) (auxstartList,auxendList) = copy.deepcopy((startList,endList)) # we clip the bad data here according to the maxmerr criteria. # badData1 = where(greater(auxmerr,maxmerr),1,0) badData2 = where(greater(abs(auxdm21),maxmerr),1,0) badData = logical_or(badData1,badData2) nTot = len(badData) nBad = sum(badData) print " ------> for these cliping parameters %6d out of %6d values would be clipped." % (nBad,nTot,) inList = [ auxxc1, auxyc1, auxxc2, auxyc2, auxdm21, auxmerr ] inListLength = len(auxxc1) (auxxc1,auxyc1,auxxc2,auxyc2,auxdm21,auxmerr,auxstartList,auxendList) = \ siftBadData(inList,auxstartList,auxendList,badData) outListLength = len(auxxc1) print " ------> actual number of data that would be clipped is %6d" % (inListLength-outListLength,) plotResiduals(gdm,"Delta magnitudes",auxdm21,auxxc1,auxyc1,xplotMin,xplotMax,yplotMin,yplotMax, new='n',hold='y') plotResiduals(gxy,"Star positions",auxyc1,auxxc1,auxyc1,-1.0,+1.0,-1.0,+1.0, new='n',hold='y') ans = askquestion('Verify clipping','Refine cliping\nparameters further?') if ans == 'no': doneClipping = 1 if not doneClipping: fields = 'maxmerr', 'sigclip', 'xplotMin', 'xplotMax', 'yplotMin', 'yplotMax' variables = maxmerr, xsigclip, xplotMin, xplotMax, yplotMin, yplotMax formFrame = Toplevel() vars = makeform(formFrame,fields,variables) Quitter(vars,formFrame,'Continue').pack(side=TOP) formFrame.bind('', (lambda event, v=vars: fetch(v))) formFrame.mainloop() formFrame.destroy() del formFrame returnList = [] for xvar in vars: returnList.append(xvar.get()) (maxmerr, xsigclip,xplotMin,xplotMax,yplotMin,yplotMax) = map(string.atof,tuple(returnList)) # we clip the bad data here according to the maxmerr criteria. # badData1 = where(greater(auxmerr,maxmerr),1,0) badData2 = where(greater(abs(auxdm21),maxmerr),1,0) badData = logical_or(badData1,badData2) nTot = len(badData) nBad = sum(badData) print " ------> initial bad data clipping: %6d out of %6d values to be clipped." % (nBad,nTot,) inList = [ auxxc1, auxyc1, auxxc2, auxyc2, auxdm21, auxmerr ] inListLength = len(auxxc1) (auxxc1,auxyc1,auxxc2,auxyc2,auxdm21,auxmerr,auxstartList,auxendList) = \ siftBadData(inList,auxstartList,auxendList,badData) outListLength = len(auxxc1) print " ------> actual number of data clipped is %6d" % (inListLength-outListLength,) badData1 = where(greater(merr,maxmerr),1,0) badData2 = where(greater(abs(dm21),maxmerr),1,0) badData = logical_or(badData1,badData2) nTot = len(badData) nBad = sum(badData) print " ------> initial bad data clipping: %6d out of %6d values to be clipped." % (nBad,nTot,) inlistLength = len(xc1) inList = [ xc1,yc1,xc2,yc2,dm21,merr ] (xc1,yc1,xc2,yc2,dm21,merr,startList,endList) = siftBadData(inList,startList,endList,badData) outlistLength = len(xc1) print " ------> actual number of data clipped is %6d" % (inListLength-outListLength,) # we then force the clipped dm21 to have zero mean # for counter in xrange(nInFiles): i1 = startList[counter] i2 = endList[counter] avgdm = sum(dm21[i1:i2+1])/len(dm21[i1:i2+1]) dm21[i1:i2+1] = dm21[i1:i2+1] - avgdm print " ------> = %e subtracted from list %5d." % (i1,i2,avgdm,counter,) # and we divide dm by the individual error estimates # normDm = dm21/merr # now we calculate the design matrix # print "\n ------> calculating design matrix" A = getDesignMatrix(xc1,yc1,xc2,yc2,merr,basisType,basisOrder) print " ------> design matrix A(%5d,%5d) obtained OK.\n" % (shape(A)[0],shape(A)[1],) # now we carry out the svd # print " ------> calculating singular value decomposition of design matrix" (u,w,vt)=MLab.svd(A) v = transpose(vt) print " ------> SVD done OK." # we calculate the threshold for validity of w # and threshold it # #print w wthresh = len(dm21)*2.2e-16*w[0] print " ------> thresholdin singular values with wthreshold = %e\n" % wthresh w = where(greater(w,wthresh),w,0.0) invw = where(greater(w,0.0),1.0/w,0.0) # now we determine the coefficients from the above matrices # sol= matrixmultiply(matrixmultiply(matrixmultiply(v,MLab.diag(invw)),transpose(u)),normDm) solt = transpose( [ sol ]) # before proceeding we plot the 2D solution just found # plotSurface2D(g2d,sol,basisType,basisOrder) # and their errors here # #sig2 = v[0,:]*0.0 #for i in xrange(v.shape[1]): # sig2 = sig2 + invw[i]*invw[i]*v[:,i]*v[:,i] #print "\n errors of parameters = " #print sqrt(sig2) # we determine the vector of residuals here # r = matrixmultiply(A,solt) - transpose( [ normDm ]) r = reshape(r, (size(r),)) # and we plot them here, just a trace through the center # of the mosaic # plotResiduals(gdm,'Delta magnitude',dm21,xc1,yc1,xplotMin,xplotMax,yplotMin,yplotMax, new='n',hold='y') plotResiduals(gres,"Fit residuals",r,xc1,yc1,xplotMin,xplotMax,yplotMin,yplotMax, new='n',hold='y') plotResiduals(gxy,"Star positions",yc1,xc1,yc1,-1.0,+1.0,-1.0,+1.0, new='n',hold='y') # and chi-square here # degreesOfFreedom = len(r)-len(sol) chisq = dot(r,r)/ degreesOfFreedom print " ------> normalized chi-square = %e (%6d d.f.) " % (chisq,degreesOfFreedom,) # we go into a loop whose pourpose is to find the # proper value for xsigclip # doneClipping = 0 while not doneClipping: (auxxc1, auxyc1, auxr, auxmerr, auxdm21) = copy.deepcopy((xc1, yc1, r, merr, dm21)) (auxstartList,auxendList) = copy.deepcopy((startList,endList)) # we clip the bad data here according to the sigclip criteria. # badData1 = where(greater(abs(auxr),xsigclip),1,0) # we clip the bad data here according to the maxmerr criteria. # badData2 = where(greater(auxmerr,maxmerr),1,0) badData3 = where(greater(abs(auxdm21),maxmerr),1,0) badData4 = logical_or(badData2,badData3) badData = logical_or(badData1,badData4) nTot = len(badData) nBad = sum(badData) print " ------> for these cliping parameters %6d out of %6d values would be clipped." % (nBad,nTot,) inList = [ auxxc1, auxyc1, auxr, auxmerr, auxdm21 ] inListLength = len(auxxc1) (auxxc1,auxyc1,auxr,auxmerr,auxdm21,auxstartList,auxendList) = siftBadData(inList,auxstartList,auxendList,badData) outListLength = len(auxxc1) print " ------> actual number of data that would be clipped is %6d" % (inListLength-outListLength,) plotResiduals(gdm,'Delta magnitude',auxdm21,auxxc1,auxyc1,xplotMin,xplotMax,yplotMin,yplotMax, new='n',hold='y') plotResiduals(gres,'Fit residuals',auxr,auxxc1,auxyc1,xplotMin,xplotMax,yplotMin,yplotMax, new='n',hold='y') plotResiduals(gxy,"Star positions",yc1,xc1,yc1,-1.0,+1.0,-1.0,+1.0, new='n',hold='y') ans = askquestion('Verify clipping','Refine cliping\nparameters further?') if ans == 'no': doneClipping = 1 if not doneClipping: fields = 'maxmerr', 'sigclip', 'xplotMin', 'xplotMax', 'yplotMin', 'yplotMax' variables = maxmerr, xsigclip, xplotMin, xplotMax, yplotMin, yplotMax formFrame = Toplevel() vars = makeform(formFrame,fields,variables) Quitter(vars,formFrame,'Continue').pack(side=TOP) formFrame.bind('', (lambda event, v=vars: fetch(v))) formFrame.mainloop() formFrame.destroy() del formFrame returnList = [] for xvar in vars: returnList.append(xvar.get()) (maxmerr, xsigclip,xplotMin,xplotMax,yplotMin,yplotMax) = map(string.atof,tuple(returnList)) # with the value of xsigclip determined above # we determine the bad data here. # badData = where(greater(abs(r),xsigclip),1,0) # and we clipp it here # nTot = len(badData) nBad = sum(badData) print " ------> bad data clipping: %6d out of %6d values to be clipped." % (nBad,nTot,) inList = [ xc1, yc1, xc2, yc2, dm21, merr ] inListLength = len(xc1) (xc1,yc1,xc2,yc2,dm21,merr,startList,endList) = siftBadData(inList,startList,endList,badData) outListLength = len(xc1) print " ------> actual number of data clipped is %6d" % (inListLength-outListLength,) if nBad <= nMaxBadPoints: done = 1 else: done = 0 if nTimes > nMaxIter: done = 1 if done == 0: done = 1 s = raw_input("\n\n ------> enter 'y' to proceed with iteration " + \ str(nTimes+1) + " > ") if s == 'y' or s == 'Y': print "\n\n\n" done = 0 # now we finish with the iterative search for a solution # we save the solution to a file # fileName = 'zpChebyshev.coefs' writePol2D(sol,basisType,basisOrder,fileName) # let us just check (btype,border,savedCoefs) = readPol2D(fileName) print btype,border # eso es to.., eso es to.., eso es todo amigos!