Index: /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py =================================================================== --- /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 19831) +++ /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 19832) @@ -35,122 +35,17 @@ # + Trends with field number, seeing etc. -def smdefaults(): - # Will these be persistent across multiple imports and going out - # of scope? Perhaps only if I do a global 'import sm'. Or it - # should work if called from something where sm is in scope? - sm.expand(1.8) - sm.lweight(5) - -# Plotting convention will be like in sm: you need to open a file, -# then you plot to it with one or multiple commands, then you close it -# to write the file. - -def smOpenPlot(filename,format='eps'): - """Issue sm device command for file of given format: - eps -> postfile - else just use given format as 'device'""" - from sm import device, erase - if filename == 'x11': - device('x11') - elif format == 'eps': - device('postfile '+filename) - else: - device(format +' '+filename) - erase() - -def smClosePlot(): - from sm import device - device('nodevice') - -def isNone(something): - return isinstance(something,type(None)) - -def smHistoPlot(vec,step=None,minbin=None,maxbin=None,nbins=None,\ - xlab=None,ylab=None,xrange=None,yrange=None,\ - box1=None,box2=None,box3=None,box4=None,\ - append=False): - """Plot a histogram, intelligently deriving bins from the given - parameters if they are given intelligently. Otherwise, silently - do nothing.""" - import sm - if isNone(minbin): - minbin = min(vec) - if isNone(maxbin): - maxbin = max(vec) - # I am assuming bins are bin centers - if not isNone(step): - nbins = (maxbin-minbin)/step+1 - if isNone(nbins): - return - histo,leftbinedges = histogram(vec,nbins,[minbin,maxbin]) - bincenters = leftbinedges + 0.5*(leftbinedges[1]-leftbinedges[0]) - - if not append: - smSetup(bincenters,histo,xrange,yrange,xlab,ylab,box1,box2,box3,box4) - sm.histogram(bincenters,histo) - -def smLinePlot(x,y,ltype=0,xlab=None,ylab=None,xrange=None,yrange=None,\ - box1=None,box2=None,box3=None,box4=None,\ - append=False): - """Make an sm scatter plot on current device. If append=True, - overplot with current limits. Otherwise, draw box box1 box2 box3 box4""" - import sm - sm.expand(1.8) - sm.lweight(5) - # Silently ignore any problems with plot generation - try: - if not append: - smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4) - sm.ltype(ltype) - sm.connect(x,y) - except: - pass - - -def smScatterPlot(x,y,ptype=41,xlab=None,ylab=None,xrange=None,yrange=None,\ - box1=None,box2=None,box3=None,box4=None,\ - append=False): - """Make an sm scatter plot on current device. If append=True, - overplot with current limits. Otherwise, draw box box1 box2 box3 box4""" - import sm - sm.expand(1.8) - sm.lweight(5) - # Silently ignore any problems with plot generation - try: - if not append: - smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4) - sm.ptype(ptype) - sm.points(x,y) - except: - pass - -def smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4): - import sm - sm.erase() - if isNone(xrange): - xrange = x - if isNone(yrange): - yrange = y - sm.limits(x,y) - smBox(box1,box2,box3,box4) - if not isNone(xlab): - sm.xlabel(xlab) - if not isNone(ylab): - sm.ylabel(ylab) - -def smBox(box1,box2,box3,box4): - from sm import box - if not isNone(box4): - box(box1,box2,box3,box4) - elif not isNone(box3): - box(box1,box2,box3) - elif not isNone(box2): - box(box1,box2) - elif not isNone(box1): - box(box1) - else: - box() - -def compIPPphoto(summaryTable,mode): +# Big (BIG) XXX: For plotting, need to make sure that ALL columns have +# "good" data in exactly the same rows, otherwise lose correspondence +# between them. E.g.: for every column, create a logical vector saying +# whether it's 'good' or not, and then plot things only for those rows +# where both are good + +plotcol_tlist = [ + ('d_mag','d_sky','scatter'), + ('d_x','d_y','scatter') + ] + + +def compIPPphoto(summaryTable,mode,plotcol_tlist=plotcol_tlist): """summaryTable: .fits table for output mode: new or append for creating summaryTable new or appending current run's output to it. @@ -178,12 +73,15 @@ 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 identical in all - # rows of the table, and output is more legible + res_hash, deltas_hash = computeStatistics(matchtable) + # Sort res_hash by its 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) + plotStatsOnefile(deltas_hash,matchtable,plotcol_tlist) newrows = numpy.rec.array(rowtuple_list,names=keylist) tabhdu = tabHDUfromRecArray(newrows) @@ -192,4 +90,27 @@ else: appendFitsTable(summaryTable,tabhdu) + +def getOutnameStatsOnefile(matchtable,kind,col1,col2=None,format='eps'): + import re + # Construct output filename + root = re.sub('(\.[sc]mf|\.fits?)$','',matchtable) + if isNone(col2): + outname = '%s_%s_%s.%s' % (root,kind,col1,format) + else: + outname = '%s_%s_%s_%s.%s' % (root,kind,col1,col2,format) + return outname + +def plotStatsOnefile(deltas_hash,matchtable,plotcol_tlist,format='eps'): + """Make diagnostic plots for a single table, based on values in deltas_hash""" + for troika in plotcol_tlist: + col1name = troika[0] + col2name = troika[1] + plottype = troika[2] + outname = getOutnameStatsOnefile(matchtable,plottype,col1name,col2name,format=format) + smOpenPlot(outname,format=format) + if plottype == 'scatter': + smScatterPlot(deltas_hash[col1name],deltas_hash[col2name],\ + xlab=col1name,ylab=col2name) + smClosePlot() def valuesKeysSortedByKeys(hash): @@ -301,5 +222,5 @@ ,'d_magerr':['psfcountserr','PSF_INST_MAG_SIG'] } - + colval_hash = {} ismag = re.compile('mag') iscounts = re.compile('counts') @@ -334,6 +255,6 @@ 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 + colval_hash[outcol] = delta avg = delta.mean() lowq,med,upq = stats_med(delta) @@ -351,5 +272,5 @@ writeTable(tablename,newprimhdu,newtab) infile_handle.close() - return outhash + return outhash,colval_hash def writeTable(filename,primhdu,tabhdu): @@ -450,6 +371,20 @@ 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""" + """Call matchByPos to match an SDSS fpObjc.fits and a PS1 bla.cmf + file.""" import pyfits + + def getOutname(SDSSfile,PS1file,sdssbandstr): + import re + # Construct output filename + fitsend = re.compile('(\.[sc]mf|\.fits?)$') + 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) + return "match_%s___%s___%s.fits" % (sdssbandstr,SDSSroot,PS1root) + filters = ['u','g','r','i','z'] # Read primary header of PS1 file to work out band @@ -480,5 +415,5 @@ ps1pos = "\'x_psf y_psf\'" - outname = getOutname(SDSSfpObjc,PS1cmf,sdssbandstr) + outname = getOutname(SDSSfpObjc,PS1cmf,sdssbandstr) matchByPos(SDSSfpObjc,PS1cmf,outname,sdsspos,ps1pos,tolerance=matchrad,\ duptag1='_sdss',duptag2='_ps1',\ @@ -511,17 +446,4 @@ return outhash -def getOutname(SDSSfile,PS1file,sdssbandstr): - import re - # Construct output filename - fitsend = re.compile('(\.cm[sf]|\.fits?)$') - 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) - 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=""): """Match two fits tables by cartesian 2-d position using stilts/tmatch2 from topcat""" @@ -542,2 +464,121 @@ print tmatch_cmd retval = subprocess.call(tmatch_cmd,shell=True) + +def smdefaults(): + # Will these be persistent across multiple imports and going out + # of scope? Perhaps only if I do a global 'import sm'. Or it + # should work if called from something where sm is in scope? + sm.expand(1.8) + sm.lweight(5) + +# Plotting convention will be like in sm: you need to open a file, +# then you plot to it with one or multiple commands, then you close it +# to write the file. + +def smOpenPlot(filename,format='eps'): + """Issue sm device command for file of given format: + eps -> postfile + else just use given format as 'device'""" + from sm import device, erase + if filename == 'x11': + device('x11') + elif format == 'eps': + device('postfile '+filename) + else: + device(format +' '+filename) + erase() + +def smClosePlot(): + from sm import device + device('nodevice') + +def isNone(something): + return isinstance(something,type(None)) + +def smHistoPlot(vec,step=None,minbin=None,maxbin=None,nbins=None,\ + xlab=None,ylab=None,xrange=None,yrange=None,\ + box1=None,box2=None,box3=None,box4=None,\ + append=False): + """Plot a histogram, intelligently deriving bins from the given + parameters if they are given intelligently. Otherwise, silently + do nothing.""" + import sm + if isNone(minbin): + minbin = min(vec) + if isNone(maxbin): + maxbin = max(vec) + # I am assuming bins are bin centers + if not isNone(step): + nbins = (maxbin-minbin)/step+1 + if isNone(nbins): + return + histo,leftbinedges = histogram(vec,nbins,[minbin,maxbin]) + bincenters = leftbinedges + 0.5*(leftbinedges[1]-leftbinedges[0]) + + if not append: + smSetup(bincenters,histo,xrange,yrange,xlab,ylab,box1,box2,box3,box4) + sm.histogram(bincenters,histo) + +def smLinePlot(x,y,ltype=0,xlab=None,ylab=None,xrange=None,yrange=None,\ + box1=None,box2=None,box3=None,box4=None,\ + append=False): + """Make an sm scatter plot on current device. If append=True, + overplot with current limits. Otherwise, draw box box1 box2 box3 box4""" + import sm + sm.expand(1.8) + sm.lweight(5) + # Silently ignore any problems with plot generation + try: + if not append: + smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4) + sm.ltype(ltype) + sm.connect(x,y) + except: + pass + + +def smScatterPlot(x,y,ptype=41,xlab=None,ylab=None,xrange=None,yrange=None,\ + box1=None,box2=None,box3=None,box4=None,\ + append=False): + """Make an sm scatter plot on current device. If append=True, + overplot with current limits. Otherwise, draw box box1 box2 box3 box4""" + import sm + sm.expand(1.8) + sm.lweight(5) + # Silently ignore any problems with plot generation + try: + if not append: + smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4) + sm.ptype(ptype) + sm.points(x,y) + except: + pass + +def smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4): + import sm + sm.erase() + if isNone(xrange): + xrange = x + if isNone(yrange): + yrange = y + sm.limits(x,y) + smBox(box1,box2,box3,box4) + if not isNone(xlab): + sm.xlabel(xlab) + if not isNone(ylab): + sm.ylabel(ylab) + +def smBox(box1,box2,box3,box4): + from sm import box + if not isNone(box4): + box(box1,box2,box3,box4) + elif not isNone(box3): + box(box1,box2,box3) + elif not isNone(box2): + box(box1,box2) + elif not isNone(box1): + box(box1) + else: + box() + +