ProFound: Stacking Images

Stacked Photometry

If you have many images of the same area (even in different bands) a photometrically driven user can create a S/N stacked image weighted using the inverse of \({sky_{RMS}}^2\). The reason we weight it this way is you actually want to weight by exposure time, which correlates with inverse variance (traditional Source Extractor weight maps are inverse variance for this reason), which is \(1/{sky_{RMS}}^2\).

To convince yourself consider you have 8 exposures with sky RMS equal to 16 originally: 4 you combine, so the sky RMS becomes \(\frac{1}{\sqrt{\frac{1}{16^2}+\frac{1}{16^2}+\frac{1}{16^2}+\frac{1}{16^2}}}=\frac{1}{\sqrt{\frac{4}{16^2}}}=8\). Now later on you want to combine the 5 images you have. Clearly the optimal weighting you can possibly get will be \(\frac{1}{\sqrt{\frac{8}{16^2}}}=\frac{8}{\sqrt{2}}=5.65\) (the direct stack of the original 8 images). With our already stacked image added to our 4 others we can achieve the same S/N by weighting our stacks by the inverse variance: \(\frac{1}{\sqrt{\frac{1}{16^2}+\frac{1}{16^2}+\frac{1}{16^2}+\frac{1}{16^2}+\frac{1}{8^2}}}=\frac{1}{\sqrt{\frac{4}{16^2}+\frac{1}{8^2}}}=\frac{1}{\sqrt{\frac{8}{16^2}}}=\frac{8}{\sqrt{2}}=5.65\).

ProFound comes with a handy function that does all the weight-watching for us called profoundMakeStack. In this vignette we will make a simulated image (with a trivially predictable improvement in S/N) and check to see that our stacking behaves as expected.

First we generate a random image with 200 stars and 200 extended sources. The value used roughly correspoond to the source densities and magnitude distributions you might expect to find in a Z-band VIKING frame (this was used to derive the image statistics).

set.seed(666)

ExamplePSF=profitMakeGaussianPSF(fwhm=5)
ExamplePSF=ExamplePSF/sum(ExamplePSF)

Ngal=200
Nstar=200

model_test=list(
    sersic=list(
        xcen=runif(Ngal,0,1000),
        ycen=runif(Ngal,0,1000),
        mag=24-rpareto(Ngal,2),
        re=rpois(Ngal,5)+runif(Ngal),
        nser=runif(Ngal,1,4),
        ang=runif(Ngal,0,180),
        axrat=runif(Ngal,0.3,1),
        box=runif(Ngal,-0.3,0.3)
    ),
    pointsource=list(
        xcen=runif(Nstar,0,1000),
        ycen=runif(Nstar,0,1000),
        mag=24-rpareto(Nstar,1.5)
    )
)

model_test$sersic$mag[model_test$sersic$mag<15]=runif(length(which(model_test$sersic$mag<15)),15,22)
model_test$pointsource$mag[model_test$pointsource$mag<15]=runif(length(which(model_test$pointsource$mag<15)),15,22)

im_test<-profitMakeModel(modellist=model_test, psf=ExamplePSF, dim=c(1000,1000), magzero = 30)$z

Now we can add noise in a realistic manner. Image-1 is the nominal depth of a single Z-band VIKING image, whereas the others are each one quarter the depth (which gives the variance, hence a factor of two in standard-deviation) so the four of them should achieve the same depth when stacked.

im_test1=im_test+rnorm(1e6,sd=sqrt(im_test))
im_test1=im_test1+rnorm(1e6,sd=10)
im_test2a=im_test+rnorm(1e6,sd=2*sqrt(im_test))
im_test2a=im_test2a+rnorm(1e6,sd=2*10)
im_test2b=im_test+rnorm(1e6,sd=2*sqrt(im_test))
im_test2b=im_test2b+rnorm(1e6,sd=2*10)
im_test2c=im_test+rnorm(1e6,sd=2*sqrt(im_test))
im_test2c=im_test2c+rnorm(1e6,sd=2*10)
im_test2d=im_test+rnorm(1e6,sd=2*sqrt(im_test))
im_test2d=im_test2d+rnorm(1e6,sd=2*10)

We now need to run ProFound on the 5 images separately:

pro_test1=profoundProFound(im_test1, skycut=1)
pro_test2a=profoundProFound(im_test2a, skycut=1)
pro_test2b=profoundProFound(im_test2b, skycut=1)
pro_test2c=profoundProFound(im_test2c, skycut=1)
pro_test2d=profoundProFound(im_test2d, skycut=1)

We can now stack the images based on the sky and sky-RMS:

stack=profoundMakeStack(image_list = list(im_test1, im_test2a, im_test2b, im_test2c, im_test2d), sky_list = list(pro_test1$sky, pro_test2a$sky, pro_test2b$sky, pro_test2c$sky, pro_test2d$sky), skyRMS_list = list(pro_test1$skyRMS, pro_test2a$skyRMS, pro_test2b$skyRMS, pro_test2c$skyRMS, pro_test2d$skyRMS), magzero_in=c(30,30,30,30,30), magzero_out=30)

Finally we can re-run ProFound to check the depth achieved. If the input standard-deviation of image-1 is 10, and each of the other 4 images is 20, then the final achieved sky-RMS should be \(\frac{1}{\sqrt{\frac{1}{10^2}+\frac{1}{20^2}+\frac{1}{20^2}+\frac{1}{20^2}+\frac{1}{20^2}}}=\frac{1}{\sqrt{\frac{4}{20^2}+\frac{1}{10^2}}}=\frac{1}{\sqrt{\frac{8}{400}}} \sim 7.07\).

First run ProFound:

pro_stack=profoundProFound(stack$image, skycut=1)

Then can check the distribution of the sky pixels directly:

maghist(pro_stack$skyRMS, grid=TRUE)
abline(v=7.07, col='red')

We look to have achieved our expected depth of a sky RMS close to 7.07!