; NAME: ; autophot.pro ; ; PURPOSE: ; Fully-automated photometry of CCD images ; ; CALLING SEQUENCE: ; autophot, imagefile, [objfile], objfile=objfile, calfile=calfile, ... ; [ see below for more] ; ; INPUTS: ; imagefile - Image file (.fits) ; objfile - File with list of object coordinates for photometry ; calfile - File with list of catalog magnitudes for reference stars ; filter - Image filter, if not specified in header ; filtkey - Filter keyword in header (default FILTER) ; radius - Aperture radius in arcseconds ; aprad - Aperture radius in pixels (?) ; starrad - Aperture radius in arcseconds for calibration stars only (for measuring flux of large science objects) ; skyradius - Sky annulus [inner,outer] radii in arcsec ; fullradius - Radius enclosing "100%" of the flux, if zeropoint given ; pos - An RA/dec pair ; raobj - Array of RA's for objects for photometry ; decobj - Array of dec's for objects for photometry ; exptime - Exposure time, if not specified in header ; expkey - Exposure time keyword (defaults EXPTIME,EXPOSURE,ELAPTIME) ; magcolnum - Specify column of magnitude values in calfile explicitly ; unccolnum - Specify column of magnitude errors in calfile explicitly ; detsigma - Minimum sigma above background to label as a detection ; limsigma - When not detected, quote limit as this many sigma ; fix - Do not recenter coordinates; appropriate for faint objects [should change this to default] ; starfix - Do not recenter stars ; zeropt - Magnitude zeropoint of image (for r=fullradius aperture) ; zairmass - does nothing ; zairterm - does nothing ; colorfilter - Check for color terms ; cboxsize - Diameter (in pixels) of reference star centering box ; printkey - Prints out this additional keyword with photometry ; fluxout - Print fluxes instead of magnitudes ; quiet - Suppress text output to console ; verbose - Print out lots of details to console ; sdig - n significant digits for mags ; calplotfile - print the calibration plot to a file (.eps) ; imcutoutfile - print the image cutout to a file (.png) ; outdata - structure containing output data ; cutmag - automatically reject calibration stars outside the two values. They still appear on plots. ; ; COMMENTS: ; Fast, automatic aperture photometry on an image ; given a list of objects for which magnitudes are desired and (optionally) ; another list of stars for which magnitudes are known. If a reference ; catalog is not available, one can be pulled over automatically from the web ; or a fixed zeropoint (with automatic aperture correction) can be used. ; ; This is a convenience routine for doing aperture photometry on an image, usually of a single object or ; small number of objects. In basic mode the program is pointed to an image, object coordinates, and a table ; of calibrated secondary standards with magnitudes. It performs photometry on the whole list and uses ; the reference-instrumental magnitude offset to calibrate the image. Basic outlier rejection is included. ; Also attempts to estimate the image zeropoint at "infinite" aperture: this does not enter into the ; normal calculation in basic mode. Can also be passed an assumed zeropoint at "infinite" aperture ; and a table of (uncalibrated) standard star positions. In this mode the aperture correction IS used to ; measure the object magnitude. If an aperture-specific zeropoint is known (stable PSF), currently this ; is specified by setting fullrad=rad. ; ; Still in testing is the ability to measure fluxes of large, bright objects accurately. ; ; EXAMPLES: ; autophot, 'image.fits', pos=['21:21:21','+00:30:11.3'], rad=1.2 ; autophot, 'image.fits', obj='objposfile.dat', cal='calfile.dat', rad=1.1, /fix ; ; ; ; ; satval doesn't actually seem to work. ; try autophot, 'azfpb120126_0124r.fits', 'grb120119a.blue.sdss.direct.cal', rad=3 ; this does a bad aperture correction; look into it pro autophot, image, objfile2, v3, calfile=calfile, objfile=objfile, radiusarcsec=radiusarcsec, skyradiusarcsec=skyradiusarcsec, fullradiusarcsec=fullradiusarcsec, nameobj=nameobj, posobj=posobj, raobj=raobj, decobj=decobj, aprad=aprad, filterin=filterin, filtkey=filtkey, exptime=exptime, expkey=expkey, magcolnum=magcolnum, unccolnum=unccolnum, quiet=quiet, silent=silent, detsigma=detsigma, limsigma=limsigma, fixedpos=fixedpos, sigma=sigma, cutmag=cutmag, minval=minval, satval=satval, invgain=invgain, zeropt=zeropt, zairmass=zairmass, zairterm=zairterm, verbose=verbose, colorfilter=colorfilter, cboxsize=cboxsize, printkey=printkey, fluxout=fluxout, sdig=sdig, plotcalfile=plotcalfile, imcutoutfile=imcutoutfile, imcutoutpix=imcutoutpix, nocutout=nocutout, imcutoutarcsec=imcutoutarcsec, outdata=outdata, noscreen=noscreen, starrad=starrad, savestarphot=savestarphot, colorterm=colorterm, table=table, starfix=starfix ; Command-line photometry of an object of interest in an image relative to a list of calibration stars. ; Check syntax and existence of files if n_elements(image) eq 0 then begin print, 'syntax: autophot, [image.fits], [calfile.cal], [targetcoordfile.pos], radius=[aperture in arcsec], ...' endif ; Reshuffle variables if old syntax used if n_elements(v3) gt 0 then begin calfile = objfile2 objfile = v3 endif else begin if n_elements(objfile2) gt 0 then objfile = objfile2 endelse ; Did user do something weird and provide arrays? if n_elements(image) gt 1 then image = image[0] if n_elements(calfile) gt 1 then calfile = calfile[0] if n_elements(objfile) gt 1 then objfile = objfile[0] ; Check that files exist if they are specified if file_test(image) eq 0 then print, 'File not found: ', image if n_elements(calfile) gt 0 then if file_test(calfile) eq 0 then print, 'File not found: ', calfile if n_elements(objfile) gt 0 then if file_test(objfile) eq 0 then print, 'File not found: ', objfile if file_test(image) eq 0 then return if n_elements(calfile) gt 0 then if file_test(calfile) eq 0 then return if n_elements(objfile) gt 0 then if file_test(objfile) eq 0 then return ; Set some defaults / move variable names if n_elements(filterin) gt 0 then filter = filterin quiet = keyword_set(quiet) if n_elements(cboxsize) eq 0 then cboxsize = 15 ; default star centering box (full)size if n_elements(posobj) eq 2 then begin raobj = posobj[0] decobj = posobj[1] endif if keyword_set(silent) then quiet = 1 ; Set some other variables from the (base) header that we need early on h = headfits(image) if n_elements(exptime) eq 0 then begin if n_elements(expkey) gt 0 then begin exptime = sxpar(h, expkey, count=ct) endif else begin exptime = sxpar(h,'EXPTIME', count=ct) ; appropriate for individual LRIS exposures if ct eq 0 then exptime = sxpar(h,'EXPOSURE', count=ct) ; appropriate for stacked exposures if ct eq 0 then exptime = sxpar(h,'ELAPTIME', count=ct) if exptime eq 0 then begin print, 'Specified exposure time is zero!' return endif endelse if ct eq 0 then begin print, 'Could not determine exposure time; setting to 1s' exptime = 1 endif endif airmass = sxpar(h,'AIRMASS', count=ct) if n_elements(satval) eq 0 then begin satval = sxpar(h,'SATURATE', count=ct) if ct eq 0 then satval = !values.f_infinity ;60000. endif if n_elements(invgain) eq 0 then invgain = !values.f_infinity ; assume no photon counting error (source noise) unless specified. if n_elements(filter) eq 0 and n_elements(magcolnum) eq 0 then begin if keyword_set(filtkey) then begin filter = clip(sxpar(h, filtkey, count=ct)) endif else begin filter = clip(sxpar(h, 'FILTER', count=ct)) if ct eq 0 then filter = clip(sxpar(h, 'FILTER1', count=ct)) if ct eq 0 then filter = clip(sxpar(h, 'FILTER2', count=ct)) if ct eq 0 then filter = clip(sxpar(h, 'FILTNAM', count=ct)) endelse if ct eq 0 then begin print, 'Cannot identify filter! Specify filter= or filtkey=.' return endif endif filter = stdfilter(filter) ; standardize filter name if n_elements(sigma) gt 0 then begin if n_elements(detsigma) eq 0 then detsigma = sigma if n_elements(limsigma) eq 0 then limsigma = sigma endif if n_elements(detsigma) eq 0 then detsigma = 1.5 if n_elements(limsigma) eq 0 then limsigma = 2 ; Read in object coordinate file, if specified if (n_elements(raobj) eq 0 or n_elements(decobj) eq 0) and n_elements(objfile) gt 0 then begin nl = countlines(objfile) if nl eq 0 then begin print, 'Object coordinate file is blank?' return endif nameobj = strarr(nl) raobj = dblarr(nl) decobj = dblarr(nl) openr, 2, objfile inline = '' nobj = 0 for l = 0, nl-1 do begin readf, 2, inline inline = clip(inline) if strlen(inline) eq 0 then continue if strmid(inline,0,1) eq '#' then continue indata = strsplit(inline,/extract) if n_elements(indata) ge 2 then begin if n_elements(indata) eq 2 then begin nameobj[nobj] = clip(nobj) rastr = indata[0] decstr = indata[1] endif if n_elements(indata) ge 3 then begin nameobj[nobj] = indata[0] rastr = indata[1] decstr = indata[2] endif if strpos(rastr,':') gt 0 then raobj[nobj] = 15*ten(rastr) else raobj[nobj] = float(rastr) if strpos(decstr,':') gt 0 then decobj[nobj] = ten(decstr) else decobj[nobj] = float(decstr) nobj += 1 endif endfor close, 2 nameobj = nameobj[0:nobj-1] raobj = raobj[0:nobj-1] decobj = decobj[0:nobj-1] endif else begin if n_elements(raobj) ne n_elements(decobj) then print, 'Warning - ra and dec arrays have different lengths.' if n_elements(nameobj) lt n_elements(raobj) then begin if n_elements(nameobj) eq 0 then nameobj = strarr(n_elements(raobj)) $ else nameobj = [nameobj, strarr(n_elements(raobj) - n_elements(nameobj))] endif if size(raobj,/tn) eq 'STRING' or size(decobj,/tn) eq 'STRING' then begin inraobj = raobj indecobj = decobj raobj = double(n_elements(inraobj)) decobj = double(n_elements(indecobj)) for i = 0, n_elements(raobj)-1 do begin if strpos(inraobj[i],':') gt 0 then raobj[i] = ten(inraobj[i])*15. if strpos(indecobj[i],':') gt 0 then decobj[i] = ten(indecobj[i]) endfor endif endelse ; If no calibration file given, download one from the internet if n_elements(calfile) eq 0 and n_elements(zeropt) eq 0 then calfile = 'auto' if n_elements(calfile) gt 0 then if calfile eq 'auto' then begin filt = strupcase(strmid(filter,0,1)) if total(filter eq ['u','U','B','g','V','r','R','i','I','z','Z','Y','J','H','K','Kp','Ks']) eq 0 then begin print, 'Nonstandard filter (', filter, '): specify zeropoint, star catalog, or standard filter.' return endif tempcat = removepath(image)+'.cat' if file_test(tempcat) eq 0 then begin if filt eq 'Y' or filt eq 'H' or filt eq 'J' or filt eq 'K' then begin printcatalog, image, '2mass', outfile=tempcat, /quiet, /unc, count=count autocat = '2MASS' endif else begin printcatalog, image, 'sdss', outfile=tempcat, /quiet, /unc, count=count autocat = 'SDSS' if count eq 0 then begin if total(filter eq ['R','I','B']) eq 0 then begin print, 'Filter (',filter,') not in USNO; specify filter=R,I,B or star catalog' return endif autocat = 'USNO' printcatalog, image, 'usno', outfile=tempcat, /quiet, /unc, count=count, magthresh=18 endif endelse if count eq 0 then begin print, 'No ',autocat,' reference stars found in field. ' ;; (FOV '+fpr()'x'+fpr()+)' return endif endif else begin autocat = tempcat endelse print, 'No calibration file given; using '+autocat calfile = tempcat endif ; Read in calibration star file, if specified ns = 0 if n_elements(calfile) gt 0 then begin if n_elements(zeropt) eq 0 then cat = readcalibration(calfile, magcol=magcolnum, unccolnum=unccolnum, filter=filter, f2=colorfilter) $ else cat = readcalibration(calfile, magcol=magcolnum, unccolnum=unccolnum, filter='') ns = n_elements(cat) if size(cat[0], /tn) ne 'STRUCT' then begin return endif colorsign = 1. if keyword_set(colorterm) then cat.mag += colorterm*colorsign*(cat.mag-cat.mag2) endif if (n_elements(raobj) eq 0 or n_elements(decobj) eq 0) and n_elements(objfile) eq 0 then begin print, 'No object coordinate file specified - calculating zeropoint only.' endif ; Turn off math warnings except = !except !except = 0 ; will stick if returns early ; Read in the FITS file and get information from the header data = mrdfits(image,0,h0, /silent) if n_elements(data) eq 1 then begin data = mrdfits(image,1,hwcs, /silent) endif else begin hwcs = h0 endelse bzero = sxpar(hwcs, 'BZERO', count=ctbzero) bscale = sxpar(hwcs, 'BSCALE', count=ctbscale) if ctbscale gt 0 then data *= sxpar(hwcs, 'BSCALE') if ctbzero gt 0 then data += sxpar(hwcs, 'BZERO') sxaddpar, hwcs,'BZERO', 0 sxaddpar, hwcs,'BSCALE', 0 ; Check for header astrometry extast, hwcs, astr, noparams, alt=alt if noparams lt 0 then begin print, 'No header astrometry for '+image+'; cannot do WCS-based photometry.' return endif ; some weird Pan-STARRS thing I don't fully understand if sxpar(hwcs,'PC001001') eq -1 then sxaddpar,hwcs, 'CDELT1',-sxpar(hwcs,'CDELT1') if sxpar(hwcs,'PC002002') eq -1 then sxaddpar,hwcs, 'CDELT2',-sxpar(hwcs,'CDELT2') pixscale = getpixelscale(hwcs) cd11 = sxpar(hwcs,'CD1_1') cd12 = sxpar(hwcs,'CD1_2') cd21 = sxpar(hwcs,'CD2_1') cd22 = sxpar(hwcs,'CD2_2') if cd21 ne 0 or cd11 ne 0 then parad = atan2(cd21, -cd11) $ ; need a differential pixel scale adjustment, in principle else parad = 0. padeg = parad * 180/!pi parity = -((cd11 * cd22) < 0 or (cd12 * cd21 > 0)*2)+1 ; either -1 or 1 if n_elements(minval) eq 0 then begin minval = (median(data)-1000) < (-1000) endif if n_elements(zeropt) gt 0 then begin if n_elements(radiusarcsec) gt 0 and n_elements(fullradiusarcsec) gt 0 then $ if radiusarcsec eq fullradiusarcsec then $ offset = 25 - zeropt ; offset is the instrumental(if 25)-real zeropoint endif ; Determine some data parameters and create pixel selection arrays nx = (size(data)) [1] ny = (size(data)) [2] datax = intarr(nx,ny) for x = 0, nx-1 do datax[x,*] = x datay = intarr(nx,ny) for y = 0, ny-1 do datay[*,y] = y wbad = where(finite(data) eq 0, ct) if ct gt 0 then data[wbad] = -!values.f_infinity ; replace nans with -Inf since can't NaN-check and limit aper simultaneously cboxr = fix(cboxsize-1)/2 if ns gt 0 then begin ; Most of the code executes only if secondary standards are specified. ; this is equivalent to there being a calfile. ; Set up some arrays to hold reference star data errmsg = strarr(ns) apmag = fltarr(ns) + !values.f_nan apunc = fltarr(ns) + !values.f_nan fullmag = fltarr(ns) + !values.f_nan fullunc = fltarr(ns) + !values.f_nan ns_field = 0 ax = fltarr(ns) + !values.f_nan ; from the fits astrometry ay = fltarr(ns) + !values.f_nan cx = fltarr(ns) + !values.f_nan ; (revised) locations of the reference stars, pix coordinates cy = fltarr(ns) + !values.f_nan kern = [[0.1,0.1,0.1],[0.1,0.2,0.1],[0.1,0.1,0.1]] kern = psf_gaussian(npix=5, fwhm=3, ndimen=2, /normalize) smdata = convol(data,kern) ; lightly smooth for centroiding ; Reference star "centroiding" (max pixel value...) for s = 0, ns-1 do begin adxy, hwcs, cat[s].ra, cat[s].dec, x, y ; this is probably in FITS xy coordinates (starts with 1?) ax[s] = x ay[s] = y if keyword_set(starfix) then begin cx[s] = ax[s] cy[s] = ay[s] continue endif rx = fix(round(x)) ry = fix(round(y)) if x lt cboxr or y lt cboxr or x gt nx-cboxr or y gt ny-cboxr or finite(x) eq 0 or finite(y) eq 0 then begin errmsg[s] = 'Off chip' continue endif cutout = data[(rx-cboxr)>0:(rx+cboxr)<(nx-1),(ry-cboxr)>0:(ry+cboxr)<(ny-1)] cutoutsm=smdata[(rx-cboxr)>0:(rx+cboxr)<(nx-1),(ry-cboxr)>0:(ry+cboxr)<(ny-1)] ; requires sky subtraction for proper centroiding cutoutx = datax[(rx-cboxr)>0:(rx+cboxr)<(nx-1),(ry-cboxr)>0:(ry+cboxr)<(ny-1)] cutouty = datay[(rx-cboxr)>0:(rx+cboxr)<(nx-1),(ry-cboxr)>0:(ry+cboxr)<(ny-1)] ;smoothcutout = convol(cutout,fltarr(3,3)+1./9) ;resampcutout = congrid(cutout, 200,200) ;window, 1, xsize=(size(resampcutout)) [1], ysize=(size(resampcutout)) [2] ;tv, cutout if max(cutout) eq min(cutout) then begin errmsg[s] = 'Blank region' continue endif mp = (where(cutoutsm eq max(cutoutsm), ct)) [0] if ct eq 0 then begin apmag[s] = 99. errmsg[s] = 'No valid pix in aperture' continue endif maxx = cutoutx[mp] maxy = cutouty[mp] cutout2 = data[(maxx-3)>0:(maxx+3)<(nx-1),(maxy-3)>0:(maxy+3)<(ny-1)] s2 = size(cutout2) xcen = total( total(cutout2, 2) * indgen(s2[1]) ) / total(cutout2) ycen = total( total(cutout2, 1) * indgen(s2[2]) ) / total(cutout2) xcen += (maxx-3)>0 ycen += (maxy-3)>0 cx[s] = xcen cy[s] = ycen ;ok = where(finite(cutout)) ;if total(cutout[ok]) eq 0 then continue ;cx[s] = total(cutout[ok]*cutoutx[ok],0) / total(cutout[ok],0) ;cy[s] = total(cutout[ok]*cutouty[ok],1) / total(cutout[ok],1) ns_field += 1 ;cntrd, data, cx[s], cy[s], xcen, ycen, 15 ; get the exact position via centroid (instead of max pixel)... ; cx[s] = xcen ; cy[s] = ycen endfor defaultfull = 0 defaultsky = 0 if n_elements(fullradiusarcsec) eq 0 then begin if n_elements(radiusarcsec) gt 0 then fullradiusarcsec = (1.0 > (radiusarcsec*4)) < 7.0 $ else fullradiusarcsec = 15. defaultfull = 1 endif fullradius = fullradiusarcsec/pixscale if n_elements(skyradiusarcsec) eq 0 then begin if n_elements(radiusarcsec) gt 0 then begin skyradiusarcsec=[5.4 > radiusarcsec, (5.4 > radiusarcsec)+1.35] if skyradiusarcsec[0] lt fullradiusarcsec then skyradiusarcsec += (fullradiusarcsec-skyradiusarcsec[0]) if skyradiusarcsec[0] lt radiusarcsec then skyradiusarcsec += (radiusarcsec-skyradiusarcsec[0]) endif else begin skyradiusarcsec=[6.0,10.0] defaultsky = 1 endelse endif skyradius = skyradiusarcsec/pixscale ; Aperture-dependent zeropoint determination if n_elements(radiusarcsec) eq 0 then begin ; Measure optimum aperture using reference stars tryrad = 1.0 + 0.1*findgen((fullradius-1.0)/0.1) bestrads = !values.f_nan + fltarr(ns) for s = 0, ns-1 do begin if finite(cx[s]) eq 0 or finite(cy[s]) eq 0 then continue aper2, data, cx[s], cy[s], testradflux, testraderr, testsky, testskyerr, invgain, tryrad, skyradius, [-1000,satval], /silent, /flux w = where(testraderr/testradflux eq min(testraderr/testradflux), ct) if ct gt 0 then bestrads[s] = tryrad[w] endfor radius = median(bestrads) radiusarcsec = radius*pixscale print, 'Calculated optimum aperture: ', fpr(radius,2.2)+'px = ', fpr(radiusarcsec,2.3)+'"' endif radius = radiusarcsec/pixscale if n_elements(starrad) eq 0 then begin starradius = radius endif else begin starradius = starrad/pixscale endelse starradiusarcsec = starradius*pixscale ; if n_elements(starradius) eq 0 and n_elements(radius) gt 0 then if fullradius lt radius then starradius = fullradius*0.5 ; strange conditonal? ; Reference star photometry for s = 0, ns-1 do begin if finite(cx[s]) eq 0 or finite(cy[s]) eq 0 then continue aper2, data, cx[s], cy[s], apfluxarr, apfluxerrarr, sky, skyerr, invgain, [starradius, fullradius], skyradius, [-1000,satval], /silent, /flux ;, /nan signoise = apfluxarr/apfluxerrarr apfluxarr /= exptime apfluxerrarr /= exptime sky /= exptime skyerr /= exptime apmagarr = 25 - 2.5*alog10(apfluxarr) apmagerrarr = -(25 - 2.5*alog10(apfluxarr+apfluxerrarr)) + apmagarr apmag[s] = (apmagarr[0,*]) [*] apunc[s] = (apmagerrarr[0,*]) [*] fullmag[s] = (apmagarr[1,*]) [*] fullunc[s] = (apmagerrarr[1,*]) [*] if apmag[s] eq 99. then errmsg[s] = 'Aper error.' endfor used = indgen(ns) isvalid = apmag gt 0 and apmag lt 50 and finite(apmag) and finite(cat.mag) wherevalid = where(isvalid, nvalid) if nvalid eq 0 then begin print, 'Found no good calibration stars. Check field overlap or saturation.' return endif catuncthresh = 0.20 > min(cat[wherevalid].unc)*1.5 ; formerly 0.20 apuncthresh = 0.20 > min(apunc[wherevalid])*1.5 if catuncthresh gt 0.20 or apuncthresh gt 0.20 then begin print, 'WARNING: No overlapping stars that are well-detected in both image and catalog.' print, 'Calibration may be detection-biased or subject to other systematics.' endif iswelldetected = cat.unc lt catuncthresh and apunc lt apuncthresh nused = 0 if n_elements(zeropt) eq 0 then begin ;w = where(apmag gt 0 and apmag lt 50 and finite(apmag) and finite(cat.mag) and cat.unc lt 0.25 and apunc lt 0.30, ct) w = where(isvalid and iswelldetected, ct) ; ct should always be >0 here sdiff = fltarr(ns) + !values.f_nan sdiff[w] = apmag[w]-cat[w].mag isvalid = intarr(ns) ; set things not well detected to no longer valid, I guess, not sure why I do this isvalid[w] = 1 nvalid = total(isvalid) if n_elements(cutmag) eq 0 then cutmag = !values.f_infinity*[-1,1] if n_elements(cutmag) eq 1 then cutmag = [-!values.f_infinity,cutmag] good = where(finite(sdiff) and cat.mag gt cutmag[0] and cat.mag lt cutmag[1]) diff = sdiff[good] used = used[good] isrejected = intarr(ns) isrejected[good] = 1 ; Iterate a few times to remove outliers offset = median(diff, /even) if n_elements(diff) gt 1 then rms = stdev(diff-offset) else rms = 0. absdev = median(abs(diff-offset),/even) if n_elements(diff) gt 2 then begin itersig = [6,5,4,4,4,4] for iter = 0, 5 do begin ;print, fpr(diff-offset,3.1) goodthresh = itersig[iter]*(absdev > 0.05) good = where(abs(diff-offset) lt goodthresh) diff = diff[good] used = used[good] offset = median(diff, /even) rms = stdev(diff-offset) absdev = median(abs(diff-offset),/even) endfor endif isused = intarr(ns) isused[used] = 1 isrejected[used] = 0 absdiffminuserror = abs(sdiff-offset)-2*sqrt(cat.unc^2+apunc^2) caution = where(isused and (absdiffminuserror gt (3*(rms > 0.025) < 0.1)), ctc) plotcol = intarr(ns) ; 0 = skipped due to cutmag, off frame, etc. plotcol[used] = 1 ; white/black if ctc gt 0 then plotcol[caution] = 3 ; yellow if total(isrejected) gt 0 then plotcol[where(isrejected)] = 2 ; red displaycols = [rgb(160,160,160), rgb(128,255,128), rgb(200,0,0), rgb(255,255,0)] plotcolrgb = displaycols[plotcol] apmag = apmag - offset ; true mag = instmag - offset fullmag = fullmag - offset ; Print/plot calibration info for user w = where(isused, nused, complement=notused) if keyword_set(verbose) then begin tstr = clip('#RA',12)+' '+clip('dec',10)+' '+clip('RA_cat',12)+' '+clip('dec_cat',10)+' '+$ clip(filter+'_cal',5)+' '+clip(filter+'_err',5)+' ' if n_elements(colorfilter) gt 0 then tstr += clip(colorfilter+'_cal',5)+' '+clip(colorfilter+'_err',5)+' ' tstr += clip(filter+'_ap',5)+' '+clip(filter+'_aperr',7)+' used diff' print, tstr coutra = dblarr(ns) coutdec = dblarr(ns) for s = 0, ns-1 do begin xyad, h, cx[s], cy[s], a, d coutra[s] = a coutdec[s] = d colorstr = '' colorvstr = '' if n_elements(colorfilter) gt 0 then colorstr = fprn(cat[s].mag2,2.2)+' '+fprn(cat[s].unc2,1.2)+' ' if n_elements(colorfilter) gt 0 then colorvstr = fpr(cat[s].mag-cat[s].mag2,2.2)+' ' ; fpr(cat[s].ra,3.6),' ',fpr(cat[s].dec,3.6),' ',fpr(cx[s],5.2),' ',fpr(cy[s],5.2), if isvalid[s] and iswelldetected[s] then $ print, fpr(a,3.7),' ',fpr(d,3.7),' ',fpr(cat[s].ra,3.7),' ',$ fpr(cat[s].dec,3.7),' ',fprn(cat[s].mag,2.2),' ',fprn(cat[s].unc,1.2)+' '+$ ' '+colorstr+fprn(apmag[s],2.2),' ',fprn(apunc[s],1.2),' ',clip(isused[s]), $ ' '+colorvstr+fpr(apmag[s]-cat[s].mag,2.2) endfor if n_elements(colorfilter) and keyword_set(noscreen) eq 0 then begin window, 1, xsize=500, ysize=350, title='Autophot: calibration plot' wp = where(cat.mag ge 0) ;where(cat.mag lt max(cat[w].mag), ct) megaplot, cat[wp].mag-cat[wp].mag2, apmag[wp]-cat[wp].mag, $ xerr=0.01+sqrt(cat[wp].unc^2 + cat[wp].unc2^2), $ yerr=0.01+sqrt(cat[wp].unc^2 + apunc[wp]^2), $ xtitle=clip(filter+'_cal')+' - '+clip(colorfilter+'_cal'), $ ytitle=clip(filter+'_meas')+' - '+clip(filter+'_cal'), color=plotcolrgb[wp], errcolor=plotcolrgb[wp], yrange=[-0.5,0.5] oplot, [-10,10], [0,0], color=rgb(128,128,128) endif endif else begin n_offfield = n_elements(cat) - ns_field if n_offfield gt 0 then offstr = ' ('+clip(n_offfield)+ ' off field)' if quiet eq 0 then print, clip(fix(total(isvalid)))+'/'+clip(ns_field)+' valid calibration stars ' endelse ; plot ranges xymin = min([cat.mag-(cat.unc < 0.2), apmag[w]-(apunc[w] < 0.2)])-0.2 xymax = max([cat.mag+(cat.unc < 0.2), apmag[w]])+0.1 if xymax-xymin lt 1.0 then begin spr = 2.0-(xymax-xymin) xymin -= spr/2. xymax += spr/2. endif np = 1 + (n_elements(plotcalfile) gt 0) for p = 1 + keyword_set(noscreen), np do begin if p eq 1 then window, 0, xsize=500, ysize=500, title='Autophot: calibration' if p eq 2 then begin psopen, plotcalfile, xsize=4, ysize=4, /inches, /encaps device, /color colors = transpose([[0,0,0],$ [150,150,150],$ [0,0,255],$ [255,255,255]]) tvlct, colors endif !p.position = [0.09, 0.08, 0.98, 0.98] if p eq 2 then !p.position = [0.14, 0.13, 0.98, 0.98] if p eq 1 then color=plotcolrgb if p eq 2 then color=1+isused plot, [cat.mag], [apmag], /nodata, /ynozero, xrange=[xymin,xymax],yrange=[xymin,xymax], /xsty, /ysty, xtitle='catalog', ytitle='measured' for i = 0, n_elements(cat)-1 do oploterror1, cat[i].mag, apmag[i], cat[i].unc+1e-3, apunc[i]+1e-3, errcolor=color[i] ; oploterror is really slow due to the cg call. ;megaplot, [cat.mag], [apmag], xerr=cat.unc+1e-3, yerr=apunc+1e-3, color=color, errcolor=color, /noerase ; uc = unique(color) ; for c = 0, n_elements(uc)-1 do begin ; wc = where(color eq uc[c]) ; oploterror, cat[wc].mag, apmag[wc], cat[wc].unc+1e-3, apunc[wc]+1e-3, color=uc[c] ; endfor oplot, [0,30], [0,30] if p eq 2 then psclose if p eq 1 and n_elements(plotcalfile) gt 0 then begin if strpos(strlowcase(plotcalfile),'.png') gt 0 then begin WRITE_PNG, plotcalfile, TVRD(/TRUE) break endif endif endfor expectedrms = sqrt(total(cat[used].unc^2 + apunc[used]^2)/n_elements(used)) ;print, 'Mag shift ', offset if quiet eq 0 then print, 'observed-catalog RMS = ', fpr(rms,2.4), ' (', fpr(expectedrms,2.4), ' expected)' colorscatter = sqrt(rms^2 - expectedrms^2) offsetunc = 1./sqrt(total(1./(cat[used].unc^2 + apunc[used]^2))) systematicunc = sqrt(offsetunc^2 + (colorscatter > 0.)^2) if quiet eq 0 then $ if rms gt expectedrms then print, 'Color scatter ', sqrt(rms^2 - expectedrms^2) ; Show cutouts of calibration stars (up to 100) if keyword_set(noscreen) eq 0 then begin nsx = (ceil(sqrt(ns)) < 9) > 2 nsy = (ceil(sqrt(ns)) < 9) > 2 ; number per side swr = floor(600/nsx/2) sws = swr*2+1 xsize = sws*nsx ysize = sws*nsy window, 4, xsize=xsize, ysize=ysize, title='Autophot: star apertures', xpos=400 ; object window is 250 polyfill, [0,0,1,1],[0,1,1,0], color=rgb(100,0,0), /norm for s = 0, (ns < 81)-1 do begin if finite(cx[s]) eq 0 or finite(cy[s]) eq 0 then continue x = cx[s] y = cy[s] sx = s mod nsx sy = s / nsy z = round(17.2*pixscale) rotdata = rot(data, -padeg, z, x+0.5, y+0.5, missing=!values.f_nan) outs = size(rotdata) ; this rotates the image again every time... rcx = outs[1]/2 rcy = outs[2]/2 ; need to pad with NaNs or something near the edge, otherwise off-center. rotdatas = rotdata[(rcx-swr)>0:(rcx+swr)<(nx-1), (rcy-swr)>0:(rcy+swr)<(ny-1)] ; zoom in on object ok = where(finite(rotdatas), ctok) rots = size(rotdatas) ;rotdata[where(finite(rotdata) eq 0)] = 0 if x lt 0 or y lt 0 or x gt nx-1 or y gt nx-1 then continue maxinap = max(data[(x-radius)>0:(x+radius)<(nx-1),(y-radius)>0:(y+radius)<(ny-1)]) maxsigma = 5. sigclipstats, rotdatas, sigmahi=3., sigmalo=3., med=datamed, stdev=datastd if maxinap gt 5*datastd then maxsigma = sqrt((maxinap/datastd) * (5.)) ; geometricmean minv = datamed-(maxsigma/5.)*datastd ; sky is always 255/(5+1)=51 maxv = maxinap ;print, median(rotdatas), minv, maxv rotdatasc = 255*sqrt((0 > (rotdatas-minv)/(maxv-minv)) < 1) tv, rotdatasc, s wait, 0.001 x0n = (sx*1.)/nsx x1n = (sx+1.)/nsx-1./600 y0n = (nsy-1-sy*1.)/nsy y1n = (nsy-1-sy+1.)/nsy-1./600 xcen = (sx*1.+0.5)/nsx ycen = (nsy-1-sy*1.+0.5)/nsy xyouts, x0n+2./600, y0n+3/600., clipzero(fpr(cat[s].mag,2.3))+' '+string(apmag[s]-cat[s].mag,'(F+5.2)'), charsize=1, /norm, color=plotcolrgb[s] plots, xcen+[-2,2]/600., ycen+[-2,2]/600., /norm, color=rgb(255,0,0) plots, xcen+[2,-2]/600., ycen+[-2,2]/600., /norm, color=rgb(255,0,0) plots, [x0n,x1n,x1n,x0n,x0n], [y0n,y0n,y1n,y1n,y0n], /norm, color=plotcolrgb[s] apcolor = rgb(255,125,125) if isused[s] eq 0 then apcolor=rgb(225,155,155) tvcircle, radius*z, xcen*xsize, ycen*ysize, apcolor, /device endfor endif ; star apertures endif else begin ; User-specified zeropoint. ; (But, also calibration table given. If no calibration stars this is still skipped.) rms = 0 isused = isvalid endelse ; Aperture correction and "true" zeropoint determination wc = where(isused and finite(fullmag) and fullunc lt 0.15, ct) if ct gt 0 then begin ; apcorr = instmag - fullmag apdiffs = apmag[wc] - fullmag[wc] apdiffuncs = sqrt(fullunc[wc]^2 + 0.02^2) ; not actually sure... if n_elements(apdiffs) ge 5 then begin gooddif = (sort(apdiffs)) [1:n_elements(apdiffs)-2] apdiffs = apdiffs[gooddif] ; remove extreme values apdiffuncs = apdiffuncs[gooddif] endif apcorr = total(apdiffs * (1/apdiffuncs^2)) / total(1/apdiffuncs^2) apcorrchisqdof = total(((apdiffs-apcorr)/apdiffuncs)^2) / (n_elements(wc)-1. > 1) apcorrunc = sqrt(total(1 / total(1/apdiffuncs^2))) * (sqrt(apcorrchisqdof) > 1.0) endif else begin if total(finite(fullmag)) gt 0 then $ print, 'No sufficiently bright stars with which to perform aperture correction.' $ else $ if quiet eq 0 then print, 'Unable to calculate aperture correction. Check sky/aperture radii.' ; aper refuses to calculate mags if r_aper > r_skyouter apcorr = !values.f_nan apcorrunc = !values.f_nan endelse ; truemag = fullinstmag - (zeropt-25) ; truemag = instmag - offset ; zeropt = 25 + fullinstmag + offset - instmag ; = 25 + offset - apcorr if n_elements(zeropt) eq 0 then begin zeropt_measured = 25 - (offset - apcorr) if starradius ne radius and radius ge fullradius then offset -= apcorr ; CRUDE CONDITIONAL endif else begin offset = 25 + apcorr - zeropt apmag = apmag - offset endelse if quiet eq 0 then begin if apcorrunc gt 0.06 then warn = ' large uncertainty' else warn = '' if starradius ne radius then snote = '(for stars)' else snote = '' if starradius ne radius then snote2 = ' - applying no aperture correction to object' else snote2 = '' if finite(apcorr) and fullradius ne radius then print, 'Aperture correction '+snote+'= ' + fpr(apcorr,2.2) + ' +/- ' + fpr(apcorrunc,0.2) + ' ['+clip(fpr(starradiusarcsec,2.1))+'" -> '+clip(fpr(fullradiusarcsec,2.1))+'"]'+warn+snote2 endif if n_elements(zeropt_measured) gt 0 then begin zeropt_meas_unc = sqrt(apcorrunc^2 + systematicunc^2) if quiet eq 0 then $ if finite(apcorr) then print, 'Measured zeropoint (r='+clip(fpr(fullradiusarcsec,2.2))+'", secz='+fpr(airmass,1.2)+') = ', fpr(zeropt_measured,2.2), ' +/-', fpr(zeropt_meas_unc,2.2) endif ; end star catalog photometry endif else begin ; No star file provided: zeropoint is assumed to be for whatever aperture was given (no aperture correction) if n_elements(radiusarcsec) eq 0 then begin print, ' No radius specified.' return endif if n_elements(skyradiusarcsec) eq 0 then skyradiusarcsec=radiusarcsec*[2,3.5] radius = radiusarcsec/pixscale skyradius = skyradiusarcsec/pixscale endelse ; Object photometry ; nphot = 0 old nobj = n_elements(raobj) if nobj gt 0 then begin objboxr = 9 kernelsize = 4*objboxr+1 kernelmid = fix(kernelsize/2) initboxr = cboxr + kernelsize krad = radius / 2.35 kernel = fltarr(kernelsize,kernelsize) for x = 0, kernelsize-1 do begin for y = 0, kernelsize-1 do begin kernel[x,y] = gaussian(x, [1, kernelmid, krad]) * gaussian(y, [1, kernelmid, krad]) endfor endfor kernel /= total(kernel) ;writefits, 'kernel.fits', kernel ;objapmag = fltarr(nobj) ;objapunc = fltarr(nobj) outdata = replicate({objname:'', filter:'', limit:0, $ mag:!values.f_nan, magunc:!values.f_nan, sysunc:!values.f_nan, maglim:0., $ flux:!values.f_nan, fluxunc:!values.f_nan, fluxlim:0.0, $ ns_field:0, nvalid:0, nused:0, $ ra:!values.f_nan, dec:!values.f_nan, radiusarcsec:0.0, $ x:!values.f_nan, y:!values.f_nan, radius:0.0}, nobj) bout = 0 for b = 0, nobj-1 do begin nphot = b ; this is a change from before and may need to be reverted. adxy, hwcs, raobj[b], decobj[b], x, y ; This does seem to convert to IDL coordinates, not DS9 coordinates. ; Integer values are centered in the middle of the respective pixels. ; Confirmed in practice and in adxy documentation. if x lt 0 or x gt nx then continue if y lt 0 or y gt ny then continue if finite(x) eq 0 or finite(y) eq 0 then continue if x eq 1 and y eq 1 then continue if keyword_set(fixedpos) eq 0 then begin rx = fix(round(x)) ry = fix(round(y)) initcutout = data[(rx-initboxr)>0:(rx+initboxr)<(nx-1),(ry-initboxr)>0:(ry+initboxr)<(ny-1)] smoothinitcutout = convol(initcutout, kernel, /edge_truncate) ; cutout = initcutout[((rx-objboxr)>0) - ((rx-initboxr)>0):(rx+objboxr)<(nx-1) - ((rx-initboxr)>0), $ ; ((ry-objboxr)>0) - ((ry-initboxr)>0):(ry+objboxr)<(ny-1) - ((ry-initboxr)>0)] smoothcutout = smoothinitcutout[((rx-objboxr)>0) - ((rx-initboxr)>0):(rx+objboxr)<(nx-1) - ((rx-initboxr)>0), $ ((ry-objboxr)>0) - ((ry-initboxr)>0):(ry+objboxr)<(ny-1) - ((ry-initboxr)>0)] cutoutx = datax[(rx-objboxr)>0:(rx+objboxr)<(nx-1),(ry-objboxr)>0:(ry+objboxr)<(ny-1)] cutouty = datay[(rx-objboxr)>0:(rx+objboxr)<(nx-1),(ry-objboxr)>0:(ry+objboxr)<(ny-1)] mp = (where(smoothcutout eq max(smoothcutout))) [0] maxx = cutoutx[mp] maxy = cutouty[mp] cutout2 = data[(maxx-3)>0:(maxx+3)<(nx-1),(maxy-3)>0:(maxy+3)<(ny-1)] s2 = size(cutout2) xcen = total( total(cutout2, 2) * indgen(s2[1]) ) / total(cutout2) ; centroid ycen = total( total(cutout2, 1) * indgen(s2[2]) ) / total(cutout2) xcen += (maxx-3)>0 ycen += (maxy-3)>0 x = xcen y = ycen endif if n_elements(offset) eq 0 then begin print, 'Cannot determine aperture correction without a star catalog' return endif ; Display a window indicating the photometry aperture and annulus if bout eq 0 and keyword_set(nocutout) eq 0 and keyword_set(noscreen) eq 0 then begin ;print, 'PA = ', padeg ;print, 'parity = ', parity ;z = 20 < round(250/(2*skyradius[1])) > 1 ; zoom factor ;wr = 50/z > (skyradius[1]+5) < 300/z ; window "radius" in image pixels ;ws = ((wr*2)+1)*z ; window size in screen pixels if n_elements(imcutoutpix) eq 0 then ws = 250 else ws = imcutoutpix if n_elements(imcutoutarcsec) eq 0 then z = ws/(2*(skyradius[1]*1.01)) else z = ws/(imcutoutarcsec/pixscale) rotdata = rot(data, -padeg, z, x+0.5, y+0.5, missing=0) ; magnified, rotated, translated data (object now at image center: nx/2, ny/2) ; The +0.5 seems necessary due to pixel center vs pixel edge reasons... verified visually for various parities. ; The IDL/DS9 array conventions otherwise don't seem to matter. if parity eq -1 then rotdata = rotate(rotdata, 5) outs = size(rotdata) if ws gt outs[1]-1 then ws = outs[1]-1 if ws gt outs[2]-1 then ws = outs[2]-1 rotdatas = rotdata[outs[1]/2-ws/2:outs[1]/2+ws/2, outs[2]/2-ws/2:outs[2]/2+ws/2] ; zoom in on object rots = size(rotdatas) ;rotdata[where(finite(rotdata) eq 0)] = 0 sigclipstats, rotdatas, sigmahi=3., sigmalo=3., med=datamed, stdev=datastd maxinap = max(data[(x-radius)>0:(x+radius)<(nx-1),(y-radius)>0:(y+radius)<(ny-1)]) maxsigma = 5. if maxinap gt 5*datastd then maxsigma = sqrt((maxinap/datastd) * (5.)) ; geometricmean minv = datamed-(maxsigma/5.)*datastd ; sky is always 255/(5+1)=51 maxv = datamed+(maxsigma)*datastd ;minv = min(rotdata) if maxinap gt maxv then maxv = datamed+min([20*datastd, maxinap]) rotdatasc = 255*((0 > (rotdatas-minv)/(maxv-minv)) < 1) window, 2, xsize=ws, ysize=ws, title='Autophot: '+(nameobj[b] eq '' ? 'target' : clip(nameobj[b])) tv, rotdatasc tvcircle, skyradius[0]*z+(abs(radius-skyradius[0]) < 1 ? 1:0), rots[1]/2, rots[2]/2, rgb(255,255,0), /device tvcircle, skyradius[1]*z, rots[1]/2, rots[2]/2, rgb(255,255,0), /device tvcircle, radius*z, rots[1]/2, rots[2]/2, rgb(255,125,125), /device if n_elements(imcutoutfile) gt 0 then WRITE_PNG, imcutoutfile, TVRD(/TRUE) endif aper2, data, x, y, apfluxarr, apfluxerrarr, sky, skyerr, invgain, radius, skyradius, [minval,satval], /silent, /flux, /nan objapflux = apfluxarr[0]/exptime objapfluxunc = apfluxerrarr[0]/exptime signoise = objapflux/objapfluxunc sky = sky[0]/exptime skyerr = skyerr[0]/exptime objapmag = 25 - 2.5*alog10(objapflux) - offset objapmagunc = -(25 - 2.5*alog10(objapflux+objapfluxunc) - offset) + objapmag xyad, hwcs, x, y, a, d limflux = (objapflux > 0) + limsigma*objapfluxunc limmag = 25 - 2.5*alog10(limflux) - offset magform = 2.2 errform = 1.2 if n_elements(sdig) gt 0 then magform = 2 + sdig/10. if n_elements(sdig) gt 0 then errform = 1 + sdig/10. if keyword_set(fluxout) then begin if signoise gt detsigma then magstr = ' = '+epr(objapflux,2) + ' +/- ' + epr(objapfluxunc,2) $ else magstr = ' > '+epr(limflux,2)+ ' ['+clip(limsigma)+' sig] ' endif else begin ; 3-sigma image limit or, if object is >1 sigma, object+2sigma. if signoise gt detsigma then magstr = ' = '+fpr(objapmag,magform) + ' +/- ' + fpr(objapmagunc,errform) $ else magstr = ' > '+fpr(limmag,magform)+ ' ['+clip(limsigma)+' sig] ' endelse keyinfo = '' for k = 0, n_elements(printkey)-1 do begin fmt = '' if strpos(printkey[k],'JD') ge 0 then fmt='(F13.5)' keyinfo += ' ' + string(sxpar(h, printkey[k]),format=fmt) endfor sepstr = ':' if keyword_set(table) or nameobj[b] eq '' then sepstr = '' outstr = nameobj[b]+keyinfo+' '+sepstr+' '+clip(filter)+ magstr+ ' (+/- '+fpr(systematicunc,errform)+')' apstr = ' (r='+clip(fpr(radiusarcsec,2.2))+'"='+clip(fpr(radius,3.2))+'px, '+clip(a)+','+clip(d)+'='+clip(x+1)+','+clip(y+1)+')' if keyword_set(table) then $ outstr = repstr(repstr(repstr(repstr(repstr(outstr,'xxx',''),'+/-',''),'(',''),')',''),',','') + ' '+clip(fpr(radiusarcsec,2.2)) $ else $ outstr = outstr + apstr if keyword_set(silent) eq 0 then print, outstr outdata[nphot].objname = nameobj[b] outdata[nphot].filter = filter outdata[nphot].limit = signoise le detsigma outdata[nphot].mag = objapmag outdata[nphot].magunc = objapmagunc if n_elements(systematicunc) gt 0 then outdata[nphot].sysunc = systematicunc outdata[nphot].maglim = limmag outdata[nphot].flux = objapflux outdata[nphot].fluxunc = objapfluxunc outdata[nphot].fluxlim = limflux if ns gt 0 then begin outdata[nphot].ns_field = ns_field outdata[nphot].nvalid = nvalid outdata[nphot].nused = nused endif outdata[nphot].ra = a outdata[nphot].dec = d outdata[nphot].radiusarcsec = radiusarcsec outdata[nphot].x = x+1 outdata[nphot].y = y+1 outdata[nphot].radius = radius nphot += 1 endfor endif ; object photometry if nphot gt 0 then outdata = outdata[0:nphot-1] else outdata = -1 if nphot eq 0 and nobj gt 0 and keyword_set(silent) eq 0 then print, 'No objects in image.' !except = except clearerr = check_math() ; hide arithmetic warnings end