# Routines for comparing IPP outputs against photometry from other # sources. # # Currently only SDSS's photo, but meant to extend to others including # simtest # # Author: Sebastian Jester jester at mpia.de # # Depends on: # # topcat with 'stilts' shell script from # http://www.star.bris.ac.uk/~mbt/topcat/ and # # pyfits from http://www.stsci.edu/resources/software_hardware/pyfits # (which in turn has dependencies) # # Also see: http://kiawe.ifa.hawaii.edu/IPPwiki/index.php/SimtestVerification # # Todo: # * rms etc. with outlier rejection - e.g. by histogram core fit. # * Plots # - Per-field: # + d_x vs. d_y # + magnitude histograms # + flux vs. sky # + magdiff vs. mag # + magdiff vs. sky # + magdiff vs. skydiff # + something about PSF sizes; use Michigan moments for photo - # M_rr_cc_psf is the size of the reconstructed PSF, and m_rr_cc is the size of the object # + magdiff, mag vs. d_x^2+d_y^2 # + histogram of chi - but meaningful only for repeat measurements? # + object density - store! # # - per-run: # + histograms of everything - recovery fractions, means, medians, quartiles # + Trends with field number, seeing etc. # + magdiff vs. stellar density # + magdiff vs. presence of bright stars (e.g. number of bright # stars above some typical S/N or counts threshold - 50% # saturation?) # # XXX Todo: # # - Merge identical plots into single .eps file # # - Make html page with jpeg/gif versions of plots # # - Need mechanism to specify default plotting ranges, since outliers # mess up axis ranges. Different quantities need different # ranges. Comparing different plots is also easier if the plot # ranges are the same. # # - towards generality, want to read these things and colname_hash and # label_l from config files # # - perhaps via the RO package and RO.DS9, add diagnostic features via # its ability to talk to ds9. E.g. load the fpC image in ds9, # overlay SDSS objects in one colour, psphot in another, have # different symbols for point vs extended source, etc. Maybe also # use the CDS Aladin API (if one is available) to be able to click # on object in image to pull up photometry. In fact, stilts should # be able to talk to the aladin tool via PLASTIC. See "Starlink User # Note 253" - somewhere on the web as sun253.html # # # # XXX Todo: labels! group by runs? different colours for different # bands; include seeing? plot things against seeing? # # Also: the plotcol lists need to be grouped into things that are to # appear on the same page. For the beginning, that's the entire list # (and that's probably quite enough anyway). # # Looking at the extended flag in PS1 vs. SDSS is done via d_pointsource. # # For running in earnest: # os.chdir('/disk1/jester/IPP/data/SDSS/stripe82/coadd/compare') # execfile('/disk1/jester/usrdevel/ipp_sj/ippTests/compIPPphoto.py') # h1,h2=compIPPphoto('test.fits','new',workdir='/IPP/data/SDSS/stripe82/coadd/compare',runlist=[1056,1755]) # # For testing on laptop: # os.chdir('/Users/jester/science/sdss/coadd/Hawaii/compare_tsObj') # execfile('/Users/jester/usrdevel/ipp_ippTests_branch/ippTests/compIPPphoto.py') # h1,h2=compIPPphoto('test.fits','new',workdir='/Users/jester/science/sdss/coadd/Hawaii/compare_tsObj',runlist=[1056,1755]) # # Here follow two specially formatted lists governing which plots are # made. The first is for plots being done on a field-by-field basis, # giving [x vector(s)],[y vector(s)], 'scatter' plot or # 'histogram. The second is for "summary" plots that plot aggregate # statistics for the entire list of fields supplied, e.g. the number # of PS1 sources in a field vs. the number of SDSS sources in a field. # # Some columns come from copyfields_list; which columns get d_ # difference columns is set by colname_hash ( see help text for # compIPPphoto() ) # # The convention is: for scatter, plot second against first; for # histogram, plot both histograms plotcol_1frame_tlist = [ (['x_psf','objc_colc'],['y_psf','objc_rowc'],'scatter') ,(['psfinstmag_sdss'],['d_mag'],'scatter') ,(['psfinstmag_sdss'],['d_sky'],'scatter') ,(['d_x'],['d_y'],'scatter') ,(['sky_ps1'],['d_mag'],'scatter') ,(['d_sky'],['d_mag'],'scatter') ,(['d_x'],['d_y'],'histogram') ,(['m_rr_cc_psf'],['d_mag'],'scatter') ,(['sky_sdss'],['sky_ps1'],'scatter') ,(['d_pointsource'],['pointsource_ps1'],'histogram') ,(['sky_sdss'],['sky_ps1'],'histogram') ,(['psfinstmag_sdss'],['psf_inst_mag'],'histogram') ] plotcol_summary_tlist = [ (['d_mag_median'],['d_sky_median'],'scatter') ,(['d_x_median'],['d_y_median'],'scatter') ,(['d_x_median'],['d_y_median'],'histogram') ,(['d_mag_median'],['N_SDSS'],'scatter') ,(['N_PS1'],['N_SDSS'],'scatter') ,(['field','field'],['d_mag_median','d_mag_mean'],'scatter') ] def mergePlots(summaryTable,plotcol_1frame_tlist=plotcol_1frame_tlist,filtlist=['u','g','r','i','z'],\ workdir='/IPP/data/SDSS/stripe82/coadd/compare',gzip=False): from subprocess import call from glob import glob import re import os os.chdir(workdir) outroot = re.sub('(\.[sc]mf|\.fits?)$','',summaryTable) for filt in filtlist: rootstr = "filestats_match_%s"%(filt) globstr = "%s_*.eps"%(rootstr) filelist = glob(globstr) outname = "merged_filestats_%s_%s.eps" % ( filt,outroot ) cmdstr = "gs -dNOPAUSE -sDEVICE=pswrite -sOutputFile=%s %s -c quit" %(outname,' '.join(filelist)) print "Making %s" %(outname) call(cmdstr,shell=True) if gzip: print "gzipping %s" %(outname) call("gzip -f %s" %(outname),shell=True) print "Deleting individual .eps files" for file in filelist: os.remove(file) def compIPPphoto(summaryTable,mode,plotcol_1frame_tlist=plotcol_1frame_tlist,\ plotcol_summary_tlist=plotcol_summary_tlist,\ cmfDir = '/IPP/data/SDSS/stripe82/coadd/prod/run_ipp_20080815/',\ workdir = '/Users/jester/science/sdss/coadd/Hawaii/compare_tsObj', runlist=[1056,1755,3388,3434,3465,4145,4192,4203,4247,4263,5052],\ copyfields_list = ['RUN','RERUN','CAMCOL','FIELD','FILTER','FWHM_MAJ','FWHM_MIN'],\ skipMatch=True): """Match a .cmf file against the corresponding fpObjc file, compute statistics (e.g. median difference of some pair of columns), plot these statistics field-by-field or aggregate basis. Which columns are subtracted from each other is hardcoded into colname_hash in computeStatistics(); the todo list for this source file includes passing colname_hash around as parameter and allowing it to be read from a config file. summaryTable: .fits table for output mode: new or append for creating summaryTable new or appending current run's output to it. plotcol_1frame_tlist: list specifying which plots should be made frame-by-frame (see top of source file) plotcol_summary_tlist: which plots should be made for aggregate quantities for the entire list of fields cmfDir: where psphot output is found; location of SDSS files is currently hardcoded in makePlan workdir: where results are stored (matched-up fits tables, plots) runlist: which SDSS runs to look at copyfields_list: which header fields to copy from the input tables (via the match table) to the output fields, so that we can plot them in the aggregate output. CAUTION: the corresponding fields need to be copied from the input cmf / fpObjc files in matchSdssPs1() - search for ps1copyfields_hash in the source - should probably split this into PS1copyfields and SDSScopyfields and then pass these down to matchSdssPs1() skipMatch: if True, assume cmf / fpObjc matching has already been done, don't repeat it, and only produce statistics and plots Steps: 1. read "plan file" with list of files to be looked at, or otherwise generate that list, e.g. by globbing (paired list fpObjc <-> cmf/smf) for every file: 2. match the files 3. Compute diagnostic columns 4. Compute statistics on diagnostic columns (begin with mean, rms, quartiles) 5. Append summary statistics to master table 6. Generate whatever plots are desired from indidvidual match tables or master table """ import pyfits,numpy if mode.lower() not in ["new","append"]: raise ValueError("Must specify mode=new or mode=append") rowtuple_list = [] column_hash={} stats_hash={} chipfile_l,fpObjc_l = makePlan(cmfDir=cmfDir,runlist=runlist) filters = [] bandindex_hash = {} for chipfile,fpObjc in zip(chipfile_l,fpObjc_l): matchtable,filter_name,bandindex = matchSdssPs1(fpObjc,chipfile,skipMatch=skipMatch) if filter_name not in filters: filters.append(filter_name) bandindex_hash[filter_name] = bandindex stats_hash, column_hash, goodrow_hash = computeStatistics(matchtable,copyfields_list = copyfields_list) # 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(stats_hash) rowtuple_list.append(vallist) plotStatsOnefile(column_hash,goodrow_hash,matchtable,plotcol_1frame_tlist,bandindex,bandid=filter_name) # return column_hash,goodrow_hash # Herre, need to check that keylist is filled i.e. files were found newrows = numpy.rec.array(rowtuple_list,names=keylist) tabhdu = tabHDUfromRecArray(newrows) if mode == 'new': tabhdu.writeto(summaryTable,clobber=True) nrows = tabhdu.header['NAXIS2'] summaryhash = hashFromPyfitsTable(tabhdu) else: newtabhdu = appendFitsTable(summaryTable,tabhdu) nrows = newtabhdu.header['NAXIS2'] summaryhash = hashFromPyfitsTable(newtabhdu) # Fake a goodrows hash sumgoodrows_hash = summaryhash.fromkeys(summaryhash.keys(),numpy.ones((nrows,),bool)) for filt in filters: print bandindex_hash[filt], filt plotStatsOnefile(summaryhash,sumgoodrows_hash,summaryTable,plotcol_summary_tlist,\ bandindex_hash[filt],bandid=filt,outroot="runstats_") mergePlots(summaryTable,plotcol_1frame_tlist=plotcol_1frame_tlist,\ filtlist=filters,workdir=workdir) return column_hash,goodrow_hash # XXX Todo: check histogram axis scaling! def hashFromPyfitsTable(pyfitsTableHDU): """Return a hash wrapping the table data so a column can be addressed as hash[colname] instead of tabledata.field(colname)""" from numpy import array outhash = {} for col in pyfitsTableHDU.columns.names: # Column names are uppercase, we want lowercase hash keys outhash[col.lower()] = array(pyfitsTableHDU.data.field(col)) return outhash def getOutnameStatsOnefile(matchtable,kind,col1_l,col2_l=None,format='eps',bandid='',useBandid=False): import re # Construct output filename root = re.sub('(\.[sc]mf|\.fits?)$','',matchtable) if useBandid and bandid != '': root = "%s_%s" %(root,bandid) if isNone(col2_l): outname = '%s_%s_%s.%s' % (root,kind,col1[0],format) else: if len(col1_l) == 1: outname = '%s_%s_%s_%s.%s' % (root,kind,col1_l[0],col2_l[0],format) else: outname = '%s_%s_%s_%s_%s_%s.%s' % (root,kind,col1_l[0],col2_l[0],col1_l[1],col2_l[1],format) return outname def plotStatsOnefile(values_hash,goodrow_hash,matchtable,plotcol_tlist,bandindex,\ format='epsport',bandid='',expand=0.8,outroot = "filestats_"): """Make diagnostic plots for a single table, based on values in values_hash""" from subprocess import call from numpy import concatenate from sm import cvar,angle import smLib import re import sm,smLib # In histograms, only plot core of distribution within these percentiles: histo_min_ntile = 0.03 histo_max_ntile = 1-histo_min_ntile angle_l = [0,45,0,45,0,90] ptype_l = [ 41, 41, 40, 40, 30,30] outname_l = [] nplots = len(plotcol_tlist) all_outname = outroot + re.sub('(\.[sc]mf|\.fits?)$','',matchtable) +".eps" smOpenPlot(all_outname,format=format) # all the explict sm. and smLib. calls should be wrapped into # functions that can be called for other plotting packages smSetTitle(re.sub('(\.[sc]mf|\.fits?)$','',matchtable)) i=0 for troika in plotcol_tlist: i+=1 Nx,Ny,winx,winy = setWindow(i,nplots) smLib.window(Nx,Ny,winx,winy) col1name_l = troika[0] col2name_l = troika[1] plottype = troika[2] outname = getOutnameStatsOnefile(matchtable,plottype,col1name_l,col2name_l,format=format,bandid=bandid) outname_l.append(outname) firstplot = True for col1name,col2name in zip(col1name_l,col2name_l): values1 = values_hash[col1name.lower()] goodrows1 = goodrow_hash[col1name.lower()] values2 = values_hash[col2name.lower()] goodrows2 = goodrow_hash[col2name.lower()] # 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)] if plottype == 'scatter': if firstplot: # Avoid plotting outliers if not re.search('objc',col2name): try: r = stats_med(values1[goodrows],[histo_min_ntile,histo_max_ntile]) [xmin,xmax] = r [ymin,ymax] = stats_med(values2[goodrows],[histo_min_ntile,histo_max_ntile]) except ValueError: print "Pffrz! ",r print col1name, goodrows, type(goodrows), len(goodrows), len(values1) # print "Huhu", outname, min(values1),max(values2),xmin,xmax smScatterPlot(values1,values2,logical=goodrows,xlab=col1name,ylab=col2name,\ # xrange=(xmin,xmax),yrange=(ymin,ymax)) xrange=None,yrange=None,expand=expand) else: smScatterPlot(values1,values2,logical=goodrows,xlab=col1name,ylab=col2name) else: angle(angle_l[i%len(angle_l)]) smScatterPlot(values1,values2,ptype=ptype_l[i%len(ptype_l)],\ logical=goodrows,append=True,expand=expand) angle(0) if plottype == 'histogram': Nbins = 20 values1 = values1[goodrows1] # This is SDSS column typically values2 = values2[goodrows2] # Avoid plotting outliers if not re.search('objc',col2name): [minbin,maxbin] = stats_med(concatenate((values1,values2)),\ [histo_min_ntile,histo_max_ntile]) else: minbin=None maxbin=None # Which one *am* I plotting first??? smHistoPlot(values2,ltype=0,nbins=Nbins,minbin=minbin,maxbin=maxbin,\ xlab=col2name+" \line 0 1000",ylab="N",expand=expand) smHistoPlot(values1,append=True,minbin=minbin,maxbin=maxbin,nbins=Nbins,\ ltype=2,expand=expand) smSetWindowTitle(col1name+" \line 2 1000") firstplot = False smClosePlot() return all_outname def smSetWindowTitle(titlestr): import sm, re sm.frelocate(0.5,1.05) sm.putlabel(5,re.sub('_','\_',titlestr)) def smSetTitle(titlestr): import sm sm.expand(0.8) sm.window(1,1,1,1) sm.location(3500,31000,3500,31000) smSetWindowTitle(titlestr) def setWindow(i,Nplots,portrait=True): """Set sm window for plotting the i-th out of Nplots. The logic is that N_x * N_y >= Nplots and N_y-N_x is 0 or 1. So Ny is the sqrt of N or the next-greatest square; Nx is either the same or one fewer if that's still enough.""" from math import sqrt if int(sqrt(Nplots))**2 == Nplots: Ny = int(sqrt(Nplots)) Nx = Ny else: Ny = int(sqrt(Nplots)+1) Nx = Ny-1 if Nx * Ny < Nplots: Nx+=1 # That's nx, ny; now count up to total Nx ny. This is like # representing i in a number system where the y-window number is # the 10s and the x-window number the 1s, but the 10s are base Ny # and the 1s are base Nx, and counting starts at 1 not 0. I.e. x = # (i % Nx) or Nx if that's 0; y = i / Nx + 1 winx = i% Nx if winx == 0: winx = Nx winy = int(i/Nx)+1 if int(i/Nx) == i/float(Nx): winy -= 1 if not portrait: tmp = Nx Nx = Ny Ny = tmp tmp = winx winx = winy winy = tmp return Nx,Ny,winx,winy def valuesKeysSortedByKeys(hash): """Return a list of hash's values and keys that is ordered by the names of the keys""" vallist = [] keylist = hash.keys() keylist.sort() for field in keylist: vallist.append(hash[field]) return vallist,keylist 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) newtabhdu = pyfits.new_table(oldtabhdulist[HDU].columns,nrows = nrows_old+len(tabHDU.data)) for col in newtabhdu.columns.names: newtabhdu.data.field(col)[nrows_old:] = tabHDU.data.field(col) oldtabhdulist.close() newtabhdu.writeto(fitsfile,clobber=True) return newtabhdu def tabHDUfromRecArray(recarr): """Generate a table HDU from a record array""" # create column names from recarray.dtype.names (NumPy style) # Assume float as data type # copy data field-wise to table HDU. import pyfits collist = [] for colname in recarr.dtype.names: if isinstance(recarr.field(colname)[0],int): format = '1J' elif isinstance(recarr.field(colname)[0],str): # This is potentially dangerous, but for now, the only # string column I am writing is the filter, i.e. 1A should # be enough already. Perhaps use the length of the entry # as column format? format = '1A' else: format = '1E' collist.append(pyfits.Column(name=colname,array=recarr.field(colname),format=format)) coldefs = pyfits.ColDefs(collist) return pyfits.new_table(coldefs) def mkRecordArray(hash): """Return a record array representing a table row from a hash; keys are column names, values are data""" newrow = numpy.rec.array([tuple(hash.values())]) newrow.dtype.names=hash.keys() # Not : # newrow = numpy.rec.array([tuple(hash.values())],names=hash.keys())? return newrow def saveHierarchHeaderList(header,label_l,value_l): """Given paired lists label_l and value_l, save the values from value_l to ESO HIERARCH header keywords with the names given in label_l""" for label,value in zip(label_l,value_l): header.update('HIERARCH %s'%(label),value) def stats_med(arr,ntiles=[0.25,0.5,0.75]): """Return a list of ntiles of array arr""" from numpy import nan,sort l = len(arr) if l == 0: return nan,nan,nan sortarr = sort(arr,kind='mergesort') outl = [] for ntile in ntiles: idx = max(int(ntile*(l-1)),0) outl.append(sortarr[idx]) return outl def computeStatistics(tablename,copyfields_list = ['RUN','RERUN','CAMCOL','FIELD','FILTER','FWHM_MAJ','FWHM_MIN']): """Compute desired statistics for offsets etc. and add them as ESO-style HIERARCH fields to the table header (8 characters are otherwise too few to contain meaningful names for the keywords). Returns a hash with all the computed statistics, keyed by their header names""" import pyfits,re,operator from numpy import sqrt,log10,log,array,isfinite label_l = ['mean','rms','median','lowerquartile','upperquartile'] filterID={'u':0,'g':1,'r':2,'i':3,'z':4} infile_handle = pyfits.open(tablename,mode='update') table = infile_handle[1] h = table.header table_data = table.data #delete h0 = infile_handle[0].header # Header names that are going to be copied to output table as # "reference columns" outhash = {} 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)) # Maybe I want to keep track of these number of "good" rows? # Count number of objects in a) SDSS, b) PS1, c) both, by # counting number of a) entries > 0 in 'id', b) non-nan entries # in IPP_IDET, c) both (Later pass the following as parameters to # be read from a config file) N_SDSS_col = 'id' N_SDSS_outcolname = 'N_SDSS' N_PS1_col = 'IPP_IDET' N_PS1_outcolname = 'N_PS1' N_both_outcolname = 'N_both' N_either_outcolname = 'N_either' has_sdss = array(table_data.field(N_SDSS_col)) > 0 has_PS1 = isfinite(array(table_data.field(N_PS1_col))) outhash[N_SDSS_outcolname] = sum(has_sdss) outhash[N_PS1_outcolname] = sum(has_PS1) outhash[N_both_outcolname] = sum(has_sdss & has_PS1) outhash[N_either_outcolname] = len(has_PS1) # Hash of output header keyword names pointing to paired lists of # columns that are to be subtracted from each other for statistics # to be generated. This lists just the columns, slicing out the # correct depth from SDSS 5-filter arrays will be done later. colname_hash = { 'd_x':['colc','X_PSF'] ,'d_y':['rowc','Y_PSF'] ,'d_sky':['sky_sdss','SKY_ps1'] ,'d_skyerr':['skyErr','SKY_SIGMA'] ,'d_pointsource':['prob_psf','pointsource_ps1'] ,'d_mag':['psfinstmag_sdss','PSF_INST_MAG'] ,'d_magerr':['psfinstmagerr_sdss','PSF_INST_MAG_SIG'] } # colname_hash should probably be passed as a parameter so it can # be read from file. ismag = re.compile('mag') iscounts = re.compile('counts') # outcoll = [] for outcol,list in colname_hash.iteritems(): SDSScolname,PS1colname=list # 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 # instmagclip.py was here delta = SDSScol_good - PS1col_good # 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) rms = sqrt(delta.var()) # Save to fits columns/header 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,[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]): outhash[label]=value # #t lines were for writing the new columns which I'm now creating in the match script #t newtab = pyfits.new_table(table.columns+pyfits.ColDefs(outcoll),header=h) #t newprimhdu = pyfits.PrimaryHDU(header=infile_handle[0].header) infile_handle.close() #t writeTable(tablename,newprimhdu,newtab) return outhash,colval_hash,goodrow_hash def writeTable(filename,primhdu,tabhdu): """Write pyfits.ColDefs thing columns to filename, clobbering existing file""" import pyfits hdulst = pyfits.HDUList([primhdu,tabhdu]) hdulst.writeto(filename,clobber=True) def goodValBool(col): """Return a bool array that contains True where col 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 goodcondition = goodValBool(col1) & goodValBool(col2) col1 = col1[goodcondition] col2 = col2[goodcondition] return col1,col2 def filterGoodVal3(col1,col2,col3): """Return a triplet 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 goodcondition = goodValBool(col1) & goodValBool(col2) & goodValBool(col3) col1 = col1[goodcondition] col2 = col2[goodcondition] col3 = col3[goodcondition] return col1,col2,col3 def matchSdssPs1(SDSSfpObjc,PS1cmf,xoff=0.5,yoff=0.5,matchrad=0.7,skipMatch=False): """Call matchByPos to match an SDSS fpObjc.fits and a PS1 bla.cmf file.""" import pyfits from numpy import nan def getOutname(SDSSfile,PS1file,sdssbandstr): import re # Construct output filename fitsend = re.compile('(\.[sc]mf|\.fits?)$') SDSSroot = fitsend.sub('',SDSSfile) PS1root = fitsend.sub('',PS1file) # Now chop off any leading path components pathroot = re.compile('.*/') SDSSroot = pathroot.sub('',SDSSroot ) PS1root = pathroot.sub('',PS1root) return "match_%s___%s___%s.fits" % (sdssbandstr,SDSSroot,PS1root) filters = ['u','g','r','i','z'] # Read primary header of PS1 file to work out band and copy other # interesting fields. These include FHWM_MAJ,MIN which used to be # called FWHM_X,Y, but are missing in some cases if the camera # stage failed. In the latter case, we can still read IQ_FW1,2. So # we first update the .cmf file to be sure that FWHM_MAJ,MIN is # present for later analysis. ensureFHWMheaderfields(PS1cmf) ps1copyfields_hash = headerfieldHash(['FWHM_MAJ','FWHM_MIN','FILTER'],PS1cmf,HDU=0) sdsscopyfields_hash = headerfieldHash(['RUN','RERUN','CAMCOL','FIELD'],SDSSfpObjc,0) sdssbandstr = ps1copyfields_hash['FILTER'] bandindex = filters.index(sdssbandstr) outname = getOutname(SDSSfpObjc,PS1cmf,sdssbandstr) if skipMatch: return outname,sdssbandstr,bandindex # Note position offset - y_sdss-y_psf1 ~ 0.6, x_sdss-x_ps1 ~ 0.35 (in 1056-r1-0192) # width of distribution is ~ 0.4 # So subtract 0.5 from SDSS numbers and match with r=0.7 should include all good matches # Arghly, need to match datatypes and 0.5 is double by default # while cols are float, hence the typecast. sdssfilt = 'addcol colc_%s colc[%s]-(float)%s' % (sdssbandstr,bandindex,xoff) sdssfilt += ';addcol rowc_%s rowc[%s]-(float)%s' % (sdssbandstr,bandindex,yoff) sdssfilt += ';addcol psfinstmag_sdss -2.5*log10(psfcounts[%s])' % (bandindex) sdssfilt += ';addcol psfinstmagerr_sdss 2.5/ln(10)*psfcountserr[%s]/psfcounts[%s]' %(bandindex,bandindex) # Filtering out objects that can't be primary, i.e. # !BRIGHT && (!BLENDED || NODEBLEND || nchild == 0) # Get flag values by left-shifting a 1 by the correct number of bits. Bit N needs to be left-shifted N times. bright = 1<<1 blended = 1<<3 nodeblend = 1<<6 binned1 = 1<<28 sdssfilt += ';select ((flags[%s]&%s)==0)&&(!((flags[%s]&%s)!=0)||((flags[%s]&%s)!=0)||nchild==0)&&((flags[%s]&%s)!=0)' % \ (bandindex,bright,bandindex,blended,bandindex,nodeblend,bandindex,binned1) # Check flag 0x4000 (extended) to compare against prob_psf in SDSS # - for some reason, 'flags' and the other integer columns get # read as float by topcat, so need a cast here ps1filt = 'addcol pointsource_ps1 (((int)flags&0x4000)==0?1:0)' sdsspos = '\'colc_%s rowc_%s\''% (sdssbandstr,sdssbandstr) ps1pos = "\'x_psf y_psf\'" matchByPos(SDSSfpObjc,PS1cmf,outname,sdsspos,ps1pos,tolerance=matchrad,\ duptag1='_sdss',duptag2='_ps1',\ filter1='\''+sdssfilt+'\'',filter2='\''+ps1filt+'\'') f=pyfits.open(outname,mode='update') h=f[1].header updateHeaderFromHash(h,sdsscopyfields_hash) updateHeaderFromHash(h,ps1copyfields_hash) f.close() return outname,sdssbandstr,bandindex def ensureFHWMheaderfields(PS1cmf): import pyfits PS1handle = pyfits.open(PS1cmf,mode='update') PS1primhead = PS1handle[0].header headfields = PS1primhead.ascardlist().keys() if not ('FWHM_MAJ' in headfields and 'FWHM_MIN' in headfields): if ('FWHM_X' in headfields and 'FWHM_Y' in headfields): PS1primhead.update('FWHM_MAJ',PS1primhead['FWHM_X']) PS1primhead.update('FWHM_MIN',PS1primhead['FWHM_Y']) else: PS1primhead.update('FWHM_MAJ',PS1primhead['IQ_FW1']) PS1primhead.update('FWHM_MIN',PS1primhead['IQ_FW2']) PS1handle.close() def updateHeaderFromHash(header,hash): """Add fields in hash to header""" for k in hash.keys(): header.update(k,hash[k]) def headerfieldHash(fieldlist,file,HDU=0): """Return a hash containing listed fits headers from file's HDU (this is like a hash slice, which doesn't seem to be possible in python?)""" import pyfits f=pyfits.open(file) h=f[HDU].header outhash = {} for field in fieldlist: outhash[field] = h[field] f.close() return outhash def matchByPos(in1,in2,out,colnames1,colnames2,tolerance=0.5,duptag1='_sdss',duptag2='_ps1',filter1="",filter2=""): """Match two fits tables by cartesian 2-d position using stilts/tmatch2 from topcat""" import subprocess # tmatch_args - sequence of parameters is # infile1, infile2, outfile, x_col, y_col, tolerance [all in pixels], duptag1,2 # Looking for best match only tmatch_args = 'in1=%s in2=%s out=%s values1=%s values2=%s params=%s fixcols=dups suffix1=%s suffix2=%s' tmatch_args += ' matcher=2d ifmt1=fits ifmt2=fits ofmt=fits-basic join=1or2 find=best' if filter1 != "": tmatch_args += " icmd1=%s" %(filter1) if filter2 != "": tmatch_args += " icmd2=%s" %(filter2) tmatch_cmd = tmatch_args % (in1,in2,out,colnames1,colnames2,tolerance,duptag1,duptag2) tmatch_cmd = os.environ.get('HOME')+'/bin/stilts tmatch2 '+tmatch_cmd print tmatch_cmd retval = subprocess.call(tmatch_cmd,shell=True) def smdefaults(): # Will these be persistent across multiple imports and going out # of scope? Perhaps only if I do a global 'import sm'. Or it # should work if called from something where sm is in scope? sm.expand(1.8) sm.lweight(3) # Plotting convention will be like in sm: you need to open a file, # then you plot to it with one or multiple commands, then you close it # to write the file. def smOpenPlot(filename,format='eps'): """Issue sm device command for file of given format: eps -> postfile epsland -> postlandfile epsport -> postportfile else just use given format as 'device'""" from sm import device, erase if filename == 'x11': device('x11') elif format == 'eps': device('postfile '+filename) elif format == 'epsland': device('postlandfile '+filename) elif format == 'epsport': device('postportfile '+filename) else: device(format +' '+filename) erase() def smClosePlot(): from sm import device device('nodevice') def isNone(something): return isinstance(something,type(None)) def smHistoPlot(vec,step=None,minbin=None,maxbin=None,nbins=None,\ xlab=None,ylab=None,xrange=None,yrange=None,\ box1=None,box2=None,box3=None,box4=None,\ ltype=0,append=False,expand=1.8): """Plot a histogram, intelligently deriving bins from the given parameters if they are given intelligently. Otherwise, silently do nothing.""" import sm from numpy import histogram if isNone(minbin): minbin = min(vec) if isNone(maxbin): maxbin = max(vec) # I am assuming bins are bin centers if not isNone(step): nbins = (maxbin-minbin)/step+1 if isNone(nbins): return histo,binedges = histogram(vec,nbins,[minbin,maxbin],new=True) bincenters = binedges + 0.5*(binedges[1]-binedges[0]) bincenters = bincenters[0:-1] sm.expand(expand) sm.ltype(ltype) if not append: logical = None smSetup(bincenters,histo,logical,xrange,yrange,xlab,ylab,box1,box2,box3,box4) sm.histogram(bincenters,histo) sm.ltype(0) return minbin,maxbin 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 line plot on current device. If append=True, overplot with current limits. Otherwise, draw box box1 box2 box3 box4""" import sm sm.expand(1.8) sm.lweight(3) # Silently ignore any problems with plot generation try: if not append: smSetup(x,y,logical,xrange,yrange,xlab,ylab,box1,box2,box3,box4) sm.ltype(ltype) sm.connect(x,y,logical) sm.ltype(0) except: pass 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,expand=1.8): """Make an sm scatter plot on current device. If append=True, overplot with current limits. Otherwise, draw box box1 box2 box3 box4""" import sm sm.expand(expand) sm.lweight(3) # Silently ignore any problems with plot generation try: if not append: smSetup(x,y,logical,xrange,yrange,xlab,ylab,box1,box2,box3,box4) sm.ptype(ptype) sm.points(x,y,logical) except: pass # May need to add erase option here which defaults to True def smSetup(x,y,logical,xrange,yrange,xlab,ylab,box1,box2,box3,box4): import sm,re sm.erase() if not isNone(logical): x=x[logical] y=y[logical] if isNone(xrange): xrange = x if isNone(yrange): yrange = y # print "Setting limits to ",xrange,yrange sm.limits(xrange,yrange) smBox(box1,box2,box3,box4) if not isNone(xlab): sm.xlabel(re.sub('_','\_',xlab)) if not isNone(ylab): sm.ylabel(re.sub('_','\_',ylab)) def smBox(box1,box2,box3,box4): from sm import box if not isNone(box4): box(box1,box2,box3,box4) elif not isNone(box3): box(box1,box2,box3) elif not isNone(box2): box(box1,box2) elif not isNone(box1): box(box1) else: box() def gaussResid(p,yx,data,blankRadius=0,blankInner=True): """return residuals after fitting two-D gaussian with parameters p = scale,x0,xsig, y0, ysig, theta, blanking pixels inside or outside blankRadius, depending on whether blankInner = True or False""" from numpy import reshape scale,x0,xsig,y0,ysig,theta = p model=scale*gauss2D(yx,x0,xsig,y0,ysig,theta) resid = reshape(data,(data.shape[0]*data.shape[1],))-reshape(model,(model.shape[0]*model.shape[1],)) if (blankRadius > 0): # set residuals within or outside blankRadius to 0, i.e. set model = data there x=yx[1] y=yx[0] if blankInner: inradius = (((x-x0)**2 + (y-y0)**2) <= blankRadius**2) else: inradius = (((x-x0)**2 + (y-y0)**2) > blankRadius**2) inradius = reshape(inradius,(inradius.shape[0]*inradius.shape[1],)) resid[inradius] = 0 # win.showArray(reshape(resid,model.shape)) return resid def fitGauss(yx,data,p0,blankrad,blankInner): """fit gaussian to data with starting parameters in p0""" from scipy.optimize import leastsq from numpy import reshape # need to add weighting by errors, if we know them plsq = leastsq(gaussResid,p0,args=(yx,data,blankrad,blankInner)) bestp = plsq[0] # print plsq[1] scale,x0,xsig,y0,ysig,theta = bestp # remove 360 degree ambivalence theta -= int(theta/180)*180 # make sure xsig is major axis if xsig < ysig: tmp=xsig xsig=ysig ysig=tmp theta += 90 # make sure theta is in [-90,90] if theta > 90: theta -= 180 bestp=(scale,x0,xsig,y0,ysig,theta) resid=reshape(gaussResid(bestp,yx,data,blankrad,blankInner),data.shape) mod=data-resid return bestp,resid,mod def gauss2D(coordgrid,x0,xsig,y0,ysig,theta=0): """compute a normalized gaussian at positions given in 2d-array coords, centered at x0, y0, with xsig, ysig for theta=0; coordinates in coordgrid given by ogrid""" from numpy import cos,sin,pi,exp cth = cos(theta/180.*pi) sth = sin(theta/180.*pi) x=coordgrid[1]-x0 y=coordgrid[0]-y0 res = exp(-0.5*(((x*cth+y*sth)/xsig)**2+((-x*sth+y*cth)/ysig)**2)) return res/(2*pi*xsig*ysig) def doGaussPsfPhotometry(infile,initialGuessList=[(2000.,1011.,2.,1176.,2.,0)],\ blankrad=9,defHalfwinsize=15,HDU=0): """Fit Gaussians to image in infile's HDU near the positions and with fluxes and sigmas given in initialGuessList""" import pyfits from numpy import nan from numarray import ones,zeros,Float32 from scipy import ogrid f_handle=pyfits.open(infile) img = f_handle[HDU].data resid_img = zeros(img.shape,Float32) fitlist = [] for initialGuessTuple in initialGuessList: scale0,x0,xsig0,y0,ysig0,theta0 = initialGuessTuple # Allow for variable blankrad leading to variable halfwinsize # (no information beyond blankrad will be used), but avoid going to 0 if blankrad > 0: halfwinsize = blankrad else: halfwinsize = defHalfwinsize subim = img[int(y0)-halfwinsize:int(y0)+halfwinsize+1,int(x0)-halfwinsize:int(x0)+halfwinsize+1] # ds9-based coords start at 1,1, while numpy's start at 0,0, # so we need to add 1 to all internally used coords to get the # fitted coordinates to match with what we see in ds9 yx=ogrid[1+int(y0)-halfwinsize:1+int(y0)+halfwinsize+1,1+int(x0)-halfwinsize:1+int(x0)+halfwinsize+1] try: par,res,bestfit=fitGauss(yx,subim,initialGuessTuple,blankrad,blankInner=False) except: par=(nan,nan,nan,nan,nan,nan) res = subim (scale,x,xsig,y,ysig,theta)=par # 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 fitlist.append(par) f_handle.close() outname = 'gaussfit_resid_r%ipix_%s' % (blankrad,infile) primHDU = pyfits.PrimaryHDU(resid_img) hdulist = pyfits.HDUList([primHDU]) hdulist.writeto(outname,clobber=True) return fitlist def doEfbComparison(simfile='simtest.1.00a.fits',efbfile='efbsim.1.00a.fits',simtab='simtest.1.00a.cmf',blankrad=9): import pyfits from numpy import array, rec import re tab = pyfits.open(simtab) tdata = tab[1].data x = array(tdata.field('x_psf')) y = array(tdata.field('y_psf')) counts = 10**(-0.4*array(tdata.field('psf_inst_mag'))) initlist = [] inparlist = [] nrows = len(tdata.field('IPP_IDET')) for i in range(nrows): initlist.append((counts[i],x[i],2.,y[i],2.,0.)) inparlist.append((counts[i],x[i],y[i])) efbparlist = doGaussPsfPhotometry(efbfile,initlist,blankrad=blankrad) simparlist = doGaussPsfPhotometry(simfile,initlist,blankrad=blankrad) savelist = [ inparlist[i] + efbparlist[i] + simparlist[i] for i in range(nrows)] saverows = rec.array(savelist,names=['counts_in','x_in','y_in',\ 'efbfit_counts','efbfit_x','efbfit_fwhm_x',\ 'efbfit_y','efbfit_fwhm_y','efbfit_theta',\ 'simfit_counts','simfit_x','simfit_fwhm_x',\ 'simfit_y','simfit_fwhm_y','simfit_theta']) for col in ['efbfit_fwhm_x','efbfit_fwhm_y','simfit_fwhm_x','simfit_fwhm_y']: saverows.field(col)[0:] *= 2.35 saveHDU = tabHDUfromRecArray(saverows) fitsend = re.compile('(\.[sc]mf|\.fits?)$') saveHDU.writeto('gaussfit_r%ipix_%s___%s' %(blankrad,fitsend.sub('',efbfile),simfile),clobber=True) def makePlan(fpObjcDir = '/IPP/data/SDSS/stripe82/coadd/input/fpObjcs/',\ cmfDir = '/IPP/data/SDSS/stripe82/coadd/prod/run_ipp_20080815/',\ lastcamcol=6,\ runlist=[1056,1755,3388,3434,3465,4145,4192,4203,4247,4263,5052],\ firstfield_h = { 1056:192, 1755:329, 3388:273, 3434:306, 3465:185, 4145:63, 4192:309, 4203:358, 4247:62, 4263:69, 5052:62 }, \ Nfields = 18): """Make paired list of which fpObjc file to compare to which IPP .cmf file. Could also think about constructing the fpObjc filename from run, camcol, field read from .cmf primary header, which replicates most of the fpC primary header.""" if cmfDir[-1] != '/': cmfDir += '/' if fpObjcDir[-1] != '/': fpObjcDir += '/' # For testing: # return ['1056-0192.421/1056-0192.421.ch.421.CHIP1.cmf'],['1056-0192.421/fpObjc-001056-1-0192.fit'] # return ['1056-0192.1/1056-0192.1.ch.1.CHIP21.cmf'],['1056-0192.1/fpObjc-001056-1-0192.fit'] from glob import glob camcollist = range(1,lastcamcol+1) sdsslist = [] ps1list = [] for run in runlist: firstfield = firstfield_h[run] for field in range(firstfield,firstfield+Nfields+1): for camcol in camcollist: sdssname = fpObjcDir + "fpObjc-%06i-%i-%04i.fit" % ( run, camcol, field) # The following should just return 5 files (1 per # filter), but just to be sure, let's check and return # the correct number of sdss files to match: globstr = cmfDir + "%i-%04i.*/%i-%04i.*.ch.*.CHIP%i.cmf" %(run,field,run,field,camcol) ps1glob = glob(globstr) for ps1name in ps1glob: ps1list.append(ps1name) sdsslist.append(sdssname) # Need to remember that CHIPx corresponds to CAMCOLx; could change # naming rule to make filter name appear in filename, too. return ps1list,sdsslist