Get the latest version of ProFound and ProFit:
Set global evaluate (basically TRUE for GitHub version and FALSE for CRAN):
First load the libraries we need:
## Loading required package: FITSio## Loading required package: magicaxis## Loading required package: RcppOne 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=readFITS(system.file("extdata", 'VISTA_K.fits', package="magicaxis"))
VST_r=readFITS(system.file("extdata", 'VST_r.fits', package="magicaxis"))
GALEX_NUV=readFITS(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!ProFound has some useful class specific diagnostic plots:
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.
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.
VST_r_warpK=magwarp(VST_r, header_out=VISTA_K$hdr)
GALEX_NUV_warpK=magwarp(GALEX_NUV, header_out=VISTA_K$hdr)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.
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:
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):
And we can check they really do compute the same results:
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):
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')]
multi$cat_grp[1:5,c('groupID', 'mag_rg')]
merge_cat[1:5,c('segID', 'mag_rt', 'origin')]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:
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:
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!
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)
segim_NUV=profoundSegimWarp(segim_in=pro_K$segim, header_in=VISTA_K$hdr, header_out=GALEX_NUV$hdr)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)And now check the diagnostics:
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')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.