Index: /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py =================================================================== --- /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 19800) +++ /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 19801) @@ -15,5 +15,23 @@ # (which in turn has dependencies) # -# +# Also see: http://kiawe.ifa.hawaii.edu/IPPwiki/index.php/SimtestVerification +# +# Todo: +# * rms etc. with outlier rejection - e.g. by histogram core fit. +# * Plots +# - Per-field: +# + d_x vs. d_y +# + magnitude histograms +# + flux vs. sky +# + magdiff vs. mag +# + magdiff vs. sky +# + magdiff vs. skydiff +# + something about PSF sizes; use Michigan moments for photo - +# M_rr_cc_psf is the size of the reconstructed PSF, and m_rr_cc is the size of the object +# + magdiff, mag vs. d_x^2+d_y^2 +# + histogram of chi - but meaningful only for repeat measurements? +# - per-run: +# + histograms of everything - recovery fractions, means, medians, quartiles +# + Trends with field number, seeing etc. def compIPPphoto(summaryTable,mode): @@ -39,21 +57,15 @@ """ import pyfits,numpy - - - chipfile_l,fpObjc_l = makePlan() - if mode.lower() not in ["new","append"]: raise ValueError("Must specify mode=new or mode=append") rowtuple_list = [] + + chipfile_l,fpObjc_l = makePlan() for chipfile,fpObjc in zip(chipfile_l,fpObjc_l): matchtable = matchSdssPs1(fpObjc,chipfile) res_hash = computeStatistics(matchtable) - # Sort by column names to make sure order is correct, and - # output is more legible - vallist = [] - keylist = res_hash.keys() - keylist.sort() - for field in keylist: - vallist.append(res_hash[field]) + # Sort by column names to make sure order is identical in all + # rows of the table, and output is more legible + vallist,keylist = valuesKeysSortedByKeys(res_hash) rowtuple_list.append(vallist) newrows = numpy.rec.array(rowtuple_list,names=keylist) @@ -62,13 +74,26 @@ tabhdu.writeto(summaryTable,clobber=True) else: - oldtabhdulist = pyfits.open(summaryTable) - nrows_old = len(oldtabhdulist[1].data) - newtabhdu = pyfits.new_table(oldtabhdulist[1].columns,nrows = nrows_old+len(tabhdu.data)) - for col in newtabhdu.columns.names: - newtabhdu.data.field(col)[nrows_old:] = tabhdu.data.field(col) - oldtabhdulist.close() - newtabhdu.writeto(summaryTable,clobber=True) - - + appendFitsTable(summaryTable,tabhdu) + +def valuesKeysSortedByKeys(hash): + """Return a list of hash's values and keys that is ordered by the names of the keys""" + vallist = [] + keylist = hash.keys() + keylist.sort() + for field in keylist: + vallist.append(hash[field]) + return vallist,keylist + +def appendFitsTable(fitsfile,tabHDU,HDU=1): + """Append rows in tabHDU to existing fits table file""" + oldtabhdulist = pyfits.open(fitsfile) + nrows_old = len(oldtabhdulist[HDU].data) + newtabhdu = pyfits.new_table(oldtabhdulist[HDU].columns,nrows = nrows_old+len(tabHDU.data)) + for col in newtabhdu.columns.names: + newtabhdu.data.field(col)[nrows_old:] = tabHDU.data.field(col) + oldtabhdulist.close() + return newtabhdu.writeto(fitsfile,clobber=True) + + def tabHDUfromRecArray(recarr): """Generate a table HDU from a record array""" @@ -84,5 +109,6 @@ # This is potentially dangerous, but for now, the only # string column I am writing is the filter, i.e. 1A should - # be enough already. + # be enough already. Perhaps use the length of the entry + # as column format? format = '1A' else: @@ -120,5 +146,5 @@ return lowerquartile,median,upperquartile -def computeStatistics(tablename): +def computeStatistics(tablename,copyfields_list = ['RUN','RERUN','CAMCOL','FIELD','FILTER','FWHM_X','FWHM_Y']): """Compute desired statistics for offsets etc. and add them as ESO-style HIERARCH fields to the table header (8 characters are @@ -136,14 +162,15 @@ # For this to work, table needs to be called match_r_... for r-band - f = pyfits.open(tablename,mode='update') - table = f[1] + infile_handle = pyfits.open(tablename,mode='update') + table = infile_handle[1] h = table.header table_data = table.data # Header names that are going to be copied to output table as # "reference columns" - copyfields_list = ['RUN','RERUN','CAMCOL','FIELD','FILTER'] outhash = {} for f in copyfields_list: - outhash[f] = h[f] + outhash[f] = h[f] + + filtname = outhash['FILTER'] # These are just the columns; need to get a slice with the correct array index later @@ -167,5 +194,5 @@ # necessary for a given column (e.g. for objc_colc it won't # be...) - SDSScol = table_data.field(SDSScolname)[:,filterID[filtnam]] + SDSScol = table_data.field(SDSScolname)[:,filterID[filtname]] PS1col = table_data.field(PS1colname) SDSScol_good,PS1col_good = filterGoodVal(SDSScol,PS1col) @@ -186,5 +213,5 @@ outcoll.append(instmagerr_col) # Compute array for internal use - SDSScounts = table_data.field('psfcounts')[:,filterID[filtnam]] + SDSScounts = table_data.field('psfcounts')[:,filterID[filtname]] SDSScol_good,PS1col_good,SDSScounts_good = filterGoodVal3(SDSScol,PS1col,SDSScounts) SDSScol_good = 2.5/log(10.)*SDSScol_good/SDSScounts_good @@ -200,9 +227,9 @@ for label,value in zip(outtablabel,[avg,rms,med,lowq,upq]): outhash[label]=value - f.close() newtab = pyfits.new_table(table.columns+pyfits.ColDefs(outcoll),header=h) - newprimhdu = pyfits.PrimaryHDU(header=f[0].header) + newprimhdu = pyfits.PrimaryHDU(header=infile_handle[0].header) writeTable(tablename,newprimhdu,newtab) + infile_handle.close() return outhash @@ -262,5 +289,8 @@ from run, camcol, field read from .cmf primary header, which replicates most of the fpC primary header.""" - + + # For testing: + return ['1056-0192.421/1056-0192.421.ch.421.CHIP1.cmf'],['1056-0192.421/fpObjc-001056-1-0192.fit'] + from glob import glob fpObjcDir = '/IPP/data/SDSS/stripe82/coadd/input/fpObjcs/' @@ -282,6 +312,4 @@ Nfields = 18 camcollist = range(1,6+1) - # For testing: - # pairlist = ['1056-0192.421/1056-0192.421.ch.421.CHIP1.cmf'],['1056-0192.421/fpObjc-001056-1-0192.fit'] sdsslist = [] ps1list = [] @@ -304,20 +332,12 @@ def matchSdssPs1(SDSSfpObjc,PS1cmf,xoff=0.5,yoff=0.5,matchrad=0.7): """Call matchByPos to match an SDSS fpObjc.fits and a PS1 bla.cmf file""" - import pyfits,re + import pyfits filters = ['u','g','r','i','z'] # Read primary header of PS1 file to work out band - ps1f = pyfits.open(PS1cmf) - h=ps1f[0].header - sdssbandstr = h['FILTER'] - ps1f.close() - - sdssf = pyfits.open(SDSSfpObjc) - h=sdssf[0].header + ps1copyfields_hash = headerfieldHash(['FWHM_X','FWHM_Y','FILTER'],PS1cmf,HDU=0) + sdssbandstr = ps1copyfields_hash['FILTER'] bandindex = filters.index(sdssbandstr) - run = h['RUN'] - rerun = h['RERUN'] - camcol=h['CAMCOL'] - field = h['FIELD'] - sdssf.close() + sdsscopyfields_hash = headerfieldHash(['RUN','RERUN','CAMCOL','FIELD'],SDSSfpObjc,0) + # Note position offset - y_sdss-y_psf1 ~ 0.6, x_sdss-x_ps1 ~ 0.35 (in 1056-r1-0192) # width of distribution is ~ 0.4 @@ -340,26 +360,48 @@ sdsspos = '\'colc_%s rowc_%s\''% (sdssbandstr,sdssbandstr) ps1pos = "\'x_psf y_psf\'" + + outname = getOutname(SDSSfpObjc,PS1cmf,sdssbandstr) + matchByPos(SDSSfpObjc,PS1cmf,outname,sdsspos,ps1pos,tolerance=matchrad,\ + duptag1='_sdss',duptag2='_ps1',\ + filter1='\''+sdssfilt+'\'',filter2='\''+ps1filt+'\'') + + f=pyfits.open(outname,mode='update') + h=f[1].header + updateHeaderFromHash(h,sdsscopyfields_hash) + updateHeaderFromHash(h,ps1copyfields_hash) + f.close() + + return outname + +def updateHeaderFromHash(header,hash): + """Add fields in hash to header""" + for k in hash.keys(): + header.update(k,hash[k]) + + +def headerfieldHash(fieldlist,file,HDU=0): + """Return a hash containing listed fits headers from file's HDU + (this is like a hash slice, which doesn't seem to be possible in + python?)""" + import pyfits + f=pyfits.open(file) + h=f[HDU].header + outhash = {} + for field in fieldlist: + outhash[field] = h[field] + return outhash + +def getOutname(SDSSfile,PS1file,sdssbandstr): + import re # Construct output filename fitsend = re.compile('(\.cm[sf]|\.fits?)$') - SDSSroot = fitsend.sub('',SDSSfpObjc) - PS1root = fitsend.sub('',PS1cmf) + SDSSroot = fitsend.sub('',SDSSfile) + PS1root = fitsend.sub('',PS1file) # Now chop off any leading path components pathroot = re.compile('.*/') SDSSroot = pathroot.sub('',SDSSroot ) PS1root = pathroot.sub('',PS1root) - outname = "match_%s___%s___%s.fits" % (sdssbandstr,SDSSroot,PS1root) - matchByPos(SDSSfpObjc,PS1cmf,outname,sdsspos,ps1pos,tolerance=matchrad,\ - duptag1='_sdss',duptag2='_ps1',\ - filter1='\''+sdssfilt+'\'',filter2='\''+ps1filt+'\'') - f=pyfits.open(outname,mode='update') - h=f[1].header - h.update('FILTER',sdssbandstr) - h.update('RUN',run) - h.update('RERUN',rerun) - h.update('CAMCOL',camcol) - h.update('FIELD',field) - f.close() - - return outname + return "match_%s___%s___%s.fits" % (sdssbandstr,SDSSroot,PS1root) + def matchByPos(in1,in2,out,colnames1,colnames2,tolerance=0.5,duptag1='_sdss',duptag2='_ps1',filter1="",filter2=""):