IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 19801


Ignore:
Timestamp:
Oct 1, 2008, 11:58:48 PM (18 years ago)
Author:
Sebastian Jester
Message:

move things into subroutines for better legibility

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py

    r19800 r19801  
    1515# (which in turn has dependencies)
    1616#
    17 #
     17# Also see: http://kiawe.ifa.hawaii.edu/IPPwiki/index.php/SimtestVerification
     18#
     19# Todo:
     20# * rms etc. with outlier rejection - e.g. by histogram core fit.
     21# * Plots
     22#   - Per-field:
     23#     + d_x vs. d_y
     24#     + magnitude histograms
     25#     + flux vs. sky
     26#     + magdiff vs. mag
     27#     + magdiff vs. sky
     28#     + magdiff vs. skydiff
     29#     + something about PSF sizes; use Michigan moments for photo -
     30#       M_rr_cc_psf is the size of the reconstructed PSF, and m_rr_cc is the size of the object
     31#     + magdiff, mag vs. d_x^2+d_y^2
     32#     + histogram of chi - but meaningful only for repeat measurements?
     33#   - per-run:
     34#     + histograms of everything - recovery fractions, means, medians, quartiles
     35#     + Trends with field number, seeing etc.
    1836
    1937def compIPPphoto(summaryTable,mode):
    … …  
    3957    """
    4058    import pyfits,numpy
    41    
    42 
    43     chipfile_l,fpObjc_l = makePlan()
    44 
    4559    if mode.lower() not in ["new","append"]:
    4660        raise ValueError("Must specify mode=new or mode=append")
    4761    rowtuple_list = []
     62
     63    chipfile_l,fpObjc_l = makePlan()
    4864    for chipfile,fpObjc in zip(chipfile_l,fpObjc_l):
    4965        matchtable = matchSdssPs1(fpObjc,chipfile)
    5066        res_hash = computeStatistics(matchtable)
    51         # Sort by column names to make sure order is correct, and
    52         # output is more legible
    53         vallist = []
    54         keylist = res_hash.keys()
    55         keylist.sort()
    56         for field in keylist:
    57             vallist.append(res_hash[field])
     67        # Sort by column names to make sure order is identical in all
     68        # rows of the table, and output is more legible
     69        vallist,keylist = valuesKeysSortedByKeys(res_hash)
    5870        rowtuple_list.append(vallist)
    5971    newrows = numpy.rec.array(rowtuple_list,names=keylist)
    … …  
    6274        tabhdu.writeto(summaryTable,clobber=True)
    6375    else:
    64         oldtabhdulist = pyfits.open(summaryTable)
    65         nrows_old = len(oldtabhdulist[1].data)
    66         newtabhdu = pyfits.new_table(oldtabhdulist[1].columns,nrows = nrows_old+len(tabhdu.data))
    67         for col in newtabhdu.columns.names:
    68             newtabhdu.data.field(col)[nrows_old:] = tabhdu.data.field(col)
    69         oldtabhdulist.close()
    70         newtabhdu.writeto(summaryTable,clobber=True)
    71 
    72 
     76        appendFitsTable(summaryTable,tabhdu)
     77
     78def valuesKeysSortedByKeys(hash):
     79    """Return a list of hash's values and keys that is ordered by the names of the keys"""
     80    vallist = []
     81    keylist = hash.keys()
     82    keylist.sort()
     83    for field in keylist:
     84        vallist.append(hash[field])
     85    return vallist,keylist
     86       
     87def appendFitsTable(fitsfile,tabHDU,HDU=1):
     88    """Append rows in tabHDU to existing fits table file"""
     89    oldtabhdulist = pyfits.open(fitsfile)
     90    nrows_old = len(oldtabhdulist[HDU].data)
     91    newtabhdu = pyfits.new_table(oldtabhdulist[HDU].columns,nrows = nrows_old+len(tabHDU.data))
     92    for col in newtabhdu.columns.names:
     93        newtabhdu.data.field(col)[nrows_old:] = tabHDU.data.field(col)
     94    oldtabhdulist.close()
     95    return newtabhdu.writeto(fitsfile,clobber=True)
     96   
     97       
    7398def tabHDUfromRecArray(recarr):
    7499    """Generate a table HDU from a record array"""
    … …  
    84109            # This is potentially dangerous, but for now, the only
    85110            # string column I am writing is the filter, i.e. 1A should
    86             # be enough already.
     111            # be enough already. Perhaps use the length of the entry
     112            # as column format?
    87113            format = '1A'
    88114        else:
    … …  
    120146    return lowerquartile,median,upperquartile
    121147
    122 def computeStatistics(tablename):
     148def computeStatistics(tablename,copyfields_list = ['RUN','RERUN','CAMCOL','FIELD','FILTER','FWHM_X','FWHM_Y']):
    123149    """Compute desired statistics for offsets etc. and add them as
    124150    ESO-style HIERARCH fields to the table header (8 characters are
    … …  
    136162    # For this to work, table needs to be called match_r_... for r-band
    137163   
    138     f = pyfits.open(tablename,mode='update')
    139     table = f[1]
     164    infile_handle = pyfits.open(tablename,mode='update')
     165    table = infile_handle[1]
    140166    h = table.header
    141167    table_data = table.data
    142168    # Header names that are going to be copied to output table as
    143169    # "reference columns"
    144     copyfields_list = ['RUN','RERUN','CAMCOL','FIELD','FILTER']
    145170    outhash = {}
    146171    for f in copyfields_list:
    147         outhash[f] = h[f]   
     172        outhash[f] = h[f]
     173
     174    filtname = outhash['FILTER']
    148175
    149176    # These are just the columns; need to get a slice with the correct array index later
    … …  
    167194        # necessary for a given column (e.g. for objc_colc it won't
    168195        # be...)
    169         SDSScol = table_data.field(SDSScolname)[:,filterID[filtnam]]
     196        SDSScol = table_data.field(SDSScolname)[:,filterID[filtname]]
    170197        PS1col = table_data.field(PS1colname)
    171198        SDSScol_good,PS1col_good = filterGoodVal(SDSScol,PS1col)
    … …  
    186213            outcoll.append(instmagerr_col)
    187214            # Compute array for internal use
    188             SDSScounts = table_data.field('psfcounts')[:,filterID[filtnam]]
     215            SDSScounts = table_data.field('psfcounts')[:,filterID[filtname]]
    189216            SDSScol_good,PS1col_good,SDSScounts_good = filterGoodVal3(SDSScol,PS1col,SDSScounts)
    190217            SDSScol_good = 2.5/log(10.)*SDSScol_good/SDSScounts_good
    … …  
    200227        for label,value in zip(outtablabel,[avg,rms,med,lowq,upq]):
    201228            outhash[label]=value
    202     f.close()
    203229   
    204230    newtab = pyfits.new_table(table.columns+pyfits.ColDefs(outcoll),header=h)
    205     newprimhdu = pyfits.PrimaryHDU(header=f[0].header)
     231    newprimhdu = pyfits.PrimaryHDU(header=infile_handle[0].header)
    206232    writeTable(tablename,newprimhdu,newtab)
     233    infile_handle.close()
    207234    return outhash
    208235
    … …  
    262289    from run, camcol, field read from .cmf primary header, which
    263290    replicates most of the fpC primary header."""
    264    
     291     
     292    # For testing:
     293    return ['1056-0192.421/1056-0192.421.ch.421.CHIP1.cmf'],['1056-0192.421/fpObjc-001056-1-0192.fit']
     294
    265295    from glob import glob
    266296    fpObjcDir = '/IPP/data/SDSS/stripe82/coadd/input/fpObjcs/'
    … …  
    282312    Nfields = 18
    283313    camcollist = range(1,6+1)
    284     # For testing:
    285     # pairlist = ['1056-0192.421/1056-0192.421.ch.421.CHIP1.cmf'],['1056-0192.421/fpObjc-001056-1-0192.fit']
    286314    sdsslist = []
    287315    ps1list = []
    … …  
    304332def matchSdssPs1(SDSSfpObjc,PS1cmf,xoff=0.5,yoff=0.5,matchrad=0.7):
    305333    """Call matchByPos to match an SDSS fpObjc.fits and a PS1 bla.cmf file"""
    306     import pyfits,re
     334    import pyfits
    307335    filters = ['u','g','r','i','z']
    308336    # Read primary  header of PS1 file to work out band
    309     ps1f = pyfits.open(PS1cmf)
    310     h=ps1f[0].header
    311     sdssbandstr = h['FILTER']
    312     ps1f.close()
    313 
    314     sdssf = pyfits.open(SDSSfpObjc)
    315     h=sdssf[0].header
     337    ps1copyfields_hash = headerfieldHash(['FWHM_X','FWHM_Y','FILTER'],PS1cmf,HDU=0)
     338    sdssbandstr = ps1copyfields_hash['FILTER']
    316339    bandindex = filters.index(sdssbandstr)
    317     run = h['RUN']
    318     rerun = h['RERUN']
    319     camcol=h['CAMCOL']
    320     field = h['FIELD']
    321     sdssf.close()
     340    sdsscopyfields_hash = headerfieldHash(['RUN','RERUN','CAMCOL','FIELD'],SDSSfpObjc,0)
     341
    322342    # Note position offset - y_sdss-y_psf1 ~ 0.6, x_sdss-x_ps1 ~ 0.35 (in 1056-r1-0192)
    323343    # width of distribution is ~ 0.4
    … …  
    340360    sdsspos = '\'colc_%s rowc_%s\''% (sdssbandstr,sdssbandstr)
    341361    ps1pos = "\'x_psf y_psf\'"
     362
     363    outname = getOutname(SDSSfpObjc,PS1cmf,sdssbandstr)   
     364    matchByPos(SDSSfpObjc,PS1cmf,outname,sdsspos,ps1pos,tolerance=matchrad,\
     365                   duptag1='_sdss',duptag2='_ps1',\
     366                   filter1='\''+sdssfilt+'\'',filter2='\''+ps1filt+'\'')
     367
     368    f=pyfits.open(outname,mode='update')
     369    h=f[1].header
     370    updateHeaderFromHash(h,sdsscopyfields_hash)
     371    updateHeaderFromHash(h,ps1copyfields_hash)
     372    f.close()
     373   
     374    return outname
     375
     376def updateHeaderFromHash(header,hash):
     377    """Add fields in hash to header"""
     378    for k in hash.keys():
     379        header.update(k,hash[k])
     380
     381
     382def headerfieldHash(fieldlist,file,HDU=0):
     383    """Return a hash containing listed fits headers from file's HDU
     384    (this is like a hash slice, which doesn't seem to be possible in
     385    python?)"""
     386    import pyfits
     387    f=pyfits.open(file)
     388    h=f[HDU].header
     389    outhash = {}
     390    for field in fieldlist:
     391        outhash[field] = h[field]
     392    return outhash
     393
     394def getOutname(SDSSfile,PS1file,sdssbandstr):
     395    import re
    342396    # Construct output filename
    343397    fitsend = re.compile('(\.cm[sf]|\.fits?)$')
    344     SDSSroot = fitsend.sub('',SDSSfpObjc)
    345     PS1root = fitsend.sub('',PS1cmf)
     398    SDSSroot = fitsend.sub('',SDSSfile)
     399    PS1root = fitsend.sub('',PS1file)
    346400    # Now chop off any leading path components
    347401    pathroot = re.compile('.*/')
    348402    SDSSroot = pathroot.sub('',SDSSroot )
    349403    PS1root = pathroot.sub('',PS1root)
    350     outname = "match_%s___%s___%s.fits" % (sdssbandstr,SDSSroot,PS1root)
    351     matchByPos(SDSSfpObjc,PS1cmf,outname,sdsspos,ps1pos,tolerance=matchrad,\
    352                    duptag1='_sdss',duptag2='_ps1',\
    353                    filter1='\''+sdssfilt+'\'',filter2='\''+ps1filt+'\'')
    354     f=pyfits.open(outname,mode='update')
    355     h=f[1].header
    356     h.update('FILTER',sdssbandstr)
    357     h.update('RUN',run)
    358     h.update('RERUN',rerun)
    359     h.update('CAMCOL',camcol)
    360     h.update('FIELD',field)
    361     f.close()
    362    
    363     return outname
     404    return "match_%s___%s___%s.fits" % (sdssbandstr,SDSSroot,PS1root)
     405   
    364406   
    365407def matchByPos(in1,in2,out,colnames1,colnames2,tolerance=0.5,duptag1='_sdss',duptag2='_ps1',filter1="",filter2=""):
Note: See TracChangeset for help on using the changeset viewer.