Index: /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py =================================================================== --- /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 20526) +++ /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 20527) @@ -62,4 +62,9 @@ # XXX Todo: labels! group by runs? different colours for different # bands; include seeing? plot things against seeing? + +# Also: the plotcol lists need to be grouped into things that are to +# appear on the same page. For the beginning, that's the entire list +# (and that's probably quite enough anyway). + plotcol_1frame_tlist = [ @@ -107,6 +112,7 @@ def compIPPphoto(summaryTable,mode,plotcol_1frame_tlist=plotcol_1frame_tlist,\ plotcol_summary_tlist=plotcol_summary_tlist,\ - cmfDir = '/IPP/data/SDSS/stripe82/coadd/prod/run_ipp_20080815/',\ - skip=True): + cmfDir = '/IPP/data/SDSS/stripe82/coadd/prod/run_ipp_20080815/',\ + copyfields_list = ['RUN','RERUN','CAMCOL','FIELD','FILTER','FWHM_MAJ','FWHM_MIN'],\ + skip=True): """summaryTable: .fits table for output mode: new or append for creating summaryTable new or appending current run's output to it. @@ -144,5 +150,5 @@ filters.append(filter_name) bandindex_hash[filter_name] = bandindex - stats_hash, column_hash, goodrow_hash = computeStatistics(matchtable) + stats_hash, column_hash, goodrow_hash = computeStatistics(matchtable,copyfields_list = copyfields_list) # Sort res_hash by its column names to make sure order is # identical in all rows of the table, and output is more @@ -284,5 +290,33 @@ os.remove(tmp_outname) return all_outname - + +def setWindow(i,Nplots): + """Set sm window for plotting the i-th out of Nplots. The logic is + that N_x * N_y >= Nplots and N_y-N_x is 0 or 1. So Ny is the sqrt + of N or the next-greatest square; Nx is either the same or one + fewer if that's still enough.""" + + from scipy import sqrt + if int(sqrt(Nplots))**2 == Nplots: + Ny = int(sqrt(Nplots)) + Nx = Ny + else: + Ny = int(sqrt(Nplots)+1) + Nx = Ny-1 + if Nx * Ny < Nplots: + Nx+=1 + + # That's nx, ny; now count up to total Nx ny. This is like + # representing i in a number system where the y-window number is + # the 10s and the x-window number the 1s, but the 10s are base Ny + # and the 1s are base Nx, and counting starts at 1 not 0. I.e. x = + # (i % Nx) or Nx if that's 0; y = i / Nx + 1 + winx = i% Nx + if winx == 0: + winx = Nx + winy = int(i/Nx)+1 + if int(i/Nx) == i/float(Nx): + winy -= 1 + return Nx,Ny,winx,winy def valuesKeysSortedByKeys(hash): @@ -379,4 +413,9 @@ # "reference columns" outhash = {} + if 'FWHM_MAJ' in copyfields_list and 'FWHM_MAJ' not in h.ascardlist().keys(): + copyfields_list.remove('FWHM_MAJ') + copyfields_list.remove('FWHM_MIN') + copyfields_list.append('FWHM_X') + copyfields_list.append('FWHM_Y') for f in copyfields_list: outhash[f] = h[f] @@ -557,5 +596,8 @@ filters = ['u','g','r','i','z'] # Read primary header of PS1 file to work out band - ps1copyfields_hash = headerfieldHash(['FWHM_MAJ','FWHM_MIN','FILTER'],PS1cmf,HDU=0) + try: + ps1copyfields_hash = headerfieldHash(['FWHM_MAJ','FWHM_MIN','FILTER'],PS1cmf,HDU=0) + except KeyError: + ps1copyfields_hash = headerfieldHash(['FWHM_X','FWHM_Y','FILTER'],PS1cmf,HDU=0) sdssbandstr = ps1copyfields_hash['FILTER'] bandindex = filters.index(sdssbandstr) @@ -832,5 +874,4 @@ """Fit Gaussians to image in infile's HDU near the positions and with fluxes and sigmas given in initialGuessList""" - # XXX Todo: Make frame with residuals. Make blank array, then insert residuals. import pyfits from numpy import nan @@ -839,5 +880,5 @@ f_handle=pyfits.open(infile) img = f_handle[HDU].data - resid_img = 1e-7*ones(img.shape,Float32) + resid_img = zeros(img.shape,Float32) fitlist = [] for initialGuessTuple in initialGuessList: @@ -860,10 +901,10 @@ res = subim (scale,x,xsig,y,ysig,theta)=par - win.showArray(res) + # win.showArray(res) resid_img[int(y0)-halfwinsize:int(y0)+halfwinsize+1,int(x0)-halfwinsize:int(x0)+halfwinsize+1] = res # print initialGuessTuple, scale,x,xsig,y,ysig,theta fitlist.append(par) f_handle.close() - outname = 'gaussfit_resid_'+infile + outname = 'gaussfit_resid_r%ipix_%s' % (blankrad,infile) primHDU = pyfits.PrimaryHDU(resid_img) hdulist = pyfits.HDUList([primHDU]) @@ -871,5 +912,5 @@ return fitlist -def doEfbComparison(simfile='simtest.1.00a.fits',efbfile='efbsim.1.00a.fits',simtab='simtest.1.00a.cmf'): +def doEfbComparison(simfile='simtest.1.00a.fits',efbfile='efbsim.1.00a.fits',simtab='simtest.1.00a.cmf',blankrad=9): import pyfits from numpy import array, rec @@ -886,6 +927,6 @@ initlist.append((counts[i],x[i],2.,y[i],2.,0.)) inparlist.append((counts[i],x[i],y[i])) - efbparlist = doGaussPsfPhotometry(efbfile,initlist) - simparlist = doGaussPsfPhotometry(simfile,initlist) + efbparlist = doGaussPsfPhotometry(efbfile,initlist,blankrad=blankrad) + simparlist = doGaussPsfPhotometry(simfile,initlist,blankrad=blankrad) savelist = [ inparlist[i] + efbparlist[i] + simparlist[i] for i in range(nrows)] @@ -900,5 +941,5 @@ saveHDU = tabHDUfromRecArray(saverows) fitsend = re.compile('(\.[sc]mf|\.fits?)$') - saveHDU.writeto('gaussfit_'+fitsend.sub('',efbfile)+'___'+simfile,clobber=True) + saveHDU.writeto('gaussfit_r%ipix_%s___%s' %(blankrad,fitsend.sub('',efbfile),simfile),clobber=True) @@ -931,5 +972,5 @@ fpObjcDir += '/' # For testing: - # return ['1056-0192.421/1056-0192.421.ch.421.CHIP1.cmf'],['1056-0192.421/fpObjc-001056-1-0192.fit'] + return ['1056-0192.421/1056-0192.421.ch.421.CHIP1.cmf'],['1056-0192.421/fpObjc-001056-1-0192.fit'] from glob import glob