Changeset 20459

Timestamp:

Oct 29, 2008, 6:35:57 AM (18 years ago)

Author:

Sebastian Jester

Message:

Add plot-merging function

File:

• branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (modified) (16 diffs)

branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py

@@ -58 +58 @@
 #
 # 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 +68 @@
     ,(['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 +82 @@
     ,(['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 +131 @@
     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 +158 @@
     # 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 +175 @@
     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 +190 @@
     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 +205 @@
         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 +272 @@
 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 +405 @@
         ,'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 +513 @@
     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 +520 @@
      
     # 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 +540 @@
         }
     Nfields = 18
-    camcollist = range(1,6+1)
+    camcollist = range(1,lastcamcol+1)
     sdsslist = []
     ps1list = []
@@ -567 +600 @@
     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 +858 @@
     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 +880 @@
         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()