Index: /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py =================================================================== --- /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 19899) +++ /branches/sj_ippTests_branch_20080929/ippTests/compIPPphoto.py (revision 19900) @@ -35,15 +35,11 @@ # + 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,12 +72,12 @@ 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,4 +86,5 @@ else: appendFitsTable(summaryTable,tabhdu) + return column_hash,goodrow_hash def getOutnameStatsOnefile(matchtable,kind,col1,col2=None,format='eps'): @@ -101,6 +98,7 @@ 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,8 +106,21 @@ 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,9 +143,9 @@ 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,11 +220,28 @@ 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,39 +250,56 @@ ,'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,5 +308,5 @@ 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,7 +315,7 @@ 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,15 +326,26 @@ 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,5 +482,5 @@ f.close() - return outname + return outname,sdssbandstr,bandindex def updateHeaderFromHash(header,hash): @@ -520,8 +576,8 @@ 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,10 +589,10 @@ 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,5 +607,5 @@ 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