sky.estimate: Function to find local sky values

In AngusWright/LAMBDAR: Measure galaxy fluxes from an arbitrary FITS image using input aperture locations and dimensions

Function to find local sky values

Description

This function calculates local sky values using a combination source masking, sigma clipping of extreme pixels and iterating the median sky values of weighted annuli that are located around the ra and dec of the object of interest. The function is designed to cope with quite complex sky problems and is robust to outliers. It returns a single sky value and error for each requested ra and dec location.

Usage

sky.estimate(data.stamp,mask.stamp=NULL,rem.mask=FALSE,
data.stamp.lims=NULL, cutlo = 0, cuthi = 100,
radweight = 1, clipiters = 5, PSFFWHMinPIX = 0,
hardlo = 3, hardhi = 10, sigma.cut = 3, cat.x=NULL, cat.y=NULL,
mpi.opts = "")
plot.sky.estimate(data.stamp, mask.stamp = NULL, rem.mask = FALSE,
data.stamp.lims = NULL, cutlo = 0, cuthi = 100,
radweight = 1, clipiters = 5, PSFFWHMinPIX = 0,
hardlo = 3, hardhi = 10, sigma.cut = 3, plot.sci =
FALSE, contams = NULL, plot.all = FALSE, path = NULL, toFile = FALSE, 
cat.id = NULL, cat.x = NULL, cat.y = NULL, res = 220, bin.lwd = 1,
est.lwd = 1.5, all.est.lwd = 1.5)

Arguments

data.stamp	numeric (n,m) array or list [[i]] of numeric [j,k] arrays; The data image stamp (or stamps) to use in sky estimation.
mask.stamp	numeric (n,m) array or list [[i]] of numeric [j,k] arrays; The sourcemask image stamp (or stamps). This must be pixel matched to data.stamp, with pixel values of 1 where the pixels are likely to be sky, and 0 where they are likely to be sources. The SourceMask output of LAMBDAR can be used.
rem.mask	logical; should the maskim be used to remove known sources.
data.stamp.lims	numeric (i,4) array; The limits of the data stamps in the data image space.
cutlo	numeric (i) vector; Vector of low radius values (in pixels) where the code will start to calculate the sky annuli around the object. Should be large enough to avoid significant object flux, i.e. a few times the flux 90 radius.
cuthi	numeric (i) vector; Vector of low radius values (in pixels) where the code will stop calculating sky annuli around the object. Should be large enough to avoid significant object flux, i.e. a few times the flux 90 radius.
radweight	numeric; What radius power-law weighting should be used to bias the sky towards sky annuli nearer to the source. Larger values weight the sky value more towards central values (default 1, so linear weighting).
clipiters	numeric; How many 3-sigma clips of the sky will be made.
PSFFWHMinPIX	numeric; The FWHM of the PSF for origim in pixels.
hardlo	numeric; The lower limit for the inner sky annuli, where the limit used is hardlo * PSFFWHMinPIX.
hardhi	numeric; The lower limit for the outer sky annuli, where the limit used is hardhi * PSFFWHMinPIX.
sigma.cut	numeric; sigma-level used for iterative clipping of sky-pixels in estimation.
cat.x	numeric (i) vector; Vector of x-pixel values for the sky locations of interest (typically these will be the positions of objects).
cat.y	numeric (i) vector; Vector of y-pixel values for the sky locations of interest (typically these will be the positions of objects).
mpi.opts	list; options to be passed to the MPI foreach backend, if it is being used.
plot.sci	logical; should all science targets be plotted?
contams	logical (i) vector; are these sources contaminants (TRUE) or science targets (FALSE). Used by plot.sci.
plot.all	logical; should all objects be plotted?
path	character; path to use for output, if images are being written to file.
toFile	logical; should figures be sent to file?
cat.id	character (n) vector; id's of the sources.
res	numeric; resolution in pixels per inch of the output PNGs.
bin.lwd	numeric; Line width of bin-lines in output figures.
est.lwd	numeric; Line width of the 'estimate' lines, for bins used in final estimates, in output figures.
all.est.lwd	numeric; Line width of the 'estimate' lines for all bins in output figures.

Value

Returns a data frame with 5 columns containing the sky value (sky), sky error (skyerr) and the number of sky annuli that have error bars encompassing the final sky (Nnearsky), the RMS of the sky pixels (skyRMS), and the Pearson's Normality Test p-value for the sky pixels (skyRMSpval). Nnearsky can be at most 10, and will be 6-7 for a well behaved sky. skyRMSpval is a good diagnostic tool for determining when the sky background is highly non-gaussian, usually indicating the presence of an unrecognised contaminant (e.g. stellar diffraction spike), and can be used to screen out poorly derived/constrained fluxes.

Author(s)

Aaron Robotham ICRAR [e-mail]

Angus H Wright ICRAR [e-mail]

Examples

#Load LAMBDAR
library(LAMBDAR)

#Load sample image
data(SDSS.sample)
#Load sample catalogue
data(ApCat.sample)
#Save sample image to File
write.fits(file="SampleImage.fits",SDSS.sample)
#Read Astrometry
astr<-read.astrometry("SampleImage.fits")
#Determine the pixel-locations of the aperture catalogue objects
xy<-ad.to.xy(ra=ApCat.sample$RAdeg,dec=ApCat.sample$DECdeg,astr.struc=astr)
xy<-xy[which(xy[,1]>=1 & xy[,2]>=1 & xy[,1]<= astr$NAXIS[1] & xy[,2]<= astr$NAXIS[2]),]
#Get the Main Object's ID
main<-which.min(sqrt((xy[,1]-astr$NAXIS[1]/2)^2+(xy[,2]-astr$NAXIS[2]/2)^2))

#Perform a sky estimate for the Primary target
skyest<-sky.estimate(cat.x=xy[main,1],cat.y=xy[main,2],
data.stamp=SDSS.sample$dat[[1]],rem.mask=FALSE)
#Show Sky Estimate
print(skyest[,c('sky','skyerr','skyRMS','Nnearsky')])
#Note that the sky is likely non-gaussian:
print(skyest[,'skyRMSpval'])

#Perform a sky estimates for purely gaussian noise
gauss<-array(rnorm(length(SDSS.sample$dat[[1]])),dim=dim(SDSS.sample$dat[[1]]))
skyest.gauss<-sky.estimate(cat.x=xy[,1],cat.y=xy[,2], cutlo=rep(0,length(xy[,1])),
cuthi=rep(100,length(xy[,1])),data.stamp=gauss,rem.mask=FALSE)
#But the p-value is still systematically low:
hist(skyest.gauss[,'skyRMSpval'],breaks=20,xlab='Sky RMS Normality p-value',
col='red',density=20,angle=-45)
#This is becuase of the sigma-clipping:
skyest.gauss.noclip<-sky.estimate(cat.x=xy[,1],cat.y=xy[,2], cutlo=rep(0,length(xy[,1])),
cuthi=rep(100,length(xy[,1])),data.stamp=gauss,rem.mask=FALSE,clipiters=0)
hist(skyest.gauss.noclip[,'skyRMSpval'],breaks=20,xlab='Sky RMS Normality p-value',
add=TRUE,col='blue',density=10)
legend('topright',legend=c("With Clipping","Without Clipping"),col=c('red','blue'),
lty=1,pch=0)

#Plot the Sky Estimate
skyest<-plot.sky.estimate(cat.x=xy[main,1],cat.y=xy[main,2],
data.stamp=gauss,rem.mask=FALSE)

AngusWright/LAMBDAR documentation built on May 12, 2022, 1:49 a.m.