Changeset 19900

Timestamp:

Oct 5, 2008, 11:37:51 PM (18 years ago)

Author:

Sebastian Jester

Message:

Use logical vectors as index to numpy arrays for selecting goodrows; allow logical vectors in plots for plotting only 'good' values; pass around all columns, computed differences and their associated logical vectors for plotting

File:

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

branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py

@@ -35 +35 @@
 #     + Trends with field number, seeing etc.
 
-# 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
+# XXX move computation of new columns into matching script
 
 plotcol_tlist = [
-    ('d_mag','d_sky','scatter'),
-    ('d_x','d_y','scatter')
+    ('d_mag','d_sky','scatter')
+    ,('d_x','d_y','scatter')
+    ,('sky_ps1','d_mag','scatter')
     ]
-
 
 def compIPPphoto(summaryTable,mode,plotcol_tlist=plotcol_tlist):
@@ -76 +72 @@
     chipfile_l,fpObjc_l = makePlan()
     for chipfile,fpObjc in zip(chipfile_l,fpObjc_l):
-        matchtable = matchSdssPs1(fpObjc,chipfile)
-        res_hash, deltas_hash = computeStatistics(matchtable)
+        matchtable,filter_name,bandindex = matchSdssPs1(fpObjc,chipfile)
+        stats_hash, column_hash, goodrow_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)
+        vallist,keylist = valuesKeysSortedByKeys(stats_hash)
         rowtuple_list.append(vallist)
-        plotStatsOnefile(deltas_hash,matchtable,plotcol_tlist)
+        plotStatsOnefile(column_hash,goodrow_hash,matchtable,plotcol_tlist,bandindex)
     newrows = numpy.rec.array(rowtuple_list,names=keylist)
     tabhdu = tabHDUfromRecArray(newrows)
@@ -90 +86 @@
     else:
         appendFitsTable(summaryTable,tabhdu)
+    return column_hash,goodrow_hash
 
 def getOutnameStatsOnefile(matchtable,kind,col1,col2=None,format='eps'):
@@ -101 +98 @@
     return outname
     
-def plotStatsOnefile(deltas_hash,matchtable,plotcol_tlist,format='eps'):
-    """Make diagnostic plots for a single table, based on values in deltas_hash"""
+def plotStatsOnefile(values_hash,goodrow_hash,matchtable,plotcol_tlist,bandindex,format='eps'):
+    """Make diagnostic plots for a single table, based on values
+    in values_hash"""
     for troika in plotcol_tlist:
         col1name = troika[0]
@@ -108 +106 @@
         plottype = troika[2]
         outname = getOutnameStatsOnefile(matchtable,plottype,col1name,col2name,format=format)
+        values1 = values_hash[col1name]
+        goodrows1 = goodrow_hash[col1name]
+        values2 = values_hash[col2name]
+        goodrows2 = goodrow_hash[col2name]
+        # Slice out depth for current filter if necessary
+        if len(values1.shape) > 1:
+            values1 = values1[:,bandindex]
+            goodrows1 = goodrows1[:,bandindex]
+        if len(values2.shape) > 1:
+            values2 = values2[:,bandindex]
+            goodrows2 = goodrows2[:,bandindex]
+        goodrows = goodrows1 & goodrows2
+        # print col1name, col2name, sum(goodrows), sum(values1 > 1e3)
+        # print values1[goodrows & (values1 > 1e3)]
         smOpenPlot(outname,format=format)
         if plottype == 'scatter':
-            smScatterPlot(deltas_hash[col1name],deltas_hash[col2name],\
-                              xlab=col1name,ylab=col2name)
+            smScatterPlot(values1,values2,logical=goodrows,xlab=col1name,ylab=col2name)
         smClosePlot()
 
@@ -132 +143 @@
     oldtabhdulist.close()
     return newtabhdu.writeto(fitsfile,clobber=True)
-    
+
         
 def tabHDUfromRecArray(recarr):
     """Generate a table HDU from a record array"""
-    # create column names from recarray._names
+    # create column names from recarray.dtype.names (NumPy style)
     # Assume float as data type
     # copy data field-wise to table HDU.
@@ -209 +220 @@
     for f in copyfields_list:
         outhash[f] = h[f]
-
     filtname = outhash['FILTER']
+
+    goodrow_hash = {}
+    colval_hash = {}
+    # Find out which entries are "good" in every column
+    for column in table.columns.names:
+        coldata = table_data.field(column)
+        # Check column dimensions:
+        # If 3D, it's one of the profile columns, and we can't plot those anyway, so just skip them.
+        if len(coldata.shape) > 2:
+            continue
+        # If 2D, it's an array, and we need to slice out the right
+        # depth
+        elif len(coldata.shape) == 2:
+            goodrow_hash[column.lower()] = goodValBool(table_data.field(column))[:,filterID[filtname]]
+            colval_hash[column.lower()] = array(table_data.field(column))[:,filterID[filtname]]
+        else:
+            goodrow_hash[column.lower()] = goodValBool(table_data.field(column))
+            colval_hash[column.lower()] = array(table_data.field(column))
 
     # These are just the columns; need to get a slice with the correct array index later
     colname_hash = {
-        'd_x':['colc','x_psf']
-        ,'d_y':['rowc','y_psf']
+        'd_x':['colc','X_PSF']
+        ,'d_y':['rowc','Y_PSF']
         ,'d_sky':['sky_sdss','SKY_ps1']
         ,'d_skyerr':['skyErr','SKY_SIGMA']
@@ -222 +250 @@
         ,'d_magerr':['psfcountserr','PSF_INST_MAG_SIG']
         }
-    colval_hash = {}
     ismag = re.compile('mag')
     iscounts = re.compile('counts')
     outcoll = []
+
+    
     for outcol,list in colname_hash.iteritems():
         SDSScolname,PS1colname=list
-        # Slice out the proper filter - may need to add another
-        # boolean dictionary that says whether or not this is
-        # necessary for a given column (e.g. for objc_colc it won't
-        # be...)
-        SDSScol = table_data.field(SDSScolname)[:,filterID[filtname]]
-        PS1col = table_data.field(PS1colname)
-        SDSScol_good,PS1col_good = filterGoodVal(SDSScol,PS1col)
+        # Slice out the proper filter if necessary. This and other
+        # array() calls are to make sure that we are working with
+        # numpy.ndarrays, not with numarray's arrays, which pyfits
+        # returns.
+        SDSScol = array(table_data.field(SDSScolname))
+        if len(SDSScol.shape) == 2:
+            SDSScol = SDSScol[:,filterID[filtname]]
+        PS1col = array(table_data.field(PS1colname))
+        bothgood_bool = goodrow_hash[SDSScolname.lower()] & goodrow_hash[PS1colname.lower()]
+        SDSScol_good = SDSScol[bothgood_bool]
+        PS1col_good = PS1col[bothgood_bool]
+        goodrow_hash[outcol] = bothgood_bool
         # Compute SDSS instrumental magnitude if necessary
-        if PS1colname == 'PSF_INST_MAG' and SDSScolname.lower() == 'psfcounts':
+        if PS1colname.lower() == 'psf_inst_mag' and SDSScolname.lower() == 'psfcounts':
             # Create column with SDSS instrumental magnitude that can
             # be written to table for diagnostic use
-            instmag_sdss_arr = -2.5*log10(table_data.field(SDSScolname))
+            sdsspsfinstmag_colname = 'psfinstmag_sdss'
+            instmag_sdss_arr = -2.5*log10(array(table_data.field(SDSScolname)))
             instmag_col = pyfits.Column(name='psfinstmag_sdss',format='5E',array=instmag_sdss_arr)
             outcoll.append(instmag_col)
+            goodrow_hash[sdsspsfinstmag_colname.lower()] = goodrow_hash[SDSScolname.lower()]
             # Compute array for internal use
             SDSScol_good = -2.5*log10(SDSScol_good)
 
-        if PS1colname == 'PSF_INST_MAG_SIG' and SDSScolname.lower() == 'psfcountserr':
+        if PS1colname.lower() == 'psf_inst_mag_sig' and SDSScolname.lower() == 'psfcountserr':
             # Create column with SDSS instrumental magnitude error
             # that can be written to table for diagnostic use
-            instmagerr_sdss_arr = 2.5/log(10.)*table_data.field('psfcountserr')/table_data.field('psfcounts')
-            instmagerr_col = pyfits.Column(name='psfinstmagerr_sdss',format='5E',array=instmagerr_sdss_arr)
+            sdsspsfinstmagerr_colname = 'psfinstmagerr_sdss'
+            instmagerr_sdss_arr = array(2.5/log(10.)*table_data.field('psfcountserr')/table_data.field('psfcounts'))
+            instmagerr_col = pyfits.Column(name=sdsspsfinstmagerr_colname,format='5E',array=instmagerr_sdss_arr)
             outcoll.append(instmagerr_col)
+            goodrow_hash[sdsspsfinstmagerr_colname.lower()] = goodrow_hash['psfcountserr'] & \
+                goodrow_hash['psfcounts']
             # Compute array for internal use
-            SDSScounts = table_data.field('psfcounts')[:,filterID[filtname]]
-            SDSScol_good,PS1col_good,SDSScounts_good = filterGoodVal3(SDSScol,PS1col,SDSScounts)
+            SDSScounts = array(table_data.field('psfcounts')[:,filterID[filtname]])
+            all3good_bool = bothgood_bool & goodrow_hash[sdsspsfinstmagerr_colname.lower()]
+            SDSScol_good = SDSScol[all3good_bool]
+            PS1col_good = PS1col[all3good_bool]
+            SDSScounts_good = SDSScounts[all3good_bool]
             SDSScol_good = 2.5/log(10.)*SDSScol_good/SDSScounts_good
+            goodrow_hash[outcol] = all3good_bool
         delta = SDSScol_good - PS1col_good
-        colval_hash[outcol] = delta
+        # Store *all* values in return hash for later plotting; good
+        # values will be filtered out during plotting
+        colval_hash[outcol] = SDSScol - PS1col
         avg = delta.mean()
         lowq,med,upq = stats_med(delta)
@@ -263 +308 @@
         outheaderlabel = [operator.add(outcol+' ',itm).upper() for itm in label_l]
         outtablabel = [operator.add(outcol+'_',itm).upper() for itm in label_l]
-        saveHierarchHeaderList(h,outheaderlabel,[avg,rms,med,lowq,upq])
+        saveHierarchHeaderList(h,outheaderlabel,[float(avg),float(rms),float(med),float(lowq),float(upq)])
         # and in return variable
         for label,value in zip(outtablabel,[avg,rms,med,lowq,upq]):
@@ -270 +315 @@
     newtab = pyfits.new_table(table.columns+pyfits.ColDefs(outcoll),header=h)
     newprimhdu = pyfits.PrimaryHDU(header=infile_handle[0].header)
+    infile_handle.close()
     writeTable(tablename,newprimhdu,newtab)
-    infile_handle.close()
-    return outhash,colval_hash
+    return outhash,colval_hash,goodrow_hash
 
 def writeTable(filename,primhdu,tabhdu):
@@ -281 +326 @@
     hdulst.writeto(filename,clobber=True)
 
-def filterGoodVal(col1,col2):
+def goodValBool(col):
+    """Return a bool array that contains True where
+    col1 doesn't have  NaNs nor SDSS "no-value" and
+    "no-error" flags -9999 and -1000."""
+    from numpy import isfinite,array
+    goodcondition = isfinite(col) 
+    for val in [-9999,-1000]:
+        goodcondition &= (col != val)
+    return goodcondition
+
+def filterGoodVal2(col1,col2):
     """Return a pair of arrays that contains those items where both
     col1 and col2 have  NaNs and SDSS "no-value" and
     "no-error" flags -9999 and -1000 filtered out."""
     from numpy import isfinite,logical_and,array
-    # I don't understand why this doesn't work...
-    # 
-    #goodcondition = logical_and(isfinite(col1),isfinite(col2))
-    #col1 = col1[goodcondition]
-    #col2_good = col2[goodcondition]
-    #return col1,col2
+    goodcondition = isfinite(col1) & isfinite(col2)
+    for val in [-9999,-1000]:
+        goodcondition &= (col1 != val)
+        goodcondition &= (col2 != val)
+    col1 = col1[goodcondition]
+    col2_good = col2[goodcondition]
+    return col1,col2
     #
     #... so I'm doing this:
@@ -426 +482 @@
     f.close()
     
-    return outname
+    return outname,sdssbandstr,bandindex
 
 def updateHeaderFromHash(header,hash):
@@ -520 +576 @@
     sm.histogram(bincenters,histo)
 
-def smLinePlot(x,y,ltype=0,xlab=None,ylab=None,xrange=None,yrange=None,\
+def smLinePlot(x,y,logical=None,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,
+    """Make an sm line plot on current device. If append=True,
     overplot with current limits. Otherwise, draw box box1 box2 box3 box4"""
     import sm
@@ -533 +589 @@
             smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4)
         sm.ltype(ltype)
-        sm.connect(x,y)
+        sm.connect(x,y,logical)
     except:
         pass
 
 
-def smScatterPlot(x,y,ptype=41,xlab=None,ylab=None,xrange=None,yrange=None,\
+def smScatterPlot(x,y,logical=None,ptype=41,xlab=None,ylab=None,xrange=None,yrange=None,\
                     box1=None,box2=None,box3=None,box4=None,\
                     append=False):
@@ -551 +607 @@
             smSetup(x,y,xrange,yrange,xlab,ylab,box1,box2,box3,box4)
         sm.ptype(ptype)
-        sm.points(x,y)
+        sm.points(x,y,logical)
     except:
         pass