IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 19832


Ignore:
Timestamp:
Oct 2, 2008, 10:13:58 AM (18 years ago)
Author:
Sebastian Jester
Message:

Very first working version of flexible plot specification - give it a plotcol_tlist, i.e. a list of tuples of plotting columns, with (column1,column2,type of plot) and it will plot them all. But realized that I need to drag all the columns around from the input table for later plotting (i.e. perhaps do plotting inside statistics computation function?), and need to add bookkeeping for which rows are 'good' in ALL the columns.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py

    r19831 r19832  
    3535#     + Trends with field number, seeing etc.
    3636
    37 def smdefaults():
    38     # Will these be persistent across multiple imports and going out
    39     # of scope? Perhaps only if I do a global 'import sm'. Or it
    40     # should work if called from something where sm is in scope?
    41     sm.expand(1.8)
    42     sm.lweight(5)
    43 
    44 # Plotting convention will be like in sm: you need to open a file,
    45 # then you plot to it with one or multiple commands, then you close it
    46 # to write the file.
    47 
    48 def smOpenPlot(filename,format='eps'):
    49     """Issue sm device command for file of given format:
    50     eps -> postfile
    51     else just use given format as 'device'"""
    52     from sm import device, erase
    53     if filename == 'x11':
    54         device('x11')
    55     elif format == 'eps':
    56         device('postfile '+filename)
    57     else:
    58         device(format +' '+filename)   
    59     erase()
    60    
    61 def smClosePlot():
    62     from sm import device
    63     device('nodevice')
    64 
    65 def isNone(something):
    66     return isinstance(something,type(None))
    67 
    68 def smHistoPlot(vec,step=None,minbin=None,maxbin=None,nbins=None,\
    69                     xlab=None,ylab=None,xrange=None,yrange=None,\
    70                     box1=None,box2=None,box3=None,box4=None,\
    71                     append=False):
    72     """Plot a histogram, intelligently deriving bins from the given
    73     parameters if they are given intelligently.  Otherwise, silently
    74     do nothing."""
    75     import sm
    76     if isNone(minbin):
    77         minbin = min(vec)
    78     if isNone(maxbin):
    79         maxbin = max(vec)
    80     # I am assuming bins are bin centers
    81     if not isNone(step):
    82         nbins = (maxbin-minbin)/step+1
    83     if isNone(nbins):
    84         return
    85     histo,leftbinedges = histogram(vec,nbins,[minbin,maxbin])
    86     bincenters = leftbinedges + 0.5*(leftbinedges[1]-leftbinedges[0])
    87 
    88     if not append:
    89         smSetup(bincenters,histo,xrange,yrange,xlab,ylab,box1,box2,box3,box4)
    90     sm.histogram(bincenters,histo)
    91 
    92 def smLinePlot(x,y,ltype=0,xlab=None,ylab=None,xrange=None,yrange=None,\
    93                     box1=None,box2=None,box3=None,box4=None,\
    94                     append=False):
    95     """Make an sm scatter plot on current device. If append=True,
    96     overplot with current limits. Otherwise, draw box box1 box2 box3 box4"""
    97     import sm
    98     sm.expand(1.8)
    99     sm.lweight(5)
    100     # Silently ignore any problems with plot generation
    101     try:
    102         if not append:
    103             smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4)
    104         sm.ltype(ltype)
    105         sm.connect(x,y)
    106     except:
    107         pass
    108 
    109 
    110 def smScatterPlot(x,y,ptype=41,xlab=None,ylab=None,xrange=None,yrange=None,\
    111                     box1=None,box2=None,box3=None,box4=None,\
    112                     append=False):
    113     """Make an sm scatter plot on current device. If append=True,
    114     overplot with current limits. Otherwise, draw box box1 box2 box3 box4"""
    115     import sm
    116     sm.expand(1.8)
    117     sm.lweight(5)
    118     # Silently ignore any problems with plot generation
    119     try:
    120         if not append:
    121             smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4)
    122         sm.ptype(ptype)
    123         sm.points(x,y)
    124     except:
    125         pass
    126 
    127 def smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4):
    128     import sm
    129     sm.erase()
    130     if isNone(xrange):
    131         xrange = x
    132     if isNone(yrange):
    133         yrange = y
    134     sm.limits(x,y)
    135     smBox(box1,box2,box3,box4)
    136     if not isNone(xlab):
    137         sm.xlabel(xlab)
    138     if not isNone(ylab):
    139         sm.ylabel(ylab)
    140 
    141 def smBox(box1,box2,box3,box4):
    142     from sm import box
    143     if not isNone(box4):
    144         box(box1,box2,box3,box4)
    145     elif not isNone(box3):
    146         box(box1,box2,box3)
    147     elif not isNone(box2):
    148         box(box1,box2)
    149     elif not isNone(box1):
    150         box(box1)
    151     else:
    152         box()
    153 
    154 def compIPPphoto(summaryTable,mode):
     37# Big (BIG) XXX: For plotting, need to make sure that ALL columns have
     38# "good" data in exactly the same rows, otherwise lose correspondence
     39# between them. E.g.: for every column, create a logical vector saying
     40# whether it's 'good' or not, and then plot things only for those rows
     41# where both are good
     42
     43plotcol_tlist = [
     44    ('d_mag','d_sky','scatter'),
     45    ('d_x','d_y','scatter')
     46    ]
     47
     48
     49def compIPPphoto(summaryTable,mode,plotcol_tlist=plotcol_tlist):
    15550    """summaryTable: .fits table for output
    15651    mode: new or append for creating summaryTable new or appending current run's output to it.
    … …  
    17873    rowtuple_list = []
    17974
     75   
    18076    chipfile_l,fpObjc_l = makePlan()
    18177    for chipfile,fpObjc in zip(chipfile_l,fpObjc_l):
    18278        matchtable = matchSdssPs1(fpObjc,chipfile)
    183         res_hash = computeStatistics(matchtable)
    184         # Sort by column names to make sure order is identical in all
    185         # rows of the table, and output is more legible
     79        res_hash, deltas_hash = computeStatistics(matchtable)
     80        # Sort res_hash by its column names to make sure order is
     81        # identical in all rows of the table, and output is more
     82        # legible
    18683        vallist,keylist = valuesKeysSortedByKeys(res_hash)
    18784        rowtuple_list.append(vallist)
     85        plotStatsOnefile(deltas_hash,matchtable,plotcol_tlist)
    18886    newrows = numpy.rec.array(rowtuple_list,names=keylist)
    18987    tabhdu = tabHDUfromRecArray(newrows)
    … …  
    19290    else:
    19391        appendFitsTable(summaryTable,tabhdu)
     92
     93def getOutnameStatsOnefile(matchtable,kind,col1,col2=None,format='eps'):
     94    import re
     95    # Construct output filename
     96    root = re.sub('(\.[sc]mf|\.fits?)$','',matchtable)
     97    if isNone(col2):
     98        outname = '%s_%s_%s.%s' % (root,kind,col1,format)
     99    else:
     100        outname = '%s_%s_%s_%s.%s' % (root,kind,col1,col2,format)
     101    return outname
     102   
     103def plotStatsOnefile(deltas_hash,matchtable,plotcol_tlist,format='eps'):
     104    """Make diagnostic plots for a single table, based on values in deltas_hash"""
     105    for troika in plotcol_tlist:
     106        col1name = troika[0]
     107        col2name = troika[1]
     108        plottype = troika[2]
     109        outname = getOutnameStatsOnefile(matchtable,plottype,col1name,col2name,format=format)
     110        smOpenPlot(outname,format=format)
     111        if plottype == 'scatter':
     112            smScatterPlot(deltas_hash[col1name],deltas_hash[col2name],\
     113                              xlab=col1name,ylab=col2name)
     114        smClosePlot()
    194115
    195116def valuesKeysSortedByKeys(hash):
    … …  
    301222        ,'d_magerr':['psfcountserr','PSF_INST_MAG_SIG']
    302223        }
    303 
     224    colval_hash = {}
    304225    ismag = re.compile('mag')
    305226    iscounts = re.compile('counts')
    … …  
    334255            SDSScol_good,PS1col_good,SDSScounts_good = filterGoodVal3(SDSScol,PS1col,SDSScounts)
    335256            SDSScol_good = 2.5/log(10.)*SDSScol_good/SDSScounts_good
    336 
    337257        delta = SDSScol_good - PS1col_good
     258        colval_hash[outcol] = delta
    338259        avg = delta.mean()
    339260        lowq,med,upq = stats_med(delta)
    … …  
    351272    writeTable(tablename,newprimhdu,newtab)
    352273    infile_handle.close()
    353     return outhash
     274    return outhash,colval_hash
    354275
    355276def writeTable(filename,primhdu,tabhdu):
    … …  
    450371
    451372def matchSdssPs1(SDSSfpObjc,PS1cmf,xoff=0.5,yoff=0.5,matchrad=0.7):
    452     """Call matchByPos to match an SDSS fpObjc.fits and a PS1 bla.cmf file"""
     373    """Call matchByPos to match an SDSS fpObjc.fits and a PS1 bla.cmf
     374    file."""
    453375    import pyfits
     376
     377    def getOutname(SDSSfile,PS1file,sdssbandstr):
     378        import re
     379        # Construct output filename
     380        fitsend = re.compile('(\.[sc]mf|\.fits?)$')
     381        SDSSroot = fitsend.sub('',SDSSfile)
     382        PS1root = fitsend.sub('',PS1file)
     383        # Now chop off any leading path components
     384        pathroot = re.compile('.*/')
     385        SDSSroot = pathroot.sub('',SDSSroot )
     386        PS1root = pathroot.sub('',PS1root)
     387        return "match_%s___%s___%s.fits" % (sdssbandstr,SDSSroot,PS1root)
     388   
    454389    filters = ['u','g','r','i','z']
    455390    # Read primary  header of PS1 file to work out band
    … …  
    480415    ps1pos = "\'x_psf y_psf\'"
    481416
    482     outname = getOutname(SDSSfpObjc,PS1cmf,sdssbandstr)   
     417    outname = getOutname(SDSSfpObjc,PS1cmf,sdssbandstr)
    483418    matchByPos(SDSSfpObjc,PS1cmf,outname,sdsspos,ps1pos,tolerance=matchrad,\
    484419                   duptag1='_sdss',duptag2='_ps1',\
    … …  
    511446    return outhash
    512447
    513 def getOutname(SDSSfile,PS1file,sdssbandstr):
    514     import re
    515     # Construct output filename
    516     fitsend = re.compile('(\.cm[sf]|\.fits?)$')
    517     SDSSroot = fitsend.sub('',SDSSfile)
    518     PS1root = fitsend.sub('',PS1file)
    519     # Now chop off any leading path components
    520     pathroot = re.compile('.*/')
    521     SDSSroot = pathroot.sub('',SDSSroot )
    522     PS1root = pathroot.sub('',PS1root)
    523     return "match_%s___%s___%s.fits" % (sdssbandstr,SDSSroot,PS1root)
    524    
    525    
    526448def matchByPos(in1,in2,out,colnames1,colnames2,tolerance=0.5,duptag1='_sdss',duptag2='_ps1',filter1="",filter2=""):
    527449    """Match two fits tables by cartesian 2-d position using stilts/tmatch2 from topcat"""
    … …  
    542464    print tmatch_cmd
    543465    retval = subprocess.call(tmatch_cmd,shell=True)
     466
     467def smdefaults():
     468    # Will these be persistent across multiple imports and going out
     469    # of scope? Perhaps only if I do a global 'import sm'. Or it
     470    # should work if called from something where sm is in scope?
     471    sm.expand(1.8)
     472    sm.lweight(5)
     473
     474# Plotting convention will be like in sm: you need to open a file,
     475# then you plot to it with one or multiple commands, then you close it
     476# to write the file.
     477
     478def smOpenPlot(filename,format='eps'):
     479    """Issue sm device command for file of given format:
     480    eps -> postfile
     481    else just use given format as 'device'"""
     482    from sm import device, erase
     483    if filename == 'x11':
     484        device('x11')
     485    elif format == 'eps':
     486        device('postfile '+filename)
     487    else:
     488        device(format +' '+filename)   
     489    erase()
     490   
     491def smClosePlot():
     492    from sm import device
     493    device('nodevice')
     494
     495def isNone(something):
     496    return isinstance(something,type(None))
     497
     498def smHistoPlot(vec,step=None,minbin=None,maxbin=None,nbins=None,\
     499                    xlab=None,ylab=None,xrange=None,yrange=None,\
     500                    box1=None,box2=None,box3=None,box4=None,\
     501                    append=False):
     502    """Plot a histogram, intelligently deriving bins from the given
     503    parameters if they are given intelligently.  Otherwise, silently
     504    do nothing."""
     505    import sm
     506    if isNone(minbin):
     507        minbin = min(vec)
     508    if isNone(maxbin):
     509        maxbin = max(vec)
     510    # I am assuming bins are bin centers
     511    if not isNone(step):
     512        nbins = (maxbin-minbin)/step+1
     513    if isNone(nbins):
     514        return
     515    histo,leftbinedges = histogram(vec,nbins,[minbin,maxbin])
     516    bincenters = leftbinedges + 0.5*(leftbinedges[1]-leftbinedges[0])
     517
     518    if not append:
     519        smSetup(bincenters,histo,xrange,yrange,xlab,ylab,box1,box2,box3,box4)
     520    sm.histogram(bincenters,histo)
     521
     522def smLinePlot(x,y,ltype=0,xlab=None,ylab=None,xrange=None,yrange=None,\
     523                    box1=None,box2=None,box3=None,box4=None,\
     524                    append=False):
     525    """Make an sm scatter plot on current device. If append=True,
     526    overplot with current limits. Otherwise, draw box box1 box2 box3 box4"""
     527    import sm
     528    sm.expand(1.8)
     529    sm.lweight(5)
     530    # Silently ignore any problems with plot generation
     531    try:
     532        if not append:
     533            smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4)
     534        sm.ltype(ltype)
     535        sm.connect(x,y)
     536    except:
     537        pass
     538
     539
     540def smScatterPlot(x,y,ptype=41,xlab=None,ylab=None,xrange=None,yrange=None,\
     541                    box1=None,box2=None,box3=None,box4=None,\
     542                    append=False):
     543    """Make an sm scatter plot on current device. If append=True,
     544    overplot with current limits. Otherwise, draw box box1 box2 box3 box4"""
     545    import sm
     546    sm.expand(1.8)
     547    sm.lweight(5)
     548    # Silently ignore any problems with plot generation
     549    try:
     550        if not append:
     551            smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4)
     552        sm.ptype(ptype)
     553        sm.points(x,y)
     554    except:
     555        pass
     556
     557def smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4):
     558    import sm
     559    sm.erase()
     560    if isNone(xrange):
     561        xrange = x
     562    if isNone(yrange):
     563        yrange = y
     564    sm.limits(x,y)
     565    smBox(box1,box2,box3,box4)
     566    if not isNone(xlab):
     567        sm.xlabel(xlab)
     568    if not isNone(ylab):
     569        sm.ylabel(ylab)
     570
     571def smBox(box1,box2,box3,box4):
     572    from sm import box
     573    if not isNone(box4):
     574        box(box1,box2,box3,box4)
     575    elif not isNone(box3):
     576        box(box1,box2,box3)
     577    elif not isNone(box2):
     578        box(box1,box2)
     579    elif not isNone(box1):
     580        box(box1)
     581    else:
     582        box()
     583
     584   
Note: See TracChangeset for help on using the changeset viewer.