ProFound: Colour Me Happy
Aaron Robotham
2023-03-15
Get the latest version of ProFound and
ProFit:
library(devtools)
install_github('asgr/ProFound')
install_github('ICRAR/ProFit')
Set global evaluate (basically TRUE for GitHub version and FALSE for
CRAN):
First load the libraries we need:
## Loading required package: Rfits
## Loading required package: magicaxis
## Loading required package: Rcpp
library(Rfits)
library(Rwcs)
It’s a Colourful Life
One of the big design considerations when developing
ProFound was to make colour photometry easy and
flexible. This vignette will discuss a couple of ways to tackle it.
First we need some images to do colour photometry on (these are taken
from the ProFound Robotham et al 2018 paper):
VISTA_K = Rfits_read_image(system.file("extdata", 'VISTA_K.fits', package="magicaxis"))
VST_r = Rfits_read_image(system.file("extdata", 'VST_r.fits', package="magicaxis"))
GALEX_NUV = Rfits_read_image(system.file("extdata", 'GALEX_NUV.fits', package="magicaxis"))
Let’s take a look at the images.
It is clear they are pretty well world coordinate system (WCS)
aligned, but have very different pixel scales. This creates challenges
for decent colour photometry, but there are a few ways to tackle it with
the tools in ProFound.
First let’s run ProFound in default blind mode on each image and just
match the nearby segments (this is the worst approach!)
pro_K=profoundProFound(VISTA_K, magzero=30)
pro_r=profoundProFound(VST_r, magzero=0)
pro_NUV=profoundProFound(GALEX_NUV, magzero=20.08) #Ugly zero point I know- see Driver et al 2016 Table 3!
## dim(image)[1]/box[1] must be >=3, box[1] modified to 28
## dim(image)[2]/box[2] must be >=3, box[2] modified to 28
## dim(image)[1]/grid[1] must be >=3, grid[1] modified to 28
## dim(image)[2]/grid[2] must be >=3, grid[2] modified to 28
ProFound has some useful class specific diagnostic
plots:
## [1] "Too many same valued pixels: turning off magmap scaling!"
Do not worry too much about the differnce in the source finding, the
point really is doing a catalogue match would return weird results for
the bright central galaxy because there are such different blind
photometry solutions. So we need to enforce some restrictions to get
better colour photometry.
Image Pixel Warping
One route is to warp the images onto a common WCS scheme and run
ProFound in either full segment or bright segment mode.
Here we will take the VIKING K-band as our target WCS image.
library(ProPane)
VST_r_warpK = propaneWarp(VST_r, keyvalues_out=VISTA_K$keyvalues)
GALEX_NUV_warpK = propaneWarp(GALEX_NUV, keyvalues_out=VISTA_K$keyvalues)
The images are now interpolated onto a common WCS, where their
surface brightness is properly maintained (so we do not gain or lose
flux). The small differences below are because the original images did
not precisely cover the K-band image WCS.
## [1] 3.321829e-06
## [1] 0
## [1] 13.16952
sum(GALEX_NUV_warpK$image)
## [1] 0
We can now easily run ProFound in matched segment
mode, turning the dilation iterations off:
pro_r_warpK=profoundProFound(VST_r_warpK, segim=pro_K$segim, magzero=0, iters=0)
pro_NUV_warpK=profoundProFound(GALEX_NUV_warpK, segim=pro_K$segim, magzero=20.08, iters=0)
And now check the diagnostics:
magplot(pro_NUV_warpK$segstats$mag-pro_K$segstats$mag, pro_r_warpK$segstats$mag-pro_K$segstats$mag, xlab='NUV-K', ylab='r-K')
We can achieve the above in a simple way using the
profoundMultiBand function, just setting iters_tot=0.
This automates many of the fiddly steps and allows for multi band
stacking (if desired):
multi=profoundMultiBand(inputlist=list(GALEX_NUV_warpK, VST_r_warpK, VISTA_K), magzero=c(20.08,0,30), iters_tot=0, detectbands='K', multibands=c('NUV','r','K'))
## *** Will use K for source detection ***
## *** Will use NUV r K for multi band photometry ***
## *** Magzero: NUV=20.08 r=0 K=30 ***
## *** Gain: not specified for any bands, so shot-noise will be ignored ***
## *** Will compute total multi band photometry ***
## *** Will compute isophotal colour multi band photometry ***
## *** Will compute grouped segment multi band photometry ***
## *** Currently processing single detection band K ***
## *** Currently processing multi band NUV ***
## * Processing total photometry *
## * Processing colour photometry *
## * Processing group photometry *
## *** Currently processing multi band r ***
## * Processing total photometry *
## * Processing colour photometry *
## * Processing group photometry *
## *** Currently processing multi band K ***
## * Processing total photometry *
## * Processing colour photometry *
## * Processing group photometry *
And we can check they really do compute the same results:
magplot(pro_NUV_warpK$segstats$mag, multi$cat_tot$mag_NUVt, grid=TRUE)
Oher Colour Outputs
You might have noticed that profoundMultiBand
mentions measuring three different types of photometry: total
photometry, bright isophotal colour photometry and grouped segment
photometry. Total photometry uses the dilated total photometry
segmentation map, whereas the isophotal photometry only uses the
brighter segim_orig (pre-dilation) map. Grouped segment photometry is a
bit more complicated- it represents the flux for groups of touching
segments (either grouped by touching segim_orig or dilated segim,
depending on the ‘groupby’ setting in
profoundMultiBand). In some cases the groups might
represent better photometry, e.g. consider the dilated grouped segments
(this is like running with tolerance=Inf):
profoundSegimPlot(multi$pro_detect$image, multi$pro_detect$group$groupim)
The three catalogues are included in the list output of
profoundMultiBand, and are named ‘cat_tot’, ‘cat_col’
and ‘cat_grp’. Mostly the grouped photometry will be the same as the
total, since mostly things are not touching pre-dilation:
magplot(multi$cat_grp$mag_rg, multi$cat_tot[match(multi$cat_grp$groupID, multi$cat_tot$segID),"mag_rt"], grid=TRUE, log='', xlab='Group Mag', ylab='Total Mag')
abline(0,1,col='red')
There are two clear outliers, one of which is the central bright
source (segment/group ID 1), which in our grouped catalogue has become
merged with its neighbouring star (segment ID 3). If we did
want to merge this together there is a handy utility function to do so,
where ‘groupID_merge’ is a vector of all the groups we prefer over the
segmented versions of the photometry (in this case just groupID 1):
merge_cat=profoundCatMerge(segstats=multi$cat_tot, groupstats=multi$cat_grp, groupsegID=multi$pro_detect$group$groupsegID, groupID_merge=1)
multi$cat_tot[1:5,c('segID', 'mag_rt')]
## segID mag_rt
## 1 1 14.70499
## 2 2 16.06213
## 3 3 15.47166
## 4 4 16.95286
## 5 5 16.30285
multi$cat_grp[1:5,c('groupID', 'mag_rg')]
## groupID mag_rg
## 1 1 14.26944
## 2 2 16.06213
## 3 4 16.95286
## 4 5 16.30285
## 5 6 20.52013
merge_cat[1:5,c('segID', 'mag_rt', 'origin')]
## segID mag_rt origin
## 1 1 14.26944 group
## 2 2 16.06213 seg
## 4 4 16.95286 seg
## 5 5 16.30285 seg
## 6 6 20.52013 seg
We gain an origin flag in this new catalogue making it clear what the
origin of the photometry is (‘group’ for groupstats, or ‘seg’ for
segstats). Notice how in the merged catalogue segment ID 3 has been
removed, since it is merged in with segment ID 1 to form groupID=1. All
this grouping information is stored in the groupsegID object:
multi$pro_detect$group$groupsegID[1,"segID"][[1]]
## [1] 1 3
We might also only want to use segments from the colour catalogue
that exist in this new merger catalogue. This is easy to do with base
R:
merge_cat_col=multi$cat_col[multi$cat_col$segID %in% merge_cat$segID,]
These now both have the same number of rows by construction, and
sources will appear on the same row, so we can plot comparisons very
easily:
magplot(merge_cat$mag_rt, merge_cat_col$mag_rc, xlab='Total Mag', ylab='Colour Mag', grid=TRUE)
abline(0,1,col='red')
As a note, there is one object that is actually fainter in its
dilated total photometry- it is pretty faint, and probably not a good
object!
Segmentation Map Warping
The alternative approach is to leave the pixels be, but warp the
segmentation map itself to fit a target WCS. Handily there is a function
to do exactly this!
segim_r=profoundSegimWarp(segim_in=pro_K$segim, header_in=VISTA_K$hdr, header_out=VST_r$hdr)
## Warning in base::":"(from, to): numerical expression has 52 elements: only the
## first used
## Warning in base::":"(from, to): numerical expression has 52 elements: only the
## first used
segim_NUV=profoundSegimWarp(segim_in=pro_K$segim, header_in=VISTA_K$hdr, header_out=GALEX_NUV$hdr)
## Warning in base::":"(from, to): numerical expression has 52 elements: only the
## first used
## Warning in base::":"(from, to): numerical expression has 52 elements: only the
## first used
We can now run ProFound with a warped segmentation
map:
pro_r_warpK2=profoundProFound(VST_r, segim=segim_r, magzero=0, iters=0)
pro_NUV_warpK2=profoundProFound(GALEX_NUV, segim=segim_NUV, magzero=20.08, iters=0)
## dim(image)[1]/box[1] must be >=3, box[1] modified to 28
## dim(image)[2]/box[2] must be >=3, box[2] modified to 28
## dim(image)[1]/grid[1] must be >=3, grid[1] modified to 28
## dim(image)[2]/grid[2] must be >=3, grid[2] modified to 28
And now check the diagnostics:
## [1] "Too many same valued pixels: turning off magmap scaling!"
Note we cannot now guarantee that we have exactly the same number of
segments since some small ones might not even cover a single pixel. This
means we need to match back by segID.
magplot(pro_r_warpK$segstats$mag[match(pro_r_warpK2$segstats$segID,pro_r_warpK$segstats$segID)], pro_r_warpK2$segstats$mag, grid=TRUE, xlab='r Image Warp / mag', ylab='r Segim Warp / mag')
magplot(pro_NUV_warpK$segstats$mag[match(pro_NUV_warpK2$segstats$segID,pro_NUV_warpK$segstats$segID)], pro_NUV_warpK2$segstats$mag, grid=TRUE, xlab='NUV Image Warp / mag', ylab='NUV Segim Warp / mag')
Better colours
You might want to allow some degree of segmentation map growth (set
iters>0) for you target photometry. This is particularly true when
the target band PSF is broader than the detection band. Also, you might
do better using the brighter part of the segmentation map returned by
ProFound, i.e. pass segim_orig rather than segim. In
this case you will need to adjust your magnitudes for the detection band
by the segstat$origfrac value, which gives the flux fraction in the
original segment compared to the final one returned. I.e. something like
profound$segstat$mag - 2.5*log10(profound$segstat$origfrac).
As you can see, there is a lot of flexibility to how colours can be
computed- either in a very static forced mode (as above), or more
dynamically to better adapt to the different characteristics of the
target band data. This approached has been used to good effect on
UV-radio data, so success should be possible with a bit of care and
thought.