# batch analyse Nicolas' Upper Sco (deep) stack tiles from pyraf import iraf from iraf import noao, digiphot, daophot, ptools, images, \ imcoords, tables, ttools, fitsutil, imutil import math, time, os, bisect, random # CASU VIRCam median parameters for tiles: zGain = 4.1775 zRon = 24.7615 zFlatErr = 0.0155 yGain = 4.1505 yRon = 21.545 yFlatErr = 0.0155 jGain = 4.3305 jRon = 25.649 jFlatErr = 0.0145 # list of tuples of VSA RICE compressed tile images to work through, with attributes: # VSA filename, gain, ron, flaterr, shapeFrac, passbandEpoch, mag range, nImages ... # required MF IDs: # Y2: 1000000043062, 1000000043058, 1000000043057, 1000000043064 # J1: 1207618, 1256433, 1252117, 1252073, # 1207574, 1249768, 1208247, 1208291, 1208345 (2nd priority) # Y1: 1206011, 1203470, 1205113, 1203586, # Z1: 1555413, 1555415 tileFiles = [\ ('v20120415_00675_st_tl.fit', yGain, yRon, yFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 10, 111),\ ('v20120417_00412_st_tl.fit', yGain, yRon, yFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 10, 131),\ ('v20120417_00442_st_tl.fit', yGain, yRon, yFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 10, 132),\ ('v20120417_00472_st_tl.fit', yGain, yRon, yFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 10, 133),\ ('v20120419_00392_st_tl.fit', yGain, yRon, yFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 10, 121),\ ] # Priority 5 Z: # ('e20170925_10017000016_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Z1', ' && MAG < -3.0 && MAG > -3.5', 20, 123),\ # ('e20170925_10017000018_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Z1', ' && MAG < -3.0 && MAG > -3.5', 20, 132),\ # ('e20170925_10017000013_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Z1', ' && MAG < -3.0 && MAG > -3.5', 20, 121),\ # ('e20170925_10017000004_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Z1', ' && MAG < -3.0 && MAG > -3.5', 20, 111),\ # ('e20170925_10017000001_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Z1', ' && MAG < -3.0 && MAG > -3.5', 20, 112),\ # ('e20170925_10017000010_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Z1', ' && MAG < -3.0 && MAG > -3.5', 20, 122),\ # ('e20170925_10017000024_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Z1', ' && MAG < -3.0 && MAG > -3.5', 20, 133),\ # ('e20170925_10017000007_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Z1', ' && MAG < -3.0 && MAG > -3.5', 20, 113),\ # ('e20170925_10017000021_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Z1', ' && MAG < -3.0 && MAG > -3.5', 20, 131),\ # Priority 5 (second try): # ('v20120422_00329_st_tl.fit', jGain, jRon, jFlatErr, 1.0, 'J1', ' && MAG < -3.0 && MAG > -3.5', 10, 131),\ # ('v20120422_00395_st_tl.fit', jGain, jRon, jFlatErr, 1.0, 'J1', ' && MAG < -3.0 && MAG > -3.5', 10, 133),\ # ('v20120421_00353_st_tl.fit', jGain, jRon, jFlatErr, 1.0, 'J1', ' && MAG < -3.0 && MAG > -3.5', 10, 121),\ # ('v20120422_00359_st_tl.fit', jGain, jRon, jFlatErr, 1.0, 'J1', ' && MAG < -3.0 && MAG > -3.5', 10, 132),\ # Priority 5 (first try!): # ('v20120420_00601_st_tl.fit', jGain, jRon, jFlatErr, 1.0, 'J1', ' && MAG < -3.0 && MAG > -3.5', 10, 111),\ # Priority 1-4: # ('e20170227_10052000001_dp_st_tl.fit', yGain, yRon, yFlatErr, 0.5, 'Y2', ' && MAG < -3.0 && MAG > -3.5', 140, 111),\ # ('e20170227_10052000002_dp_st_tl.fit', yGain, yRon, yFlatErr, 0.5, 'Y2', ' && MAG < -3.0 && MAG > -3.5', 140, 112),\ # ('e20170227_10052000003_dp_st_tl.fit', yGain, yRon, yFlatErr, 0.5, 'Y2', ' && MAG < -3.0 && MAG > -3.5', 140, 122),\ # ('e20170227_10052000004_dp_st_tl.fit', yGain, yRon, yFlatErr, 0.5, 'Y2', ' && MAG < -3.0 && MAG > -3.5', 140, 121),\ # ('v20120420_00631_st_tl.fit', jGain, jRon, jFlatErr, 1.0, 'J1', ' && MAG < -3.0 && MAG > -3.5', 10, 112),\ # ('v20120428_00641_st_tl.fit', jGain, jRon, jFlatErr, 1.0, 'J1', ' && MAG < -3.0 && MAG > -3.5', 10, 123),\ # ('v20120428_00671_st_tl.fit', jGain, jRon, jFlatErr, 1.0, 'J1', ' && MAG < -3.0 && MAG > -3.5', 10, 113),\ # ('v20120430_00561_st_tl.fit', jGain, jRon, jFlatErr, 1.0, 'J1', ' && MAG < -3.0 && MAG > -3.5', 10, 122),\ # ('v20120416_00645_st_tl.fit', yGain, yRon, yFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 10, 122),\ # ('v20120416_00719_st_tl.fit', yGain, yRon, yFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 10, 123),\ # ('v20120418_00449_st_tl.fit', yGain, yRon, yFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 10, 113),\ # ('v20120419_00422_st_tl.fit', yGain, yRon, yFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 10, 112),\ # ouch ! First run of this script had this typo in it, so rerun......V......... # ('e20130416_1001700010_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 20, 122),\ # ('e20130416_1001700016_dp_st_tl.fit', zGain, zRon, zFlatErr, 1.0, 'Y1', ' && MAG < -3.0 && MAG > -3.5', 20, 123),\ #] ######################################################################################## iraf.datapars.fwhmpsf=2.6 iraf.datapars.scale=1. iraf.findpars.threshold=2.5 iraf.datapars.gain="" iraf.datapars.ccdread="" iraf.photpars.apertures="3." iraf.photpars.zmag=25.0 iraf.daopars.function='penny2' iraf.daopars.varorder=2 iraf.daopars.nclean=10 iraf.daopars.psfrad=11. iraf.daopars.proferr=5.0 iraf.daopars.fitrad= 2.0 * iraf.datapars.fwhmpsf nReqPsfStars = 100 # number of required PSF stars def filterPsfCandidates(allDetectionsCooFile, psfCandidatesCooFile, radius): # open and read all detections: xyList = [] f = open(allDetectionsCooFile, mode='r') line = f.readline() while line.startswith('#'): line = f.readline() n = 0 while line != '': #print line n = n + 1 items = line.split() x = float(items[0]) y = float(items[1]) xyList.append((x,y)) line = f.readline() f.close() print 'Read ',n,' lines from ',allDetectionsCooFile # sort on increasing Y: xyList.sort(key=lambda r: r[1]) keys = [r[1] for r in xyList] # now open up the candidates file, copy the header into a new one and examine each # entry for proximity of detection(s) within the specified radius: f = open(psfCandidatesCooFile, mode='r') outputFileName = psfCandidatesCooFile.replace('.coo','_iso.coo') g = open(outputFileName, mode='w') h = open(psfCandidatesCooFile.replace('.coo','_rej.coo'), mode='w') d2lim = radius * radius line = f.readline() while line.startswith('#'): g.write(line) h.write(line) line = f.readline() n = 0 nrej = 0 nkept = 0 while line != '': n = n + 1 items = line.split() x = float(items[0]) y = float(items[1]) # search the list for this object: idx = bisect.bisect_left(keys, y) inxt = 1 dNearest = radius # search forwards: if idx + inxt < len(xyList): while idx + inxt < len(xyList) and xyList[idx+inxt][1] - y <= dNearest: dx = xyList[idx+inxt][0] - x dy = xyList[idx+inxt][1] - y d = math.sqrt(dx*dx + dy*dy) if d < dNearest: dNearest = d idxNearest = idx+inxt inxt = inxt + 1 # search backwards: inxt = 1 if idx - inxt >= 0: while idx - inxt >= 0 and y - xyList[idx-inxt][1] <= dNearest: dx = xyList[idx-inxt][0] - x dy = xyList[idx-inxt][1] - y d = math.sqrt(dx*dx + dy*dy) if d < dNearest: dNearest = d idxNearest = idx-inxt inxt = inxt + 1 if dNearest >= radius: g.write(line) nkept = nkept + 1 else: h.write(line) nrej = nrej + 1 line = f.readline() print 'Rejected ',nrej,' PSF candidates with proximal sources and kept ', nkept f.close() g.close() h.close() return outputFileName, nkept for tileFile in tileFiles: # set up IRAF for this image from the attributes in the list of tuples: nImage = tileFile[7] gain = tileFile[1] ron = tileFile[2] flatErr = tileFile[3] magRange = tileFile[6] shapeFraction = tileFile[4] tileId = tileFile[8] #iraf.images.imgets(image=tileFile,param='OBJECT') #object = str(imgets.value) #print 'OBJECT keyword is ',object # image dependent paramaters: iraf.datapars.epadu = gain * nImage iraf.datapars.readnoise = ron / math.sqrt(nImage) iraf.daopars.flaterr = flatErr / math.sqrt(nImage) # inflate: funpack originFName = tileFile[0] newName = originFName.replace('.fit','.fits.fz') os.rename(originFName, newName) iraf.fitsutil.funpack(images=newName) fitsName = newName.replace('.fz','') # split into lo/hi parts: subSetFiles = [fitsName.replace('.fits','_lo.fits'), fitsName.replace('.fits','_hi.fits')] iraf.images.imutil.imcopy(input=fitsName+'[1000:11763,200:7850]', output=subSetFiles[0]) iraf.images.imutil.imcopy(input=fitsName+'[1000:11763,7800:15449]', output=subSetFiles[1]) for fileName in subSetFiles: print '\n\n Analysing '+fileName # find standard deviation of pixel values: startTime = time.time() s = iraf.imstat(fileName, fields='mean,stddev', lower=100, upper=10000, nclip=10, Stdout=1) print 'IMSTAT took ', int(time.time() - startTime), ' seconds.' stddev=float(s[1].split()[1]) iraf.datapars.sigma=stddev meanSky=float(s[1].split()[0]) print 'Mean sky level (ADU): ',meanSky,' +/- ',stddev # off we go: DAOFIND: startTime = time.time() daofindOutputFile = fileName.replace('.fits', '_daofind.coo') iraf.noao.digiphot.daophot.daofind(image=fileName, output=daofindOutputFile, \ verify='no', update='no', verbose='no') print 'DAOFIND took ', int(time.time() - startTime), ' seconds.' # work out the ranges for SHARPNESS, SROUND and GROUND: s = iraf.tables.ttools.tstat(intable=daofindOutputFile, column=4, Stdout=1) medSharp=float(s[2].split()[3]) sdevSharp=float(s[2].split()[2]) s = iraf.tables.ttools.tstat(intable=daofindOutputFile, column=5, Stdout=1) medSRound=float(s[2].split()[3]) sdevSRound=float(s[2].split()[2]) s = iraf.tables.ttools.tstat(intable=daofindOutputFile, column=6, Stdout=1) medGRound=float(s[2].split()[3]) sdevGRound=float(s[2].split()[2]) # filter DAOFIND output on sharpness and roundness to create an initial point # source list from which to get the PSF candidates pointSourceFile = daofindOutputFile.replace('daofind','pntsrc') expression = 'SHARPNESS>'+str(medSharp - sdevSharp*shapeFraction) + \ ' && SHARPNESS<'+str(medSharp + sdevSharp*shapeFraction) + \ ' && SROUND>'+str(medSRound - sdevSRound*shapeFraction) + \ ' && SROUND<'+str(medSRound + sdevSRound*shapeFraction) + \ ' && GROUND>'+str(medGRound - sdevGRound*shapeFraction) + \ ' && GROUND<'+str(medGRound + sdevGRound*shapeFraction) + magRange print 'Filtering DAOFIND detection list with ' + expression iraf.noao.digiphot.ptools.pdump(infile=daofindOutputFile, \ fields='XCENTER,YCENTER,MAG,SHARPNESS,SROUND,GROUND,ID', expr=expression, \ headers='yes', parameters='yes', Stdout=pointSourceFile) # do the proximity filtering: newPointSourceFile, nCand = filterPsfCandidates(daofindOutputFile, pointSourceFile, 2.0*iraf.daopars.psfrad) if nCand < 50: print 'ERROR: found too few psf candidates: ',nCand continue # do the aperture photometry on both full file and the initial selection of point sources photOutputFile = daofindOutputFile.replace('coo','mag') pointSourcePhotOutputFile = newPointSourceFile.replace('coo','mag') iraf.noao.digiphot.daophot.phot(image=fileName, coords=newPointSourceFile, \ output=pointSourcePhotOutputFile, cache='no', \ inter='no', verify='no', update='no', verbose='no') # edit down the previous list for a list of PSF stars, avoiding any where higher sky may # indicate presence of nearby star(s) and/or proximity to very bright star etc. psfCandidatesTmpFile = pointSourcePhotOutputFile.replace('mag','tmp') iraf.noao.digiphot.ptools.pdump(infile=pointSourcePhotOutputFile, fields='ID,XCENTER,YCENTER,MAG,MSKY', \ expr='MSKY<'+str(meanSky + stddev), headers='yes', parameters='yes', Stdout=psfCandidatesTmpFile) # take 1-in-N stars in order to leave a random selection of 100: psfCandidatesFile = psfCandidatesTmpFile.replace('tmp','pst') thresh = float(nReqPsfStars) / float(nCand) f = open(psfCandidatesTmpFile, mode='r') g = open(psfCandidatesFile, mode='w') line = f.readline() nRandomlySelected = 0 while line.startswith('#'): g.write(line) line = f.readline() while line != '': if random.random() < thresh or nCand < 1.2*nReqPsfStars: g.write(line) nRandomlySelected = nRandomlySelected + 1 line = f.readline() f.close() g.close() print 'Randomly selected ',nRandomlySelected,' PSF stars.' # pdump them so we can tvmark to check later iraf.noao.digiphot.ptools.pdump(infile=psfCandidatesFile, fields='XCENTER,YCENTER', \ expr='yes', headers='no', parameters='no', \ Stdout=psfCandidatesFile.replace('.pst','_filtered.coo')) #continue # TODO remove after checking candidate selection numbers etc. ! # do the big aperture photometry run startTime = time.time() iraf.noao.digiphot.daophot.phot(image=fileName, coords=daofindOutputFile, \ output=photOutputFile, cache='no', \ inter='no', verify='no', update='no', verbose='no') print 'PHOT took ',int(time.time() - startTime), ' seconds.' # now do the PSF using the filtered point source list psfImageFile = psfCandidatesFile.replace('pst','psf.fits') opstFile = pointSourcePhotOutputFile.replace('mag','opst') groupFile = psfCandidatesFile.replace('pst','psg') startTime = time.time() iraf.noao.digiphot.daophot.psf(image=fileName, \ photfile=pointSourcePhotOutputFile, psfimage=psfImageFile, pstfile=psfCandidatesFile, \ opstfile=opstFile, groupfile=groupFile, inter='no', verify='no', update='no', verbose='yes') print 'PSF took ',int(time.time() - startTime), ' seconds.' # see the psf iraf.noao.digiphot.daophot.seepsf(psfimage=psfImageFile, \ image=psfImageFile.replace('psf','seepsf'), dimension=100) # NSTAR & SUBSTAR: nstarFile = psfCandidatesFile.replace('pst','nst') rejFile = psfCandidatesFile.replace('pst','nrj') iraf.noao.digiphot.daophot.nstar(image=fileName, groupfile=groupFile, psfimage=psfImageFile, \ nstarfile=nstarFile, rejfile=rejFile, verify='no', update='no', verbose='yes') subImageFile = psfImageFile.replace('psf','psfsub') iraf.noao.digiphot.daophot.substar(image=fileName, photfile=nstarFile, \ exfile='', psfimage=psfImageFile, subimage=subImageFile, verify='no', update='no', verbose='no') # finally ALLSTAR: allStarFile = daofindOutputFile.replace('coo','als') allStarRejFile = daofindOutputFile.replace('coo','arj') startTime = time.time() iraf.noao.digiphot.daophot.allstar(image=fileName, photfile=photOutputFile, \ psfimage=psfImageFile, allstarfile=allStarFile, rejfile=allStarRejFile, \ subimage=subImageFile.replace('pntsrc_iso','daofind'), verify='no', update='no', verbose='no') print 'ALLSTAR took ',int(time.time() - startTime), ' seconds.' # pdump the attributes into a temp file for wcsctran: iraf.noao.digiphot.ptools.pdump(infile=allStarFile, fields='XCENTER,YCENTER,MAG,MERR,SHARPNESS,CHI', \ expr='yes', headers='no', parameters='no', Stdout=allStarFile.replace('als','pdp')) # apply WCS to the XY coords to get RA & Dec: iraf.images.imcoords.wcsctran(input=allStarFile.replace('als','pdp'), \ output=allStarFile.replace('als','tmp'), image=fileName, \ inwcs='logical', outwcs='world', columns='1 2', formats='%11.7f %11.7f') # write in headers, including tile number, lo/hi and filter/epoch in column name !!! hilo = 'hi' if fileName.endswith('_lo.fits'): hilo='lo' suffix = '_' + str(tileFile[5]) + '_' + str(tileFile[8]) + '_' + hilo header = '#' + \ ' RA' + suffix + \ ' Dec' + suffix + \ ' psfMag' + suffix + \ ' psfMagErr' + suffix + \ ' sharp' + suffix + \ ' chi2' + suffix + '\n' f = open(allStarFile.replace('als','tmp'), mode = 'r') g = open(allStarFile.replace('als','rdc'), mode = 'w') f.readline() f.readline() f.readline() g.write(header) line = f.readline() while line != '': g.write(line) line = f.readline() f.close() g.close() # Test one all the way through: TODO comment out !!! #exit('Successfuly completed first run.') exit('Successfully completed.')