Ctrl+K

PSF photometry of NGC 2210

PSF photometry of NGC 2210#


“Astronomers never seem to want to do anything easy.”

- Peter B. Stetson (1987)

In this tutorial we will conduct some basic PSF photometry on a DESI Legacy Survey image of the globular cluster NGC 2210. We will use the photutils.psf module to perform PSF photometry on the image. Below image shows each frame of \(grz\) band and their false-color image of the NGC 2210. The aim for this tutorial is to get \(grz\) magnitudes and colors of stars inside this field, by PSF photometry, and to get their color-magnitude diagram for the investigation of the properties of NGC 2210.

Reading lists#

For the general understanding of the processes regarding to the astronomical photometric measurement, one can consult following references.

If you are not familiar with handling the astronomical data, the following tutorials may be helpful. Especially, the SNU_AOclass repository deveoloped by Yoonsoo P. Bach, which is the TA lecture note for astronomical observation class for SNU undergraduate students, will provide the most of the things you need, if you concern your lack of backgrounds.

PSF photometry - an overview#

PSF photometry is the process of measuring the brightness of the stars by numerically fitting their light profile across the CCD image frame. For the stars isolated in the frame, the aperture photometry suffices to measure their flux and uncertainty, without any heavy computations of optimization process, and without the possibility of failure in the model to represent the stellar profile inside the frame. As Stetson (1987) mentioned, the aperture photmetry is “the most easiest and potentially most accurate method for obtaining instrumental photometry for bright, uncrowded stars, because it involves simple counting of detected photons.” (Stetson, 1987) However, for stars inside the crowded field, such as open clusters, globular clusters, and Magellanic Clouds, it is impossible to measure the brightness of blended stars by just summing their photon counts inside some apertures. Thus we have to adopt numerical fitting technique to model the brightness profile of the blended stars simultaneously, under the assumption of the linearity of the detector, and extract the total flux of a star from the model.

PSF photometry is processed in the following steps:

  1. Subtract background light from the image

  2. Find stars inside the frame

  3. Define a model for stellar light profile (point spread function; PSF)

  4. Group the blended stars which have to be simultaneously fitted

  5. Fit the light profiles of the star groups

  6. Examine the residual (data - model) image to decide whether or not to …

    • change the parameters for finding stars

    • reconstruct the PSF model

    • resize fitting window

    • iterate additional extraction in the residual

    • etc.

This is the outline of this tutorial regarding to the above steps.

  1. We will subtract the background of each frame using sep package, which is developed to replicate Source Extractor (Bertin, 1996). Here we use sep.Background to subtract the background.

  2. We will first find stars with simple method find_peaks from photutils.detection. This simple selection of stars will provide the stars for modeling PSF. For better modeling of PSF, one should carefully pick the bright isolated stars, but here we just adopt some simple selection criteria to pick the stars, expecting that the adjacent stars will be eventually averaged out.

  3. With the stars we picked, we derive the effective PSF (EPSF), using EPSFBuilder from photutils.psf. We use photutils.psf.extract_stars to make EPSFStars object out of the table of stars we obtained by find_peaks, then feed this object to EPSFBuilder to get the EPSF model. For better model, we can subtract the adjacent stars in the window of EPSFStars by fit them using the EPSF model, and rederive the EPSF from the ‘cleaned’ EPSFStars.

  4. +5. BasicPSFPhotometry from photutils.psf will subtract background, find stars, group them, and fit the stellar groups to conduct PSF photometry. What we have to do is to select the background estimator (MMMBackground, which is the same with the one in DAOPHOT), star finder (DAOStarFinder; in here we input the FWHM measured from aperture photometry of stars found by find_peak, to better performance for finding stars), fitter (LevMarLSQFitter), etc.

  5. From BasicPSFPhotometry we can also obtain the residual image. In this image we can examine the goodness of fit, and identify some failures of the model (Stetson 1987):

    • images of stars which were in the frame but were not recognized by the star finder, perhaps due to having been blended with brighter companions

    • stars whose images contained defective pixels

    • images of galaxies which are almost but not quite, stellar

    • other similar departures from a perfect fit may be recognized and flagged

  6. We will iterate 1-6. over all bands (\(grz\)) to get PSF photometry in different bands.

  7. We match the stars with PSF photometry in each band. The match will be done by their positions. Here we just used the pixel coordinates, since the image we are dealing with is the cutout image from Legacy Survey so the images are aligned already. However, in more general case this is often not the case, so the match by celestial coordinates using WCS information should be better, otherwise you need to align the images first.

  8. We will obtain colors of the matched stars and draw color-magnitude diagram of them, which are in the field of the NGC 2210.

  9. In order to compare our PSF photometry with other results, we will use the photometry from the Legacy Survey catalog and the homoegeous photometric catalog provided by Stetson et al. (2019). To compare with the result of the latter, we will convert the DECam (the telescope used for DESI Legacy Survey) \(gr\) filter system to Johnson-Cousins \(BV\) system.

Caveat#

The DESI Legacy Survey image we adopt here is already standard-calibrated with each pixel in the unit of nanomaggy. In reality, you need to convert the instrumental flux into standard flux, with the instrumental flux of standard stars.

This tutorial is constructed with the limit of time, so this may not be complete and self-consistent, so if you have a question for this notebook ask me in person or send me an email to bahkhyeonguk@gmail.com.

%config InlineBackend.figure_format = 'retina'

No.

Software

Version

0

Python

3.12.2 64bit [Clang 14.0.6 ]

1

IPython

8.20.0

2

OS

macOS 13.4.1 arm64 arm 64bit

3

numpy

1.26.4

4

matplotlib

3.8.0

5

astropy

5.3.4

6

photutils

1.11.0

 
import copy
from pathlib import Path
from matplotlib import patches
from matplotlib import pyplot as plt
import numpy as np
import sys
from time import time, ctime
import warnings
from astropy.io import fits
from astropy.table import Table, hstack, vstack
from astropy.modeling.fitting import LevMarLSQFitter
from astropy.stats import sigma_clipped_stats, gaussian_fwhm_to_sigma, SigmaClip
from astropy.table import Table
from astropy.nddata import NDData, CCDData
from astropy.utils.exceptions import AstropyWarning
from astropy.visualization import simple_norm, ZScaleInterval
from photutils.aperture import CircularAperture, ApertureStats
from photutils.detection import DAOStarFinder, find_peaks
from photutils.psf import extract_stars, EPSFBuilder
from photutils.psf import DAOPhotPSFPhotometry, BasicPSFPhotometry
from photutils.psf import DAOGroup
from photutils.psf.epsf_stars import EPSFStars
from photutils.psf.utils import _extract_psf_fitting_names
from photutils.psf.utils import get_grouped_psf_model
from photutils.detection import IRAFStarFinder
from photutils.background import MMMBackground, MADStdBackgroundRMS
# from photutils.segmentation import make_source_mask
from photutils.segmentation import SegmentationImage
from photutils.segmentation import detect_threshold, detect_sources
from photutils.utils import circular_footprint
from mpl_toolkits.axes_grid1 import make_axes_locatable
import sep

sys.setrecursionlimit(10**7)
# plt.rcParams["font.family"] = 'Times New Roman'
plt.rcParams["font.size"] = 15
plt.rcParams["text.usetex"] = False
# plt.rcParams["mathtext.fontset"] = 'cm'

warnings.simplefilter('ignore', category=AstropyWarning)

WD = Path('.')
DATADIR = WD/'data'/'proj2'
OUTDIR = DATADIR/'out'

if not OUTDIR.exists():
    OUTDIR.mkdir()

1. Data Reading and Background Subtraction#

def zscale_imshow(ax, img, vmin=None, vmax=None, **kwargs):
    if vmin==None or vmax==None:
        interval = ZScaleInterval()
        _vmin, _vmax = interval.get_limits(img)
        if vmin==None:
            vmin = _vmin
        if vmax==None:
            vmax = _vmax
    im = ax.imshow(img, vmin=vmin, vmax=vmax, origin='lower', **kwargs)
    
    return im

def read_sci_data(fpath, show=True, newbyteorder=False, bkgsubtract=False):
    ccd_raw = CCDData.read(fpath, unit='nmgy') # be careful about the unit!
    ccd = CCDData.read(fpath, unit='nmgy')
    
    data = ccd.data
    
    data_mean, data_med, data_std = sigma_clipped_stats(data, sigma=2.)
    ccd.mask = data < data_mean - 3*data_std 
    ccd.data[ccd.mask] = data_med
    ccd.data -= data_med
    ccd.mask = (ccd.mask) & (data > 10.)
    if newbyteorder:
        ccd.data = ccd.data.byteswap().newbyteorder()
        
    n = len(ccd.data)
    if bkgsubtract:
        bkgsubs = []
        for i in range(n):
            bkgsub = background_substraction(ccd.data[i])
            bkgsubs.append(bkgsub)
        ccd.data = np.array(bkgsubs)
    
    if show:
        fig, axes = plt.subplots(1, n, figsize=(5*n, 5))
        for i in range(n):
            im = zscale_imshow(axes[i], ccd.data[i], cmap='gray')
            divider = make_axes_locatable(axes[i])
            cax = divider.append_axes('right', size='5%', pad=0.05)
            fig.colorbar(im, cax=cax, orientation='vertical')
            
            band = ccd.header[f'BAND{i}']
            
            axes[i].text(0.98, 0.02, f'${band}$-band',
                         transform=axes[i].transAxes,
                         va='bottom', ha='right', c='k',
                         bbox=dict(boxstyle='round', facecolor='w', alpha=0.9))
            
        plt.tight_layout()
        
    return ccd, ccd_raw

def background_substraction(img, box=100, filt=3, show=True):
    sigma_clip = SigmaClip(sigma=3.0, maxiters=10)
    threshold = detect_threshold(img, nsigma=3.0, sigma_clip=sigma_clip)
    segment_img = detect_sources(img, threshold, npixels=10)
    footprint = circular_footprint(radius=1)
    mask = segment_img.make_source_mask(footprint=footprint)
    # mask = SegmentationImage.make_source_mask(img, nsigma=2, npixels=1, dilate_size=3)
    bkg_sep = sep.Background(img.byteswap().newbyteorder(), 
                             mask=mask, bw=box, bh=box, fw=filt, fh=filt)
    bkgsub_sep = img - bkg_sep.back()
    
    if show:
        fig, axs = plt.subplots(2, 3, figsize=(10, 8))
        axs[1, 1].axis("off")
        
        data2plot = [
            dict(ax=axs[0, 0], arr=img,              title="Original data"),
            dict(ax=axs[0, 1], arr=bkg_sep.back(),   title=f"bkg (filt={filt:d}, box={box:d})"),
            dict(ax=axs[0, 2], arr=bkgsub_sep,       title="bkg subtracted"),
            dict(ax=axs[1, 0], arr=mask,             title="Mask"),
            # dict(ax=axs[1, 1], arr=bkg_sep.background_mesh,   title="bkg mesh"),
            dict(ax=axs[1, 2], arr=10*bkg_sep.rms(), title="10 * bkg RMS")
        ]
        
        for dd in data2plot:
            im = zscale_imshow(dd['ax'], dd['arr'])
            cb = fig.colorbar(im, ax=dd['ax'], orientation='horizontal')
            # cb.ax.set_xticklabels([f'{x:.3f}' for x in cb.get_ticks()])#, rotation=45)
            dd['ax'].set_title(dd['title'])
        
        plt.tight_layout()
    
    return bkgsub_sep

ccd, ccd_raw = read_sci_data(DATADIR/'NGC2210.grz.fits', bkgsubtract=True)
../_images/bb13bac535fe445e3abffb6fecd3685a3c4cdb45a84e2b749584b0a0e9e23046.png ../_images/2205cbf0939069b78e89edb76c32305e47fe74d3355f0e4fc709eb3af8aa8e28.png ../_images/399d97b2438993c88159a60acbe770e5c1f0c6dab50e1fe45d71db57f88692bb.png ../_images/f3f3e2c15b95c4b172150ce7280e413eab132554af53cfa5b985804a64eecd50.png

2. Selection of Stars for the PSF Estimation#

bands = ['g', 'r', 'z']
i = 0
band = bands[i]
data = ccd.data[i]
mask = ccd.mask[i]

bkgrms = MADStdBackgroundRMS()
std = bkgrms(data)
thres = 3*std

find_mask = np.zeros_like(data, dtype=bool)
# find_mask[250:650, 250:650] = True
# ll, hh = 0, 200
ll, hh = 250, 650
find_mask[ll:hh, ll:hh] = True
peaks_tbl = find_peaks(data, threshold=thres, mask=find_mask)
peaks_tbl['peak_value'].info.format = '%.8g'

print('number of peaks:', len(peaks_tbl))
peaks_tbl[:5].show_in_notebook()

# select stars within the cutout of specified size, which must be odd
size = 21 # typically ~4*FWHM. user should consider the PSF size for selecting this value
hsize = (size - 1) / 2
x = peaks_tbl['x_peak']  
y = peaks_tbl['y_peak']  
bound = ((x > hsize) & (x < (data.shape[1] -1 - hsize)) &
         (y > hsize) & (y < (data.shape[0] -1 - hsize)))
bright = (peaks_tbl['peak_value'] > 0.5) & (peaks_tbl['peak_value'] < 10.)

bbmask = (bound & bright)
bx, by = x[bbmask], y[bbmask]

isolated = [False if np.count_nonzero(np.sqrt((bx-xi)**2+(by-yi)**2)<size) > 1
            else True for xi, yi in zip(bx, by)]

mask_stars = isolated
# mask_stars = np.ones_like(bx, dtype=bool)

stars_tbl = Table()
stars_tbl['x'] = bx[mask_stars]  
stars_tbl['y'] = by[mask_stars]

print('number of stars picked:', len(stars_tbl))
stars_tbl[:5].show_in_notebook()

nddata = NDData(data=data, mask=ccd.mask[i])

# plot the locations of our selected stars
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
im = zscale_imshow(ax, data, cmap='gray')
ax.plot(stars_tbl['x'], stars_tbl['y'], '.r')
rect = patches.Rectangle((ll,ll), hh-ll, hh-ll, ec='tab:red', ls='--',
                         hatch='/', fc='none', lw=1)
ax.add_patch(rect)

number of peaks: 3781
number of stars picked: 107

<matplotlib.patches.Rectangle at 0x176215400>
../_images/0e05203c8ac74f7a562f03171418dafa9bfc8a1ac63aa213ace96c82a811b23f.png 
# extract cutouts of our selected stars
stars = extract_stars(nddata, stars_tbl, size=size)

# plot cutouts
nrows = 5
ncols = 5
fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(10, 10),
                       squeeze=True)
ax = ax.ravel()
for i in range(nrows * ncols):
    norm = simple_norm(stars[i], 'log', percent=99.0)
    ax[i].imshow(stars[i], norm=norm, origin='lower', cmap='viridis')

# derive fwhm from aperture statistics
cap = CircularAperture(stars.center_flat, size/5.)
apstat = ApertureStats(data, cap)

fwhm_mean, fwhm_med, fwhm_std = sigma_clipped_stats(apstat.fwhm.value,
                                                    sigma=3.0)
../_images/08f88b2dd1f33332beda2d5ef3ad78b25a4455b2653975d3bc60d62dc0924594.png

3. Derive Effective PSF#

def get_epsf(stars, band,
             oversampling=4, maxiters=10, smoothing_kernel='quadratic',
             show=True, ncols=None, figsize=(10, 15)):
    
    if show:
        if ncols == None:
            n, m = figsize
            ncols = n * np.round(np.sqrt(len(stars)/n/m)).astype(int)
        nrows = int(np.ceil(len(stars)/ncols))
        fig, axs = plt.subplots(nrows=nrows, ncols=ncols, figsize=figsize,
                               squeeze=True)
        axs = axs.ravel()
        for ax in axs:
            ax.axis('off')
        for i in range(len(stars)):
            norm = simple_norm(stars[i], 'log', percent=99.)
            axs[i].imshow(stars[i], norm=norm, origin='lower', cmap='viridis')
        fig.text(0.99,0.01,f'{band}-band',
                  ha='right', va='bottom', fontsize=15)
        # fig.suptitle(f'Picked stars to model PSF ({band})')
        plt.tight_layout(pad=0.3)
        
    epsf_builder = EPSFBuilder(oversampling=oversampling, maxiters=maxiters,
                               progress_bar=False,
                               smoothing_kernel=smoothing_kernel)  
    epsf, fitted_stars = epsf_builder(stars)
    
    if show:
        plt.figure()
        norm = simple_norm(epsf.data, 'log', percent=99.)
        plt.imshow(epsf.data, norm=norm, origin='lower', cmap='viridis')
        plt.title(f'Effective PSF ({band})')
        plt.colorbar()
        
    return epsf

epsf0 = get_epsf(stars, band)
../_images/6fa393ec87393a6c6ec7518b700915d4ebd031325952e88842c808f03578f1eb.png ../_images/4f01c83e2d69cae74d394ecaaa16ce5bc17c8002962b71017bd34de4afa1ae85.png

SUBSTAR - subtract adjacent stars and reconstruct EPSF#

# substar procedure
# subtract adjacent star profiles using the epsf model we obtained.

def get_substars(stars, psf_model, thres, fwhm, band, peakmax=None,
            show=True, ncols=10, figsize=(10, 15)):
    substars_list = []
    for star in stars:
        finder = DAOStarFinder(threshold=thres, fwhm=fwhm,
                               roundhi=5.0, roundlo=-5.0,
                               sharplo=0.0, sharphi=2.0, peakmax=peakmax)
        found = finder(star.data)
        found.rename_columns(['xcentroid', 'ycentroid', 'flux'],
                             ['x_0', 'y_0', 'flux_0'])
        if len(found) > 1:
            fitter = LevMarLSQFitter()
            y, x = np.indices(star.shape)
            xname, yname, fluxname = _extract_psf_fitting_names(psf_model)
            pars_to_set = {'x_0': xname, 'y_0': yname, 'flux_0': fluxname}
            group_psf = get_grouped_psf_model(psf_model, found, pars_to_set)
            fit_model = fitter(group_psf, x, y, star.data)
            
            if hasattr(star, 'cutout_center'):
                xc, yc = star.cutout_center
            elif hasattr(star, 'cutout_center_flat'):
                xc, yc = star.cutout_center_flat[0]
            else:
                raise AttributeError
            idx = np.argmin((found['x_0']-xc)**2+(found['y_0']-yc)**2) # crude
            cent_model = fit_model[idx]
            
            substar_data = star.data - fit_model(x, y) + cent_model(x, y)
            mask = np.ones_like(substar_data, dtype=bool)
            hsize = substar_data.shape[0]//2
            psize = int(fwhm)*2
            mask[hsize-psize:hsize+psize, hsize-psize:hsize+psize] = False
            _, med, sig = sigma_clipped_stats(substar_data[mask], sigma=2.0)
            # substar_data[substar_data < med-2*sig] = med
            substar_data = substar_data - med
            
            substar = copy.deepcopy(star)
            substar._data = substar_data
        else:
            substar = star
            
        substars_list.append(substar)
    
    substars = EPSFStars(substars_list)
    return substars

substars = get_substars(stars, epsf0, 2*std, fwhm_mean, band)
epsf1 = get_epsf(substars, band)
../_images/5fa435b48ea781ec990ba636f42dae4d9fededcf8d855c9fbbbe271ac2e416dc.png ../_images/54efeb838cba789b9e14fe30065fbbbae6fe9cd2120a021ab9eab27119a37757.png 
substars2 = get_substars(substars, epsf0, 1*std, fwhm_mean, band)
epsf2 = get_epsf(substars2, band)
../_images/49042ef11c439b7e93a2d161be23fefd184bea495c2fa5d97b99f5e311fce3be.png ../_images/53c21405e9bd7e2994aa5716e8881fd87dbb287868904a2068db546e4b9558c7.png

4. PSF Photometry with the EPSF#

def get_psfphot(data, mask, psf_model, fwhm, psf_size, sigma_thres=3.,
                peakmax=None):
    # set up the photometry object
    bkgrms = MADStdBackgroundRMS() 
    std = bkgrms(data) # background rms
    
    threshold = sigma_thres * std # threshold for star finding
    crit_separation = 2.0 * fwhm # critical separation for star grouping
    
    daogroup = DAOGroup(crit_separation) # grouping algorithm
    mmm_bkg = MMMBackground() # bkg estimator
    fitter = LevMarLSQFitter() # psf model fitter
    
    fsize = int(np.ceil(psf_size))
    fitshape = fsize if fsize%2 == 1 else fsize +1 # size of the fitting region
    
    # find stars for the initial guess
    daofind = DAOStarFinder(threshold=threshold,
                            fwhm=fwhm, #sigma_psf * gaussian_sigma_to_fwhm,
                            roundhi=5.0, roundlo=-5.0,
                            sharplo=0.0, sharphi=2.0, peakmax=peakmax)
    stars = daofind(data, mask=mask)
    
    # define the photometry object
    photometry = BasicPSFPhotometry(finder=daofind,
                                    group_maker=daogroup,
                                    bkg_estimator=mmm_bkg,
                                    psf_model=psf_model,
                                    fitter=fitter,
                                    aperture_radius=1.5*fwhm,
                                    # niters=1,
                                    fitshape=fitshape)

    # rename columns to match the names in the table
    stars.rename_columns(['xcentroid','ycentroid','flux'],
                         ['x_0', 'y_0', 'flux_0'])

    start = time() # to record the time
    print('start time -', ctime(start))
    
    # run the photometry
    result_tab = photometry.do_photometry(image=data, mask=mask,
                                          init_guesses=stars,
                                          progress_bar=True)
    residual_image = photometry.get_residual_image()
    
    finish = time()
    print('finish time -', ctime(finish))
    print(f'elapsed time - {(finish-start)/60:.2f} min')

    # discard stars outside the image
    shape = data.shape
    result_tab = discard_stars_outside(shape, result_tab)

    return result_tab, photometry


def discard_stars_outside(shape, result_tab):
    ny, nx = shape
    xin = np.logical_and(result_tab['x_fit'] > 0 , result_tab['x_fit'] < nx-1)
    yin = np.logical_and(result_tab['y_fit'] > 0 , result_tab['y_fit'] < ny-1)
    isin = np.logical_and(xin, yin)
    return result_tab[isin]

# set fwhm, psf_model, psf_size, sigma_thres, peakmax
fwhm = fwhm_med # median fwhm of selected stars
sigma_psf = fwhm * gaussian_fwhm_to_sigma # sigma of psf
psf_model = epsf0 # epsf model - just using the epsf without iteration
psf_size = fwhm * 4 # size of psf model to fit
sigma_thres = 5.0 # threshold for star detection
peakmax = 10.0 # maximum peak value for the detection (prevent saturated stars)

# get photometry by BasicPSFPhotometry
res = get_psfphot(data, mask, psf_model, fwhm, psf_size,
                  sigma_thres=sigma_thres, peakmax=peakmax)

result_tab, photometry = res # extract result table and photometry object
residual_image = photometry.get_residual_image() # residual_image = data-model

# save results
hdu = fits.PrimaryHDU(residual_image)
hdu.writeto(OUTDIR/f'res_{band}_{sigma_thres:.1f}.fits', overwrite=True)
result_tab.write(OUTDIR/f'result_tab_{band}_{sigma_thres:.1f}.csv',
                 format='csv', overwrite=True)

start time - Tue Apr  2 21:20:56 2024

finish time - Tue Apr  2 21:21:32 2024
elapsed time - 0.60 min

def show_psfphot_result(data, band, result_tab, residual_image):
    fig, axs = plt.subplots(1, 3, figsize=(10, 4))
    interval = ZScaleInterval()
    vmin, vmax = interval.get_limits(data)
    
    zscale_imshow(axs[0], data, vmin, vmax)
    zscale_imshow(axs[1], data - residual_image, vmin, vmax)
    zscale_imshow(axs[2], residual_image, vmin, vmax)
    axs[0].plot(result_tab['x_fit'], result_tab['y_fit'], marker='+', c='r',
                ms=1, ls='None')
    axs[2].text(0.95, 0.05, f'{band}-band', ha='right', va='bottom',
                   c='w', fontweight='bold', transform=axs[2].transAxes)
    axs[0].set_title('Original')
    axs[1].set_title('Model')
    axs[2].set_title('Residual')
    axs = axs.ravel()
    for ax in axs:
        ax.axis('off')
    
    plt.tight_layout()
    
    return fig

# plot results - data, model, residual
fig = show_psfphot_result(data, band, result_tab, residual_image)
../_images/49a795695b8dd1e36bdf58ae56cbb8d8142c7d620d78947e9ef39a4cb47de151.png

5. Iterations Over All Bands#

# Let's do the same for all bands
bands = ['g', 'r', 'z']
for i in range(1, len(bands)):
    band = bands[i]
    data = ccd.data[i]
    mask = ccd.mask[i]

    print(f'### Band: {band}')
    bkgrms = MADStdBackgroundRMS()
    std = bkgrms(data)
    thres = 3*std

    find_mask = np.zeros_like(data, dtype=bool)
    ll, hh = 250, 650
    find_mask[ll:hh, ll:hh] = True
    peaks_tbl = find_peaks(data, threshold=thres, mask=find_mask)
    peaks_tbl['peak_value'].info.format = '%.8g'

    # select stars within the cutout of specified size
    size = 21
    hsize = (size - 1) / 2
    x = peaks_tbl['x_peak']  
    y = peaks_tbl['y_peak']  
    bound = ((x > hsize) & (x < (data.shape[1] -1 - hsize)) &
            (y > hsize) & (y < (data.shape[0] -1 - hsize)))
    bright = (peaks_tbl['peak_value'] > 0.5) & (peaks_tbl['peak_value'] < 10.)

    bbmask = (bound & bright)
    bx, by = x[bbmask], y[bbmask]

    isolated = [False if np.count_nonzero(np.sqrt((bx-xi)**2+(by-yi)**2)<size) > 1
                else True for xi, yi in zip(bx, by)]

    mask_stars = isolated
    # mask_stars = np.ones_like(bx, dtype=bool)

    stars_tbl = Table()
    stars_tbl['x'] = bx[mask_stars]  
    stars_tbl['y'] = by[mask_stars]

    # print(stars_tbl)

    nddata = NDData(data=data, mask=ccd.mask[i])

    # plot the locations of our selected stars
    fig, ax = plt.subplots(1, 1, figsize=(10, 10))
    im = zscale_imshow(ax, data, cmap='gray')
    ax.plot(stars_tbl['x'], stars_tbl['y'], '.r')
    rect = patches.Rectangle((ll,ll), hh-ll, hh-ll, ec='tab:red', ls='--',
                            hatch='/', fc='none', lw=1)
    ax.add_patch(rect)

    # extract cutouts of our selected stars
    stars = extract_stars(nddata, stars_tbl, size=size)

    # derive fwhm from aperture statistics
    cap = CircularAperture(stars.center_flat, size/5.)
    apstat = ApertureStats(data, cap)

    fwhm_mean, fwhm_med, fwhm_std = sigma_clipped_stats(apstat.fwhm.value,
                                                        sigma=3.0)

    epsf0 = get_epsf(stars, band)
#     substars = get_substars(stars, epsf0, 2*std, fwhm_mean, band)
#     epsf1 = get_epsf(substars, band)
#     substars2 = get_substars(substars, epsf0, 1*std, fwhm_mean, band)
#     epsf2 = get_epsf(substars2, band)

    fwhm = fwhm_med
    sigma_psf = fwhm * gaussian_fwhm_to_sigma
    psf_model = epsf0
    psf_size = fwhm * 4
    sigma_thres = 5.
    peakmax=10.0

    res = get_psfphot(data, mask, psf_model, fwhm, psf_size,
                    sigma_thres=sigma_thres, peakmax=peakmax)

    result_tab, photometry = res
    residual_image = photometry.get_residual_image()

    hdu = fits.PrimaryHDU(residual_image)
    hdu.writeto(OUTDIR/f'res_{band}_{sigma_thres:.1f}.fits', overwrite=True)
    result_tab.write(OUTDIR/f'result_tab_{band}_{sigma_thres:.1f}.csv', format='csv', overwrite=True)

    show_psfphot_result(data, band, result_tab, residual_image)

### Band: r
start time - Tue Apr  2 21:21:39 2024

finish time - Tue Apr  2 21:22:18 2024
elapsed time - 0.67 min
### Band: z
start time - Tue Apr  2 21:22:24 2024

finish time - Tue Apr  2 21:22:43 2024
elapsed time - 0.32 min
../_images/92ee2be1216d5b756f06a9aceea5880750387c32e2d634e46a9e2ce98e56d93f.png ../_images/78904649e4b37ceb371af7e1f49a7c698bb955d6a63be7cba4cff6ac6c2bc875.png ../_images/65006ab37709c8143fd5396fb250cadffb2ce5851703320cf9c020a3f7825756.png ../_images/3e78d29cab339c0c17e1761899d2a417108120040c88a37ff497eb5d0253648d.png ../_images/2b924241a26db031180525f7a29202a0319e1b559980916bdf8eaa24aadf2379.png ../_images/16425af74a972ce695fc55076275cf61b7806fc7a330e902668134924d82a5f8.png ../_images/6a98f451b0d1f2df44f124b8dfe5e1154740a5b214865b31d12037263e487e60.png ../_images/2a6281deb6da89e9107bbda1d8ef80f7329b65640df70377186d0d5715499706.png 
# now let's bring the results together to match the stars from each band

sigma_thres = 5.0 # threshold for detection
tbls = []
for band in bands:
    tbls.append(Table.read(OUTDIR/f'result_tab_{band}_{sigma_thres:.1f}.csv'))

gtbl, rtbl, ztbl = tbls # extract the tables for each band

gtbl[:5].show_in_notebook()
# gtbl.show_in_notebook(display_length=10)
Table length=5
idxidx_0y_0sharpnessroundness1roundness2npixskypeakflux_0maggroup_idx_fity_fitflux_fitflux_uncx_0_uncy_0_unc
05851.45783192174191.24319228578171840.3622949160971678-0.43358245491981506-0.6376099061779756250.00.059649612754583362.6161646842956543-1.04416269701830135851.54952530214151.08945252819079391.57604421781860760.0206640364499113550.037530656401620320.04052453853709505
16867.37020331452471.72594938168105960.37489165251501855-0.463954895734787-0.4848240826943701250.00.13188531994819644.749582290649414-1.69163854151326826867.03462007063291.7659351982005464.0966441611949530.065579410683678020.044352856787225310.05080049441169921
27345.71174798059632.3251912338019090.4317357941222555-0.069928802549839020.07785421898175726250.00.15546333789825444.834078311920166-1.7107842042147177345.709658596607262.3191079800336213.905638608228380.016763198306766790.0123747570018171890.01151253322746038
38207.113030205389153.0587893116540890.4055527432918673-0.11554194986820221-0.04160403307501191250.00.706046342849731422.39531707763672-3.37539303914345328207.164547742643633.05723980025385317.7818739340906580.084215979729182380.0140420463583294640.013455827169577786
49256.4446457831383.2654800362278650.40959385970212470.07777648419141770.062090565073566986250.00.0446284115314483641.279773235321045-0.26783255815603899256.74697147189613.4773002757703751.3929294438785880.059015710677229110.121798575666866440.10961134586058725

6. Match the Results of Each Band#

# crude and slow match of the stars from each band
# for more sophisticated matching, see the exgalcosutils.catalog.match_catalogs function
# (https://github.com/hbahk/exgalcosutils/blob/d8dc13f05237bc6c21e187d3465f4b7a12cb3469/exgalcosutils/catalog.py#L16)

def dist(x1,y1,x2,y2):
    """
    Calculate the distance between two points.
    """
    return np.sqrt((x1-x2)**2+(y1-y2)**2)

def match_stars(reftab, objtab, table_names, verbose=False,
                objcard=['x_fit', 'y_fit', 'flux_fit', 'group_id'],
                refcard=['x', 'y'], rlim=1):
    """
    Return a table matched by positions of stars in each band.
    """
    tab_list = []
    for star in reftab: # for loop is slow, but it's ok for now
        xdata, ydata = objtab[objcard[0]], objtab[objcard[1]]
        match = dist(xdata, ydata, star[refcard[0]], star[refcard[1]]) < rlim
        tab = objtab[match]
        if len(tab) == 0:
            if verbose:
                print(star['ID'], ': nothing matched with the given standard star')
        else:
            if len(tab) > 1:
                if verbose:
                    print(star['ID'], ': duplication occured in the pixel match process')
                    print(tab[objcard[2], objcard[3]])
                tab.sort(objcard[2], reverse=True)
                tab = tab[0]
            tab_list.append(hstack([star, tab], table_names=table_names))
    
    matched_tab = vstack(tab_list)
    return matched_tab

print(f'g-band: {len(gtbl)} stars, r-band: {len(rtbl)} stars, z-band: {len(ztbl)} stars')

g-band: 2727 stars, r-band: 3309 stars, z-band: 1645 stars

# let's set the reference band to r-band, since it has the most stars.
# we will match the stars in the other bands to the r-band stars.
rlim = 3.0 # matching radius in pixel
grtbl = match_stars(rtbl, gtbl, ['r', 'g'], refcard=['x_fit', 'y_fit'], rlim=rlim)

grtbl[:5].show_in_notebook()
# grtbl.show_in_notebook(display_length=10)
Table length=5
idxid_rx_0_ry_0_rsharpness_rroundness1_rroundness2_rnpix_rsky_rpeak_rflux_0_rmag_rgroup_id_rx_fit_ry_fit_rflux_fit_rflux_unc_rx_0_unc_ry_0_unc_rid_gx_0_gy_0_gsharpness_groundness1_groundness2_gnpix_gsky_gpeak_gflux_0_gmag_ggroup_id_gx_fit_gy_fit_gflux_fit_gflux_unc_gx_0_unc_gy_0_unc_g
04216.25052360824571.29468353942724160.31965956104172266-0.34563782811164856-1.3672357013964573250.00.052184060215950011.208444356918335-0.205566645051835064207.394444733230132.936881552048430626.6021331731796380.40456092367224230.0456921616345036740.028435517423113728207.113030205389153.0587893116540890.4055527432918673-0.11554194986820221-0.04160403307501191250.00.706046342849731422.39531707763672-3.37539303914345328207.164547742643633.05723980025385317.7818739340906580.084215979729182380.0140420463583294640.013455827169577786
16851.46035553550271.2107791833192410.38046402004755103-0.45898571610450745-0.4806057426938176250.00.083797596395015722.6913681030273438-1.07493275256670346851.44574147679041.02328703232961041.9106907674503760.0246535882868330540.033063537852142630.040465068211044875851.45783192174191.24319228578171840.3622949160971678-0.43358245491981506-0.6376099061779756250.00.059649612754583362.6161646842956543-1.04416269701830135851.54952530214151.08945252819079391.57604421781860760.0206640364499113550.037530656401620320.04052453853709505
27867.30552122195731.693754032974010.4262485769040658-0.42902904748916626-0.4098428273941409250.00.175122618675231934.529127597808838-1.64003639027201997866.94015746106341.71323854308845674.64389049691466750.090062663345332720.051443768342771710.051049957808310126867.37020331452471.72594938168105960.37489165251501855-0.463954895734787-0.4848240826943701250.00.13188531994819644.749582290649414-1.69163854151326826867.03462007063291.7659351982005464.0966441611949530.065579410683678020.044352856787225310.05080049441169921
38345.61064400320752.30983777647632940.4469803029838169-0.14323532581329346-0.01388777237022494250.00.239335477352142335.749619483947754-1.89909775890231248345.58493318078532.29755320900552775.2299773916055470.0242958290545595180.0112125175766376580.0125343048589684347345.71174798059632.3251912338019090.4317357941222555-0.069928802549839020.07785421898175726250.00.15546333789825444.834078311920166-1.7107842042147177345.709658596607262.3191079800336213.905638608228380.016763198306766790.0123747570018171890.01151253322746038
49206.979029205891353.15899369643324060.44485950183475437-0.19737784564495087-0.06359504616173453250.01.340695023536682128.900096893310547-3.6522482470435819205.173442588456064.0632311763814946.30973452191956950.199661253248714840.0940200699752180.083533566493505348207.113030205389153.0587893116540890.4055527432918673-0.11554194986820221-0.04160403307501191250.00.706046342849731422.39531707763672-3.37539303914345328207.164547742643633.05723980025385317.7818739340906580.084215979729182380.0140420463583294640.013455827169577786
# let's check how well the matching went by plotting the difference of positions of the stars
# if the rlim=3.0 is too large, then try smaller value and repeat the matching process
# rlim ~< FWHM is recommended
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(111)
ax.plot(grtbl['x_fit_r']-grtbl['x_fit_g'], grtbl['y_fit_r']-grtbl['y_fit_g'], '.k', alpha=0.5)
ax.set_xlabel('$\Delta x$ (pixel)')
ax.set_ylabel('$\Delta y$ (pixel)')

fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111)
rgbimg = fits.getdata(DATADIR/'NGC2210.grz.rgb.fits')
ax.imshow(rgbimg, origin='lower')
ax.plot(grtbl['x_fit_r'], grtbl['y_fit_r'], 'xr', alpha=0.5)
ax.plot(grtbl['x_fit_g'], grtbl['y_fit_g'], 'xb', alpha=0.5)

<>:7: SyntaxWarning: invalid escape sequence '\D'
<>:8: SyntaxWarning: invalid escape sequence '\D'
<>:7: SyntaxWarning: invalid escape sequence '\D'
<>:8: SyntaxWarning: invalid escape sequence '\D'
/var/folders/c3/0f663kw90xldstvklyr5nq9c0000gn/T/ipykernel_6485/2619397371.py:7: SyntaxWarning: invalid escape sequence '\D'
  ax.set_xlabel('$\Delta x$ (pixel)')
/var/folders/c3/0f663kw90xldstvklyr5nq9c0000gn/T/ipykernel_6485/2619397371.py:8: SyntaxWarning: invalid escape sequence '\D'
  ax.set_ylabel('$\Delta y$ (pixel)')

[<matplotlib.lines.Line2D at 0x28abb49e0>]
../_images/be7085e37e045fc28e42d946e6df35ab58bf1bfc16da26d56432b1261ef59740.png ../_images/a74287151980108c06a3beec25c1fca35f7a6cd5f47c3ce8b203bbf30a4a0b9a.png

7. Color-Magnitude Diagram#

# then let's draw the color-magnitude diagram

def flux2mag(flux, zp=22.5): # the default zeropoint is set to 22.5 for the nanomaggy unit
    return -2.5*np.log10(flux) + zp

gmag = flux2mag(grtbl['flux_fit_g'])
rmag = flux2mag(grtbl['flux_fit_r'])

fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(121)
ax.scatter(gmag-rmag, rmag, marker='.', s=1.5, c='k', alpha=0.5)
ax.set_xlabel('$g-r$')
ax.set_ylabel('$r$')
ax.invert_yaxis()

# let's plot the color-magnitude diagram focusing on the non-outliers
ax2 = fig.add_subplot(122)
ax2.scatter(gmag-rmag, rmag, marker='.', s=5, c='k', alpha=0.5)
ax2.set_xlabel('$g-r$')
ax2.set_ylabel('$r$')
ax2.invert_yaxis()
ax2.set_xlim(-1, 1.5)
ax2.set_ylim(23, 15)

/var/folders/c3/0f663kw90xldstvklyr5nq9c0000gn/T/ipykernel_6485/3953004836.py:4: RuntimeWarning: invalid value encountered in log10
  return -2.5*np.log10(flux) + zp

(23.0, 15.0)
../_images/36ab45e9e90cdc7f609ba2f9cdd2ce77aeb2170a0d239eb6e7603cf6daa11066.png

8. Comparison with the References#

# let's compare the color-magnitude diagram with the legacy survey catalog
vertices = ccd.wcs.celestial.calc_footprint()
ralo, rahi = vertices[:,0].min(), vertices[:,0].max()
declo, dechi = vertices[:,1].min(), vertices[:,1].max()

lgs_url = 'https://www.legacysurvey.org/viewer/ls-dr10/cat.fits?'
cat_bbox = f'ralo={ralo:.4f}&rahi={rahi:.4f}&declo={declo:.4f}&dechi={dechi:.4f}'

query_url = lgs_url + cat_bbox
cathdu = fits.open(query_url)

cat = Table(cathdu[1].data)
cat[:5].show_in_notebook()
# cat.show_in_notebook(display_length=10)
Table length=5
idxreleasebrickidbricknameobjidbrick_primarymaskbitsfitbitstyperadecra_ivardec_ivarbxbydchisqebvmjd_minmjd_maxref_catref_idpmrapmdecparallaxpmra_ivarpmdec_ivarparallax_ivarref_epochgaia_phot_g_mean_maggaia_phot_g_mean_flux_over_errorgaia_phot_g_n_obsgaia_phot_bp_mean_maggaia_phot_bp_mean_flux_over_errorgaia_phot_bp_n_obsgaia_phot_rp_mean_maggaia_phot_rp_mean_flux_over_errorgaia_phot_rp_n_obsgaia_phot_variable_flaggaia_astrometric_excess_noisegaia_astrometric_excess_noise_siggaia_astrometric_n_obs_algaia_astrometric_n_good_obs_algaia_astrometric_weight_algaia_duplicated_sourcegaia_a_g_valgaia_e_bp_min_rp_valgaia_phot_bp_rp_excess_factorgaia_astrometric_sigma5d_maxgaia_astrometric_params_solvedflux_gflux_rflux_iflux_zflux_w1flux_w2flux_w3flux_w4flux_ivar_gflux_ivar_rflux_ivar_iflux_ivar_zflux_ivar_w1flux_ivar_w2flux_ivar_w3flux_ivar_w4fiberflux_gfiberflux_rfiberflux_ifiberflux_zfibertotflux_gfibertotflux_rfibertotflux_ifibertotflux_zapflux_gapflux_rapflux_iapflux_zapflux_resid_gapflux_resid_rapflux_resid_iapflux_resid_zapflux_blobresid_gapflux_blobresid_rapflux_blobresid_iapflux_blobresid_zapflux_ivar_gapflux_ivar_rapflux_ivar_iapflux_ivar_zapflux_masked_gapflux_masked_rapflux_masked_iapflux_masked_zapflux_w1apflux_w2apflux_w3apflux_w4apflux_resid_w1apflux_resid_w2apflux_resid_w3apflux_resid_w4apflux_ivar_w1apflux_ivar_w2apflux_ivar_w3apflux_ivar_w4mw_transmission_gmw_transmission_rmw_transmission_imw_transmission_zmw_transmission_w1mw_transmission_w2mw_transmission_w3mw_transmission_w4nobs_gnobs_rnobs_inobs_znobs_w1nobs_w2nobs_w3nobs_w4rchisq_grchisq_rrchisq_irchisq_zrchisq_w1rchisq_w2rchisq_w3rchisq_w4fracflux_gfracflux_rfracflux_ifracflux_zfracflux_w1fracflux_w2fracflux_w3fracflux_w4fracmasked_gfracmasked_rfracmasked_ifracmasked_zfracin_gfracin_rfracin_ifracin_zngood_gngood_rngood_ingood_zanymask_ganymask_ranymask_ianymask_zallmask_gallmask_rallmask_iallmask_zwisemask_w1wisemask_w2psfsize_gpsfsize_rpsfsize_ipsfsize_zpsfdepth_gpsfdepth_rpsfdepth_ipsfdepth_zgaldepth_ggaldepth_rgaldepth_igaldepth_znea_gnea_rnea_inea_zblob_nea_gblob_nea_rblob_nea_iblob_nea_zpsfdepth_w1psfdepth_w2psfdepth_w3psfdepth_w4wise_coadd_idwise_xwise_ylc_flux_w1lc_flux_w2lc_flux_ivar_w1lc_flux_ivar_w2lc_nobs_w1lc_nobs_w2lc_fracflux_w1lc_fracflux_w2lc_rchisq_w1lc_rchisq_w2lc_mjd_w1lc_mjd_w2lc_epoch_index_w1lc_epoch_index_w2sersicsersic_ivarshape_rshape_r_ivarshape_e1shape_e1_ivarshape_e2shape_e2_ivar
010000216060928m6923632True8192129PSF92.78757007578037-69.1295936521012662271860000000.099642075000000.01867.96483453.933143652.3 .. 0.00.07265162557283.3132659066759268.07856435075GE52797426863998853121.31090091.4812818-0.0707614352.38216115.7709455.5941812016.020.15112192.302127020.43511818.27049019.8973912.5708550False1.12946750.91282725000.0False0.00.01.09004870.91672444318.2985118.9451598.4879568.4328670.00.00.00.0877.3331718.1057268.1623590.842420.00.00.00.06.46511846.9689026.61270956.5697916.46464736.96839436.6122286.5693122.498601 .. 18.7982183.2165787 .. 27.843212.270678 .. 25.204992.5010936 .. 27.251497-0.094051026 .. 6.24606-0.1272347 .. 11.613896-0.034183294 .. 8.227281-0.013881785 .. 9.485866-0.094051026 .. 6.24606-0.1272347 .. 11.613896-0.034183294 .. 8.227281-0.013881785 .. 9.4858664302.66 .. 26.0971492988.5078 .. 31.5650921799.0396 .. 9.936828550.3695 .. 2.775580.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.47616625 .. 2.26276230.34674054 .. 4.482998-1.2738854 .. -13.678328-1.0870681 .. -34.9204330.46386555 .. 2.0955890.34032154 .. 4.3984013-1.2854753 .. -13.803032-1.0955067 .. -35.1557189.92067 .. 7.04280921.144802 .. 1.57573070.3320874 .. 0.0246302650.012552651 .. 0.00093064260.80649020.86513460.89894980.92216270.987763170.992467160.998388650.99939126434600001.94906972.69097571.07279011.03822850.00.00.00.00.00083586220.00049063960.00309495140.00229435720.00.00.00.00.00540391870.00.000288884560.172913270.99999991.00000011.01.0434500020480000001.4234461.2622231.5232911.46771552335.55131391.1311308.954497.814591479.8228810.94653198.5214759.838894.569283.54202565.3087114.8301644.569363.54208525.30879264.07281160.00.00.00.00.00.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00 .. 00 .. 00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00 .. 160 .. 160.00.00.00.00.00.00.00.0
110000216060928m6923638True8192129PSF92.78808919227178-69.1281897030704531441444000000.057867800000000.01865.4283473.224674006.31 .. 0.00.0725688957283.3132659066759268.07856435075GE52797427551302063362.6800222.46930810.69054791.45550073.08795863.27135852016.020.50525166.318427921.080729.514463020.03784611.6501010False0.79706230.24483252000.0False0.00.01.07425821.1786243954.2647637.32602078.6706229.57978250.00.00.00.0212.29167783.91205264.0349789.638540.00.00.00.03.31784085.6993956.7454497.4527453.3155165.69540176.7407237.4475231.3301907 .. 18.1853562.6378582 .. 26.6799472.2484665 .. 24.3603572.7981932 .. 25.01609-0.10030402 .. 3.5931304-0.11436806 .. 8.013494-0.10821519 .. 4.957432-0.06148577 .. 4.714269-0.10030402 .. 3.5931304-0.11436806 .. 8.013494-0.10821519 .. 4.957432-0.06148577 .. 4.714269927.5121 .. 21.6090373344.3213 .. 31.582361766.0658 .. 10.62638545.5569 .. 2.88994650.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.8132416 .. 4.69107530.7068993 .. 6.130646-2.7616856 .. -20.86119-0.3840622 .. -17.7332970.7977932 .. 4.11314630.6985612 .. 5.8887625-2.7767415 .. -21.121397-0.40319905 .. -18.03723787.94504 .. 6.90451720.665148 .. 1.56572160.3306132 .. 0.0246387680.0126051605 .. 0.00093475490.80668770.865277350.89905890.92224780.987777050.99247570.99839050.999392434600001.09341123.00500541.78607440.84815310.00.00.00.00.00282636210.00115974610.00348275340.00213750080.00.00.00.00.224189650.06.584817e-070.188449160.999181.01.01.0334520480020480000001.4234461.2622231.5232911.4677155273.656771391.1311308.954497.81459163.1253810.94653198.5214759.838894.55897333.54165985.3082784.8297160.534461143.54172065.30836064.07243440.00.00.00.00.00.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00 .. 00 .. 00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00 .. 160 .. 160.00.00.00.00.00.00.00.0
210000216060928m6923641True8192129PSF92.78849191182563-69.1379376102810747808930000000.073617224000000.01863.4283339.2837124841.4 .. 0.00.0731291757283.3132659066759268.07856435075GE52797426864162091520.25219831.4458102-0.0161461682.1834772.59723734.19035632016.020.36864175.9475927520.5099838.335049019.74629615.8874170False0.00.0000.0False0.00.01.3834041.0020562316.3670558.385989.3691338.9164710.00.00.00.01034.6625752.45465265.209584.309540.00.00.00.04.9505366.52029667.2847226.93276644.95030986.51999957.284396.93245031.7929345 .. 83.369242.8661578 .. 136.463152.3365793 .. 152.700652.6167917 .. 169.7976-0.1991595 .. 8.943009-0.26675004 .. 15.530283-0.17375015 .. 10.050449-0.16860354 .. 12.066279-0.1991595 .. 8.943009-0.26675004 .. 15.530283-0.17375015 .. 10.050449-0.16860354 .. 12.0662795254.3286 .. 51.785013218.978 .. 26.6044621767.5148 .. 9.691271486.22397 .. 2.59522680.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.07.5557785 .. 55.458234.264436 .. 33.605732-3.201238 .. -26.795685-11.259108 .. -30.2129617.5557785 .. 55.458234.264436 .. 33.605732-3.201238 .. -26.795685-11.259108 .. -30.21296140.17313 .. 4.20128715.625595 .. 1.32962570.33043286 .. 0.0245830360.012395947 .. 0.0009228850.80535090.86431120.898320560.92167160.987683240.99241780.99837810.99938726434600006.7988439.8141073.51140331.02281020.00.00.00.00.0146359660.0138848120.041220430.0386567820.00.00.00.00.00.00.000932372470.241418910.999999941.00.99999991.0434500020480000001.4234461.2622231.5232911.46771552335.55131391.1311308.954497.814591479.8228810.94653198.5214759.838894.57539133.56046685.29258444.82790144.575473.56052735.29266834.0709050.00.00.00.00.00.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00 .. 00 .. 00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00 .. 160 .. 160.00.00.00.00.00.00.00.0
310000216060928m6923644True8192129PSF92.78882908793277-69.1485166501719713858959000000.016943370000000.01861.74783193.922643837.08 .. 0.00.0737943957283.3132659066759268.07856435075GE52796957510081269761.30033350.5740968-0.745068130.71306241.36276950.96396982016.020.81724799.49190517221.0264439.564173020.6406827.57290840False0.00.0000.0False0.00.01.09338961.7364924954.30152134.62954854.45630844.39703560.00.00.00.0926.8981926.74927251.7507594.462150.00.00.00.03.3495443.60497523.47007513.42392023.34950693.60493523.47003653.4238821.3143821 .. 31.1908051.6728723 .. 54.665531.1416255 .. 58.5412251.3034669 .. 67.23993-0.032942876 .. 3.879783-0.046171136 .. 9.670687-0.032774374 .. 7.0183578-0.012071089 .. 11.460372-0.032942876 .. 3.879783-0.046171136 .. 9.670687-0.032774374 .. 7.0183578-0.012071089 .. 11.4603722076.3816 .. 6.5437463999.7795 .. 29.482911733.6534 .. 9.708909577.8104 .. 2.40309480.0 .. 0.096318650.0 .. 0.00.0 .. 0.00.0 .. 0.02.2776866 .. 42.421421.4591917 .. 27.440006-4.0022893 .. 4.625561-4.2769136 .. 41.475432.2776866 .. 42.421421.4591917 .. 27.440006-4.0022893 .. 4.625561-4.2769136 .. 41.4754364.83757 .. 4.724855418.852259 .. 1.37274840.32796517 .. 0.0244123380.012417219 .. 0.00092655170.803766550.86316550.89744480.92098810.98757190.992349150.99836330.99938166433600000.81985781.61311671.00322060.78441960.00.00.00.00.0071301650.0071517290.033385680.0232076740.00.00.00.00.203303890.000128919840.000949505830.190420520.999882940.970037760.80944580.99981916433500020480000001.4234461.2622231.55494651.46771552335.55131391.1311266.3957597.814591479.8228810.94653173.542659.838894.5759483.55020625.23828274.83654744.57609223.58011985.40936474.0793530.00.00.00.00.00.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00 .. 00 .. 00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00 .. 160 .. 160.00.00.00.00.00.00.00.0
410000216060928m6923647True8192129PSF92.78952269793871-69.138795899437814409476000000.028886644000000.01858.38173327.491559808.633 .. 0.00.0731751857283.3132659066759268.07856435075GE5279742686416267392-0.89152880.405525150.13063940.54060980.857314941.73573122016.020.747555114.86680616221.655887.766733020.3970835.16768840False0.01.6806109e-15000.0False0.00.00.89545151.9916974953.4685655.5902456.67087657.4199220.00.00.00.01355.5065872.0004274.7689573.059460.00.00.00.02.6932724.3407155.1798045.7614232.69371494.34142835.18065555.762371.0480822 .. 73.5212862.0099716 .. 127.483531.782606 .. 145.575912.078548 .. 163.38074-0.041219242 .. 9.795995-0.07765031 .. 16.660069-0.039027747 .. 11.232489-0.013793198 .. 13.646062-0.041219242 .. 9.795995-0.07765031 .. 16.660069-0.039027747 .. 11.232489-0.013793198 .. 13.6460627059.858 .. 53.4692733730.3472 .. 26.9583931835.2968 .. 10.0860815482.82944 .. 2.6782230.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.06.2135167 .. 59.957583.531377 .. 35.15216-0.92712647 .. -27.9808-7.43345 .. -41.604986.2135167 .. 59.957583.531377 .. 35.15216-0.92712647 .. -27.9808-7.43345 .. -41.6049840.100945 .. 4.04084515.733147 .. 1.31511030.33136323 .. 0.0245763680.01244733 .. 0.00092347660.80524120.86423190.898260.921624360.987675550.99241310.9983770.9993869434600001.52229481.7520351.12881280.65002070.00.00.00.00.0238207760.0195424580.063130070.0459956270.00.00.00.04.8863744e-090.06.8382506e-050.333077131.00.99999990.999999940.99999994434400020480000001.4234461.2622231.5232911.46771552335.55131391.1311308.954477.652081479.8228810.94653198.5214749.089194.57634833.55925875.29005864.82933764.57642753.559325.2901423.23273160.00.00.00.00.00.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00 .. 00 .. 00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00.0 .. 0.00 .. 160 .. 160.00.00.00.00.00.00.00.0
gcat = flux2mag(cat['flux_g']) # nanomaggy to magnitude
rcat = flux2mag(cat['flux_r'])

xcat, ycat = ccd.wcs.celestial.all_world2pix(cat['ra'], cat['dec'], 0)
cat['x'] = xcat
cat['y'] = ycat

fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(121)
ax.imshow(rgbimg, origin='lower')
ax.plot(xcat, ycat, 'xr')
ax.plot(grtbl['x_fit_g'], grtbl['y_fit_g'], '+b')
ax.set_xlabel('x (pixel)')
ax.set_ylabel('y (pixel)')
ax.set_xlim(0, len(rgbimg[0]))

ax2 = fig.add_subplot(122)
ax2.scatter(gmag-rmag, rmag, marker='.', s=5, c='k', alpha=0.5)
ax2.scatter(gcat-rcat, rcat, marker='+', s=5, c='r', alpha=0.5)
ax2.set_xlabel('$g-r$')
ax2.set_ylabel('$r$')
ax2.invert_yaxis()
ax2.set_xlim(-1, 1.5)
ax2.set_ylim(23, 15)

/var/folders/c3/0f663kw90xldstvklyr5nq9c0000gn/T/ipykernel_6485/3953004836.py:4: RuntimeWarning: invalid value encountered in log10
  return -2.5*np.log10(flux) + zp

(23.0, 15.0)
../_images/9b91286c64e523f05e38e7fc8f36b84c07b001b1b443bae0130131d83bc3bdef.png 
# let's match the stars in the catalog with the stars of our photometry
cmtbl_r = match_stars(rtbl, cat, ['our', 'cat'], refcard=['x_fit', 'y_fit'],
                    objcard=['x', 'y', 'flux_r'], rlim=rlim)

# and compare the magnitude
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(111)
crmag = flux2mag(cmtbl_r['flux_r']) # catalog r-band magnitude
ormag = flux2mag(cmtbl_r['flux_fit']) # our r-band magnitude
ax.plot(crmag, ormag, '.k', alpha=0.5)
ax.set_ylim([14.5, 21.5])
ax.set_xlim([14.5, 21.5])
ax.plot([14, 22], [14, 22], ls='--', c='tab:red', zorder=0)
ax.set_xlabel('catalog r-band magnitude')
ax.set_ylabel('our r-band magnitude')

diff = ormag - crmag
offset, scatter = np.nanmean(diff), np.nanstd(diff)

ax.text(0.95, 0.05, f'offset = {offset:.2f}\n scatter = {scatter:.2f}',
        transform=ax.transAxes, va='bottom', ha='right')

/var/folders/c3/0f663kw90xldstvklyr5nq9c0000gn/T/ipykernel_6485/3953004836.py:4: RuntimeWarning: invalid value encountered in log10
  return -2.5*np.log10(flux) + zp

Text(0.95, 0.05, 'offset = 0.16\n scatter = 0.56')
../_images/ecdfccb9d00c22f9c234a0f29efb619f30d37eb821e6936e6bd832f761db62fd.png 
# doing the same thing in gz-band
cmtbl_g = match_stars(gtbl, cat, ['our', 'cat'], refcard=['x_fit', 'y_fit'],
                      objcard=['x', 'y', 'flux_g'], rlim=rlim)
cmtbl_z = match_stars(ztbl, cat, ['our', 'cat'], refcard=['x_fit', 'y_fit'],
                      objcard=['x', 'y', 'flux_z'], rlim=rlim)

fig, axes = plt.subplots(1, 3, figsize=(15,5))
cgmag = flux2mag(cmtbl_g['flux_g']) # catalog g-band magnitude
ogmag = flux2mag(cmtbl_g['flux_fit']) # our g-band magnitude
czmag = flux2mag(cmtbl_z['flux_z']) # catalog z-band magnitude
ozmag = flux2mag(cmtbl_z['flux_fit']) # our z-band magnitude
for ax, cmag, omag, band in zip(axes, [cgmag, crmag, czmag],
                                [ogmag, ormag, ozmag], ['g', 'r', 'z']):
    ax.plot(cmag, omag, '.k', alpha=0.5)
    ax.set_ylim([14.5, 22])
    ax.set_xlim([14.5, 22])
    ax.plot([14, 22], [14, 22], ls='--', c='tab:red', zorder=0)
    ax.set_xlabel(f'catalog {band}-band magnitude')
    ax.set_ylabel(f'our {band}-band magnitude')
    
    diff = omag - cmag
    offset, scatter = np.nanmean(diff), np.nanstd(diff)

    ax.text(0.95, 0.05, f'offset = {offset:.2f}\n scatter = {scatter:.2f}',
            transform=ax.transAxes, va='bottom', ha='right')

/var/folders/c3/0f663kw90xldstvklyr5nq9c0000gn/T/ipykernel_6485/3953004836.py:4: RuntimeWarning: invalid value encountered in log10
  return -2.5*np.log10(flux) + zp
../_images/88a3295a8268b82450e47f27330604cde0103f40122302cebe3d50180c001dd4.png

There are systematic offsets between cataloged magnitude and our PSF magnitude, for each band. What made these offsets? Who derived more correct result? There are many possible reasons of the offsets one can think of:

  • Difference in aperture correction process between our method and Tractor, which is the tool for extracting PSF magnitude in DESI Legacy Survey catalog.

  • Difference in PSF model - only a percent level error in PSF model can 10~100% errors of the faint star near the bright stars. The systematic offset in PSF model can derive the systematic offsets in the magnitude.

  • Difference in the background subtraction

  • And many more -

One possible way to resolve the cause of the offset is compare the PSF photometry at the un-crowded field, to find out how the crowding affects the resulting magnitudes. Another way is simulate a star inside the image and repeat the photometry over, and compare the input with the magnitudes from the PSF photometry, which is called the artificial star test.

I will leave this task to the reader.

To compare with another reference with the \(BV\) magnitude, we need a transformation relation of our \(gr\) magnitude into \(BV\) magnitude.

Transformation from DECam \(griz\) to Johnson-Cousins (\(UBV R_C I_C\)) system

\[B = g_{\rm DES} + 0.371 (g-r)_{\rm DES} + 0.197 \ \ [-0.5 < (g-r)_{\rm DES} ≤ 0.2] \ \ ({\rm RMS}: 0.022 {\rm mag})\]
\[B = g_{\rm DES} + 0.542 (g-r)_{\rm DES} + 0.141 \ \ [ 0.2 < (g-r)_{\rm DES} ≤ 0.7] \ \ ({\rm RMS}: 0.017 {\rm mag})\]
\[B = g_{\rm DES} + 0.454 (g-r)_{\rm DES} + 0.200 \ \ [ 0.7 < (g-r)_{\rm DES} ≤ 1.8] \ \ ({\rm RMS}: 0.059 {\rm mag})\]
\[V = g_{\rm DES} - 0.465 (g-r)_{\rm DES} - 0.020 \ \ [-0.5 < (g-r)_{\rm DES} ≤ 0.2] \ \ ({\rm RMS}: 0.012 {\rm mag})\]
\[V = g_{\rm DES} - 0.496 (g-r)_{\rm DES} - 0.015 \ \ [ 0.2 < (g-r)_{\rm DES} ≤ 0.7] \ \ ({\rm RMS}: 0.011 {\rm mag})\]
\[V = g_{\rm DES} - 0.445 (g-r)_{\rm DES} - 0.062 \ \ [ 0.7 < (g-r)_{\rm DES} ≤ 1.8] \ \ ({\rm RMS}: 0.024 {\rm mag})\]
def Bmag(gmag, rmag):
    B = np.zeros_like(gmag)
    gmr = gmag - rmag
    blu = (gmr > -0.5) & (gmr <= 0.2)
    grn = (gmr > 0.2) & (gmr <= 0.7)
    red = (gmr > 0.7) & (gmr <= 1.8)
    B[blu] = gmag[blu] + 0.371*gmr[blu] + 0.197
    B[grn] = gmag[grn] + 0.542*gmr[grn] + 0.141
    B[red] = gmag[red] + 0.454*gmr[red] + 0.200
    return B

def Vmag(gmag, rmag):
    V = np.zeros_like(gmag)
    gmr = gmag - rmag
    blu = (gmr > -0.5) & (gmr <= 0.2)
    grn = (gmr > 0.2) & (gmr <= 0.7)
    red = (gmr > 0.7) & (gmr <= 1.8)
    V[blu] = gmag[blu] - 0.465*gmr[blu] - 0.020
    V[grn] = gmag[grn] - 0.496*gmr[grn] - 0.015
    V[red] = gmag[red] - 0.445*gmr[red] - 0.062
    return V

B, V = Bmag(gmag, rmag), Vmag(gmag, rmag)

# let's compare the color-magnitude diagram of ours to the reference
# reference: Peter Stetson's stellar standards and software - Homogeneous photometry
# one can download the data from https://www.cadc.hia.nrc.gc.ca/en/community/STETSON/homogeneous/

# read the reference data
reftbl = Table.read(DATADIR/'stetson'/'NGC2210.dat', format='ascii.fixed_width_two_line')

# plot the color-magnitude diagram
fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111)
ax.scatter(reftbl['B']-reftbl['V'], reftbl['V'], marker='.', s=1, c='k', alpha=0.3, label='Stetson')
ax.scatter(B-V, V, marker='.', s=5, c='r', alpha=0.5, label='Our result')
ax.set_ylim([14.5, 22])
ax.set_xlim([-1, 2])
ax.set_xlabel('$B-V$')
ax.set_ylabel('$V$')
ax.invert_yaxis()
ax.legend(markerscale=5)

<matplotlib.legend.Legend at 0x176307f20>
../_images/da743e0322fa1da5eda59d0cbf67e86ecd1232db59798f374a5b1ccf2c647f19.png

From the color-magnitude diagram we obtained, which properties of this system will you be able to find out?

  • Properties of stars inside this field

  • Properties of the globular cluster NGC 2210

  • Is our photometry reliable?

  • etc

Further Readings#

  • Stetson (1987PASP…99…191S), DAOPHOT: A Computer Program for Crowded-Field Stellar Photometry

    Seminal paper about how to find stars and do photometry from a CCD image over crowded field

  • Anderson & King (2000PASP…112.1360A), Toward High-Precision Astrometry with WFPC2. I. Deriving an Accurate Point-Spread Function

    This work describes how the effective point-spread function is estimated.

  • Stetson et al. (2019MNRAS.485.3042S), Homogeneous photometry - VII. Globular clusters in the Gaia era

    Modern homogeneous photometry given by Peter Stetson, especially for the globular clusters, by conducting the PSF photometry.

  • Baumgardt et al. (2023MNRAS.521.3991B), Evidence for a bottom-light initial mass function in massive star clusters

    This reference may gave you some idea about how member stars in a cluster can be identified.

Contents