Index: /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py =================================================================== --- /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 20458) +++ /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 20459) @@ -58,4 +58,9 @@ # # Convention: for scatter, plot second against first; for histogram, plot both histograms + + +# XXX Todo: labels! group by runs? different colours for different +# bands; include seeing? plot things against seeing? + plotcol_1frame_tlist = [ (['d_sky'],['d_mag'],'scatter') @@ -63,5 +68,7 @@ ,(['sky_ps1'],['d_mag'],'scatter') ,(['psfinstmag_sdss'],['d_mag'],'scatter') + ,(['psfinstmag_sdss'],['d_sky'],'scatter') ,(['m_rr_cc_psf'],['d_mag'],'scatter') + ,(['sky_sdss'],['sky_ps1'],'scatter') ,(['psfinstmag_sdss'],['psf_inst_mag'],'histogram') ,(['sky_sdss'],['sky_ps1'],'histogram') @@ -75,7 +82,22 @@ ,(['d_x_median'],['d_y_median'],'scatter') ,(['d_mag_median'],['N_SDSS'],'scatter') - ,(['d_mag_median','d_mag_mean'],['field','field'],'scatter') + ,(['field','field'],['d_mag_median','d_mag_mean'],'scatter') ] +def mergePlots(plotcol_1frame_tlist=plotcol_1frame_tlist,filtlist=['u','g','r','i','z'],workdir='/IPP/data/SDSS/stripe82/coadd/compare'): + from subprocess import call + from glob import glob + import os + os.chdir(workdir) + for tup in plotcol_1frame_tlist: + t1,t2,t3=tup + for filt in filtlist: + globstr = getOutnameStatsOnefile("match_%s_*"%(filt),t3,t1,t2) + filelist = glob(globstr) + outname = getOutnameStatsOnefile("merged_%s" %(filt),t3,t1,t2) + cmdstr = "gs -dNOPAUSE -sDEVICE=pswrite -sOutputFile=%s %s -c quit" %(outname,' '.join(filelist)) + print "Making %s" %(outname) + call(cmdstr,shell=True) + call("gzip %s" %(outname),shell=True) def compIPPphoto(summaryTable,mode,plotcol_1frame_tlist=plotcol_1frame_tlist,\ @@ -109,6 +131,11 @@ stats_hash={} chipfile_l,fpObjc_l = makePlan() + filters = [] + bandindex_hash = {} for chipfile,fpObjc in zip(chipfile_l,fpObjc_l): matchtable,filter_name,bandindex = matchSdssPs1(fpObjc,chipfile,skip=skip) + if filter_name not in filters: + filters.append(filter_name) + bandindex_hash[filter_name] = bandindex stats_hash, column_hash, goodrow_hash = computeStatistics(matchtable) # Sort res_hash by its column names to make sure order is @@ -131,6 +158,10 @@ # Fake a goodrows hash sumgoodrows_hash = summaryhash.fromkeys(summaryhash.keys(),numpy.ones((nrows,),bool)) - plotStatsOnefile(summaryhash,sumgoodrows_hash,summaryTable,plotcol_summary_tlist,bandindex) + for filt in filters: + plotStatsOnefile(summaryhash,sumgoodrows_hash,summaryTable,plotcol_summary_tlist,bandindex_hash[filt],bandid=filt) + mergePlots(plotcol_1frame_tlist=plotcol_1frame_tlist,filtlist=filters) return column_hash,goodrow_hash + +# XXX Todo: check histogram axis scaling! def hashFromPyfitsTable(pyfitsTableHDU): @@ -144,8 +175,10 @@ return outhash -def getOutnameStatsOnefile(matchtable,kind,col1_l,col2_l=None,format='eps'): +def getOutnameStatsOnefile(matchtable,kind,col1_l,col2_l=None,format='eps',bandid=''): import re # Construct output filename root = re.sub('(\.[sc]mf|\.fits?)$','',matchtable) + if bandid != '': + root = "%s_%s" %(root,bandid) if isNone(col2_l): outname = '%s_%s_%s.%s' % (root,kind,col1[0],format) @@ -157,9 +190,10 @@ return outname -def plotStatsOnefile(values_hash,goodrow_hash,matchtable,plotcol_tlist,bandindex,format='eps'): +def plotStatsOnefile(values_hash,goodrow_hash,matchtable,plotcol_tlist,bandindex,format='eps',bandid=''): """Make diagnostic plots for a single table, based on values in values_hash""" from numpy import concatenate from sm import cvar,angle + import re # In histograms, only plot core of distribution within these percentiles: histo_min_ntile = 0.03 @@ -171,5 +205,5 @@ col2name_l = troika[1] plottype = troika[2] - outname = getOutnameStatsOnefile(matchtable,plottype,col1name_l,col2name_l,format=format) + outname = getOutnameStatsOnefile(matchtable,plottype,col1name_l,col2name_l,format=format,bandid=bandid) smOpenPlot(outname,format=format) firstplot = True @@ -238,4 +272,5 @@ def appendFitsTable(fitsfile,tabHDU,HDU=1): """Append rows in tabHDU to existing fits table file""" + import pyfits oldtabhdulist = pyfits.open(fitsfile) nrows_old = len(oldtabhdulist[HDU].data) @@ -370,6 +405,4 @@ ,'d_skyerr':['skyErr','SKY_SIGMA'] ,'d_pointsource':['prob_psf','pointsource_ps1'] - # ,'d_mag':['psfcounts','PSF_INST_MAG '] - # ,'d_magerr':['psfcountserr','PSF_INST_MAG_SIG'] ,'d_mag':['psfinstmag_sdss','PSF_INST_MAG'] ,'d_magerr':['psfinstmagerr_sdss','PSF_INST_MAG_SIG'] @@ -480,5 +513,5 @@ return array(good1_l),array(good2_l),array(good3_l) -def makePlan(): +def makePlan(lastcamcol=6): """Make paired list of which fpObjc file to compare to which IPP .cmf file. Could also think about constructing the fpObjc filename @@ -487,5 +520,5 @@ # 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 @@ -507,5 +540,5 @@ } Nfields = 18 - camcollist = range(1,6+1) + camcollist = range(1,lastcamcol+1) sdsslist = [] ps1list = [] @@ -567,4 +600,6 @@ sdssfilt += ';select !((flags[%s]&1<<1)!=0)&&(!((flags[%s]&1<<3)!=0)||((flags[%s]&1<<6)!=0)||nchild==0)' % \ (bandindex,bandindex,bandindex) + # add check for BINNED1 !!! bit28 in flags1 + # Check flag 0x4000 (extended) to compare against prob_psf in SDSS # - for some reason, 'flags' and the other integer columns get @@ -823,5 +858,4 @@ img = f_handle[HDU].data resid_img = 1e-7*ones(img.shape,Float32) - print resid_img.shape,type(resid_img),img.shape,type(img) fitlist = [] for initialGuessTuple in initialGuessList: @@ -846,5 +880,5 @@ 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 + # print initialGuessTuple, scale,x,xsig,y,ysig,theta fitlist.append(par) f_handle.close()