Index: /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py =================================================================== --- /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 19787) +++ /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 19787) @@ -0,0 +1,382 @@ +# Routines for comparing IPP outputs against photometry from other +# sources. +# +# Currently only SDSS's photo, but meant to extend to others including +# simtest +# +# Author: Sebastian Jester jester at mpia.de +# +# Depends on: +# +# topcat with 'stilts' shell script from +# http://www.star.bris.ac.uk/~mbt/topcat/ and +# +# pyfits from http://www.stsci.edu/resources/software_hardware/pyfits +# (which in turn has dependencies) +# +# + +def compIPPphoto(summaryTable,mode): + """summaryTable: .fits table for output + mode: new or append for creating summaryTable new or appending current run's output to it. + + Steps: + + 1. read "plan file" with list of files to be looked at, or otherwise + generate that list, e.g. by globbing (paired list fpObjc <-> cmf/smf) + + for every file: + 2. match the files + + 3. Compute diagnostic columns + + 4. Compute statistics on diagnostic columns (begin with mean, rms, quartiles) + + 5. Append summary statistics to master table + + 6. Generate whatever plots are desired from indidvidual match tables + or master table + """ + 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 = [] + 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]) + rowtuple_list.append(vallist) + newrows = numpy.rec.array(rowtuple_list,names=keylist) + tabhdu = tabHDUfromRecArray(newrows) + if mode == 'new': + 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) + + +def tabHDUfromRecArray(recarr): + """Generate a table HDU from a record array""" + # create column names from recarray._names + # Assume float as data type + # copy data field-wise to table HDU. + import pyfits + collist = [] + for colname in recarr.dtype.names: + if isinstance(recarr.field(colname)[0],int): + format = '1J' + elif isinstance(recarr.field(colname)[0],str): + # This is potentially dangerous, but for now, the only + # string column I am writing is the filter, i.e. 1A should + # be enough already. + format = '1A' + else: + format = '1E' + collist.append(pyfits.Column(name=colname,array=recarr.field(colname),format=format)) + coldefs = pyfits.ColDefs(collist) + return pyfits.new_table(coldefs) + +def mkRecordArray(hash): + """Return a record array representing a table row from a hash; + keys are column names, values are data""" + newrow = numarray.records.array([tuple(hash.values())]) + newrow.dtype.names=hash.keys() + # Not : + # newrow = numpy.rec.array([tuple(hash.values())],names=hash.keys())? + return newrow + +def saveHierarchHeaderList(header,label_l,value_l): + """Given paired lists label_l and value_l, save the values from + value_l to ESO HIERARCH header keywords with the names given in + label_l""" + + for label,value in zip(label_l,value_l): + header.update('HIERARCH %s'%(label),value) + +def stats_med(array): + from numpy import nan,sort + sortarray = sort(array,kind='mergesort') + l = len(array) + if l == 0: + return nan,nan,nan + lowerquartile = sortarray[0.25*l] + median = sortarray[0.5*l] + upperquartile= sortarray[0.75*l] + return lowerquartile,median,upperquartile + +def computeStatistics(tablename): + """Compute desired statistics for offsets etc. and add them as + ESO-style HIERARCH fields to the table header (8 characters are + otherwise too few to contain meaningful names for the + keywords). Returns a hash with all the computed statistics, keyed + by their header names""" + import pyfits,re,operator + from numpy import sqrt,log10,log + + label_l = ['mean','rms','median','lowerquartile','upperquartile'] + # Hash of output header keyword names pointing to paired lists of + # columns that are to be subtracted from each other for statistics + # to be generated. + filterID={'u':0,'g':1,'r':2,'i':3,'z':4} + # For this to work, table needs to be called match_r_... for r-band + + f = pyfits.open(tablename,mode='update') + table = f[1] + h = table.header + table_data = table.data + run = h['RUN'] + rerun = h['RERUN'] + camcol = h['CAMCOL'] + field = h['FIELD'] + filtnam = h['FILTER'] + + # These are just the columns; need to get a slice with the correct array index later + colname_hash = { + 'd_x':['colc','x_psf'] + ,'d_y':['rowc','y_psf'] + ,'d_mag':['psfcounts','PSF_INST_MAG'] + ,'d_magerr':['psfcountserr','PSF_INST_MAG_SIG'] + ,'d_sky':['sky_sdss','SKY_ps1'] + ,'d_skyerr':['skyErr','SKY_SIGMA'] + ,'d_pointsource':['prob_psf','pointsource_ps1'] + } + + outhash = {'RUN':run,'RERUN':rerun,'CAMCOL':camcol,'FIELD':field,'FILTER':filtnam} + ismag = re.compile('mag') + iscounts = re.compile('counts') + outcoll = [] + for outcol,list in colname_hash.iteritems(): + SDSScolname,PS1colname=list + # Slice out the proper filter - may need to add another + # boolean dictionary that says whether or not this is + # necessary for a given column (e.g. for objc_colc it won't + # be...) + SDSScol = table_data.field(SDSScolname)[:,filterID[filtnam]] + PS1col = table_data.field(PS1colname) + SDSScol_good,PS1col_good = filterGoodVal(SDSScol,PS1col) + # Compute SDSS instrumental magnitude if necessary + if PS1colname == 'PSF_INST_MAG' and SDSScolname.lower() == 'psfcounts': + # Create column with SDSS instrumental magnitude that can + # be written to table for diagnostic use + instmag_sdss_arr = -2.5*log10(table_data.field(SDSScolname)) + instmag_col = pyfits.Column(name='psfinstmag_sdss',format='5E',array=instmag_sdss_arr) + outcoll.append(instmag_col) + # Compute array for internal use + SDSScol_good = -2.5*log10(SDSScol_good) + if PS1colname == 'PSF_INST_MAG_SIG' and SDSScolname.lower() == 'psfcountserr': + # Create column with SDSS instrumental magnitude error + # that can be written to table for diagnostic use + instmagerr_sdss_arr = 2.5/log(10.)*table_data.field('psfcountserr')/table_data.field('psfcounts') + instmagerr_col = pyfits.Column(name='psfinstmagerr_sdss',format='5E',array=instmagerr_sdss_arr) + outcoll.append(instmagerr_col) + # Compute array for internal use + SDSScounts = table_data.field('psfcounts')[:,filterID[filtnam]] + SDSScol_good,PS1col_good,SDSScounts_good = filterGoodVal3(SDSScol,PS1col,SDSScounts) + SDSScol_good = 2.5/log(10.)*SDSScol_good/SDSScounts_good + delta = SDSScol_good - PS1col_good + avg = delta.mean() + lowq,med,upq = stats_med(delta) + rms = sqrt(delta.var()) + # Save to fits columns/header + outheaderlabel = [operator.add(outcol+' ',itm).upper() for itm in label_l] + outtablabel = [operator.add(outcol+'_',itm).upper() for itm in label_l] + saveHierarchHeaderList(h,outheaderlabel,[avg,rms,med,lowq,upq]) + # and in return variable + 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) + writeTable(tablename,newprimhdu,newtab) + return outhash + +def writeTable(filename,primhdu,tabhdu): + """Write pyfits.ColDefs thing columns to filename, clobbering + existing file""" + import pyfits + hdulst = pyfits.HDUList([primhdu,tabhdu]) + hdulst.writeto(filename,clobber=True) + +def filterGoodVal(col1,col2): + """Return a pair of arrays that contains those items where both + col1 and col2 have NaNs and SDSS "no-value" and + "no-error" flags -9999 and -1000 filtered out.""" + from numpy import isfinite,logical_and,array + # I don't understand why this doesn't work... + # + #goodcondition = logical_and(isfinite(col1),isfinite(col2)) + #col1 = col1[goodcondition] + #col2_good = col2[goodcondition] + #return col1,col2 + # + #... so I'm doing this: + good1_l = [] + good2_l = [] + for v1,v2 in zip(col1,col2): + if isfinite(v1) and isfinite(v2) and v1 not in [-9999,-1000] and v2 not in [-9999,-1000]: + good1_l.append(v1) + good2_l.append(v2) + return array(good1_l),array(good2_l) + +def filterGoodVal3(col1,col2,col3): + """Return a triplet of arrays that contains those items where both + col1 and col2 have NaNs and SDSS "no-value" and + "no-error" flags -9999 and -1000 filtered out.""" + from numpy import isfinite,array + # + #... so I'm doing this: + good1_l = [] + good2_l = [] + good3_l = [] + for i in range(len(col1)): + v1 = col1[i] + v2 = col2[i] + v3 = col3[i] + if isfinite(v1) and isfinite(v2) and isfinite(v3) and \ + v1 not in [-9999,-1000] and v2 not in [-9999,-1000] and \ + v3 not in [-9999,-1000]: + good1_l.append(v1) + good2_l.append(v2) + good3_l.append(v3) + return array(good1_l),array(good2_l),array(good3_l) + +def makePlan(): + """Make paired list of which fpObjc file to compare to which IPP + .cmf file. Could also think about constructing the fpObjc filename + from run, camcol, field read from .cmf primary header, which + replicates most of the fpC primary header.""" + + from glob import glob + fpObjcDir = '/IPP/data/SDSS/stripe82/coadd/input/fpObjcs/' + cmfDir = '/IPP/data/SDSS/stripe82/coadd/prod/run_ipp/' + runlist = [1056,1755,3388,3434,3465,4145,4192,4203,4247,4263,5052] + firstfield_h = { + 1056:192, + 1755:329, + 3388:273, + 3434:306, + 3465:185, + 4145:63, + 4192:309, + 4203:358, + 4247:62, + 4263:69, + 5052:62 + } + 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 = [] + for run in runlist[0:1]: + firstfield = firstfield_h[run] + for field in range(firstfield,firstfield+Nfields+1): + for camcol in camcollist: + sdssname = fpObjcDir + "fpObjc-%06i-%i-%04i.fit" % ( run, camcol, field) + # The following should just return 5 files (1 per + # filter), but just to be sure, let's check and return + # the correct number of sdss files to match: + ps1glob = glob(cmfDir + "%i-%04i.*/%i-%04i.*.ch.*.CHIP%i.cmf" %(run,field,run,field,camcol)) + for ps1name in ps1glob: + ps1list.append(ps1name) + sdsslist.append(sdssname) + # Need to remember that CHIPx corresponds to CAMCOLx; could change + # naming rule to make filter name appear in filename, too. + return ps1list,sdsslist + +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 + 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 + bandindex = filters.index(sdssbandstr) + run = h['RUN'] + rerun = h['RERUN'] + camcol=h['CAMCOL'] + field = h['FIELD'] + sdssf.close() + # 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 + # So subtract 0.5 from SDSS numbers and match with r=0.7 should include all good matches + # Arghly, need to match datatypes and 0.5 is double by default + # while cols are float, hence the typecast. + sdssfilt = 'addcol colc_%s colc[%s]-(float)%s' % (sdssbandstr,bandindex,xoff) + sdssfilt += ';addcol rowc_%s rowc[%s]-(float)%s' % (sdssbandstr,bandindex,yoff) +# sdssfilt += ';addcol psfinstmag_sdss -2.5*log10(psfcounts)' +# sdssfilt += ';addcol psfinstmagerr_sdss 2.5/ln(10)*psfcountserr/psfcounts' + # Filtering out objects that can't be primary, i.e. + # !BRIGHT && (!BLENDED || NODEBLEND || nchild == 0) + sdssfilt += ';select !((flags[%s]&1<<1)!=0)&&(!((flags[%s]&1<<3)!=0)||((flags[%s]&1<<6)!=0)||nchild==0)' % \ + (bandindex,bandindex,bandindex) + # Check flag 0x4000 (extended) to compare against prob_psf in SDSS + # - for some reason, 'flags' and the other integer columns get + # read as float by topcat, so need a cast here + ps1filt = 'addcol pointsource_ps1 (((int)flags&0x4000)==0?1:0)' + + sdsspos = '\'colc_%s rowc_%s\''% (sdssbandstr,sdssbandstr) + ps1pos = "\'x_psf y_psf\'" + # Construct output filename + fitsend = re.compile('(\.cm[sf]|\.fits?)$') + SDSSroot = fitsend.sub('',SDSSfpObjc) + PS1root = fitsend.sub('',PS1cmf) + # 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 + +def matchByPos(in1,in2,out,colnames1,colnames2,tolerance=0.5,duptag1='_sdss',duptag2='_ps1',filter1="",filter2=""): + """Match two fits tables by cartesian 2-d position using stilts/tmatch2 from topcat""" + import subprocess + # tmatch_args - sequence of parameters is + # infile1, infile2, outfile, x_col, y_col, tolerance [all in pixels], duptag1,2 + # Looking for best match only + tmatch_args = 'in1=%s in2=%s out=%s values1=%s values2=%s params=%s fixcols=dups suffix1=%s suffix2=%s' + tmatch_args += ' matcher=2d ifmt1=fits ifmt2=fits ofmt=fits-basic join=1or2 find=best' + if filter1 != "": + tmatch_args += " icmd1=%s" %(filter1) + if filter2 != "": + tmatch_args += " icmd2=%s" %(filter2) + + tmatch_cmd = tmatch_args % (in1,in2,out,colnames1,colnames2,tolerance,duptag1,duptag2) + tmatch_cmd = os.environ.get('HOME')+'/bin/stilts tmatch2 '+tmatch_cmd + + print tmatch_cmd + retval = subprocess.call(tmatch_cmd,shell=True)