<jupyter_start><jupyter_text>Photometry with Photutils Photutils is an affiliated package of Astropy to provide tools for detecting and performing photometry of astronomical sources.It contains tools for:* **Background Estimation (photutils.background)*** **Source Detection (photutils.detection)*** Grouping Algorithms* **Aperture Photometry (photutils.aperture)*** PSF Photometry (photutils.psf)* PSF Matching (photutils.psf.matching)* Image Segmentation (photutils.segmentation)* Centroids (photutils.centroids)* Morphological Properties (photutils.morphology)* Geometry Functions (photutils.geometry)* Datasets (photutils.datasets)* Utility Functions (photutils.utils)The content of the current exercise has been adapted from this [Python course held at ICE](https://github.com/Python4AstronomersAndParticlePhysicists/PythonWorkshop-ICE) in 2017. As usual, the first step is to import the required packages.<jupyter_code>import warnings
import os.path
from glob import glob

from astropy import wcs
from astropy.io import fits

# FITS files are rather old and do not conform to FITS standard
warnings.filterwarnings("ignore", category=wcs.FITSFixedWarning)

%run plotfits.ipynb # Load plotting function<jupyter_output><empty_output><jupyter_text>Astrometry Images of M67 have no astrometric calibration in header and, therefore, pixel positions cannot be directly transformed to sky coordinates. The usual way to calibrate a science images is by using [Astrometry.net](http://astrometry.net/) [(Lang et al., 2010)](https://ui.adsabs.harvard.edu/abs/2010AJ....139.1782L/abstract). However, the installation procedure in Windows is quite complex and requires a quite large amount of time. Therefore, the instructions to astrometrically calibrate the M67 images are shown here, but only the final producs provided in resources will be used.  The commands to obtain the astrometric calibration of the M67 images are:```cd resources/red_datasolve-field bf_m67_b_1.fit --config ../resources/astrometry/astrometry.cfgsolve-field bf_m67_b_2.fit --config ../resources/astrometry/astrometry.cfgsolve-field bf_m67_b_3.fit --config ../resources/astrometry/astrometry.cfgsolve-field bf_m67_y_1.fit --config ../resources/astrometry/astrometry.cfgsolve-field bf_m67_y_2.fit --config ../resources/astrometry/astrometry.cfgsolve-field bf_m67_y_3.fit --config ../resources/astrometry/astrometry.cfgcp *.wcs ../resources/astrometry```The result of this execution already provides astrometrically calibrated FITS images (named `*.new`), but in order to illustrate the use of WCS functionality, the `*.wcs` files will be used. These files contain the tranformation from pixel positions to celestial coordinates in a format called [World Coordinate System](https://docs.astropy.org/en/stable/wcs/) (WCS), a standard for FITS files. The goal is to introduce this information into the clean M67 images.<jupyter_code>for wcs_calib in glob(os.path.join('..', 'resources', 'astrometry', '*.wcs')):
    basename = os.path.basename(wcs_calib).replace('wcs', 'fit')
    clean_name = os.path.join('..','resources', 'red_data', basename)
    calibrated_name = clean_name.replace('bf_', 'bfa_')

    wcs_header = fits.open(wcs_calib)[0].header
    clean_image = fits.open(clean_name)
    
    clean_image[0].header.extend(wcs_header, update=True)
    clean_image.writeto(calibrated_name, overwrite=True)

plot(calibrated_name)<jupyter_output><empty_output><jupyter_text>Once astrometric calibration is introduced in the images, some information can be retrieved:<jupyter_code>import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord

wcs_info = wcs.WCS(clean_image[0].header)
pixel_scaling = wcs.utils.proj_plane_pixel_scales(wcs_info) * u.degree
pixel_scale = np.mean(pixel_scaling).to(u.arcsec)
print('Pixel scale is:', pixel_scale)

corners = SkyCoord(wcs_info.calc_footprint(), unit=(u.degree, u.degree))
corners.to_string('hmsdms')<jupyter_output><empty_output><jupyter_text>Photometry Image statistics Before performing any photometry we need to have a first guess of the image background properties of our science images of M67. Let's get the basic statistics of it. For that we will need to remove the sources using a sigma clipping routine:<jupyter_code>from astropy.stats import sigma_clipped_stats, median_absolute_deviation
from astropy.io import fits
from astropy import table

rows = []
for filename in glob(os.path.join('..', 'resources', 'red_data', 'bfa_*.fit')):
    print('Reading', filename)
    extension = fits.open(filename)[0]
    values = list(sigma_clipped_stats(extension.data, sigma=3.0, maxiters=5))
    rows.append([filename] + values + [extension.data.T])

columns = ['Filename', 'Mean', 'Median', 'Deviation', 'Image']
statistics = table.Table(rows=rows, names=columns)
statistics.sort('Filename')
statistics.show_in_notebook()<jupyter_output><empty_output><jupyter_text>Object detectionTo detect the sources inside a astronomical image [Photutils](https://photutils.readthedocs.io/en/stable/) provides several implementations. We will use `DAOStarFinder`:<jupyter_code>from photutils import DAOStarFinder
sigma_detection = 5.0
fwhm = 1.5 * u.arcsec / pixel_scale

sources = []
for image in statistics:
    filename = image['Filename']
    std = image['Deviation']
    median = image['Median']
    daofind = DAOStarFinder(fwhm=fwhm.value, threshold=sigma_detection*std)
    photometry = daofind(image['Image'] - median)
    photometry['Filename'] = [filename] * len(photometry)
    sources.append(photometry)

catalogs = table.vstack(sources).group_by('Filename')
catalogs[:10].show_in_notebook()<jupyter_output><empty_output><jupyter_text>In order to validate the result, the data of the first image can be plotted:<jupyter_code>%matplotlib widget
import matplotlib.pyplot as plt

filename = catalogs.groups.keys[0][0]

plot(filename)

sources = catalogs.groups[0]
plt.scatter(sources['xcentroid'], sources['ycentroid'], alpha=0.5, color='limegreen')<jupyter_output><empty_output><jupyter_text>Some of the objects detected are not real stars, but CCD artifacts. The best way to identify them is to evaluate their `sharpness` and `roundness`. We are going to select only the stars with `sharpness` and `roundness` around the median values:<jupyter_code>valid = []
for filename, sources in zip(catalogs.groups.keys, catalogs.groups):
    sharpness = sources['sharpness']
    roundness = sources['roundness1']

    sharp_med = np.median(sharpness)
    sharp_mas = median_absolute_deviation(sharpness)
    round_med = np.median(roundness)
    round_mas = median_absolute_deviation(roundness)

    mask = (sharpness - sharp_med)**2 < 9 * sharp_mas**2 
    mask &= (roundness - round_med)**2 < 9 * round_mas**2
    mask &= sources['peak'] < 50000

    valid += list(mask)

    print(f'Found {len(sources[mask])} of {len(sources)} '
          f'valid stars in {filename[0]}')

catalogs.add_column(valid, index=1, name='valid')
    
fig = plt.figure()
plt.scatter(sharpness, roundness)
plt.scatter(sharpness[mask], roundness[mask], color='red')
plt.show(fig)<jupyter_output><empty_output><jupyter_text>To illustrate what has been done, the first image can be plotted again:<jupyter_code>plot(filename[0])

plt.scatter(sources['xcentroid'], sources['ycentroid'], alpha=0.5, color='limegreen')
plt.scatter(sources['xcentroid'][~mask], sources['ycentroid'][~mask], alpha=0.5, color='red')
catalogs[:10].show_in_notebook()<jupyter_output><empty_output><jupyter_text>As it can be seen, stars not fulfilling the criteria are considered invalid. An important aspect to be addressed is that the images looks properly calibrated, but the photometry table does not contain the sky positions of the stars. To solve that `wcs` can be used:<jupyter_code>coords = None
for filename, sources in zip(catalogs.groups.keys, catalogs.groups):
    wcs_transform = wcs.WCS(fits.open(filename[0])[0].header)
    new_pos = SkyCoord.from_pixel(sources['ycentroid'], sources['xcentroid'],
                                  wcs_transform)

    if coords:
        coords = coords.insert(len(coords), new_pos)
    else:
        coords = new_pos

if 'ra' in catalogs.colnames:
    del catalogs['ra']
if 'dec' in catalogs.colnames:
    del catalogs['dec']

catalogs.add_column(coords.ra, index=2, name='ra')
catalogs.add_column(coords.dec, index=3, name='dec')

catalog_name = os.path.join('..', 'resources', 'red_data', 
                            'instrumental_photometry.txt')
catalogs.write(catalog_name, format='ascii.fixed_width', overwrite=True)

catalogs[:10].show_in_notebook()<jupyter_output><empty_output><jupyter_text>PSF Modelling (Not given in the course, to be completed)We assumed that the image had a typical value of 1.5" per pixel. But we can make a more accurate estimation by fitting the pixels to a moffat profile around a detected star. First, let's select the source we want to use for PSF modelling. For that, we are going to select the brightest **valid** star:<jupyter_code>psf_sources = catalogs[:0]
for filename, sources in zip(catalogs.groups.keys, catalogs.groups):
    sources.sort('mag')
    selected = sources[sources['valid']][0]
    psf_sources.add_row(selected)
    
psf_sources.show_in_notebook()<jupyter_output><empty_output><jupyter_text>To simplify the problem, we turn the 2D profile into a 1D distance array from the center of each pixel to the centroid of the source estimated by DAOPhot. As we intend to find the profile of the source, we need to remove the possible sky background that lies behind:<jupyter_code>flux_counts = []
pixel_distance = []

# Median bkg subtracted image
bkg_subtracted_image = image - median

x_cen = int(isource['xcentroid'])
y_cen = int(isource['ycentroid'])

# Pixels around detection loop
analysis_radius = 10
for x in range(x_cen - analysis_radius, x_cen + analysis_radius):
    for y in range(y_cen - analysis_radius, y_cen + analysis_radius):
        flux_counts.append(((bkg_subtracted_image[y][x]) / isource['peak']))
        pixel_distance.append(np.linalg.norm((isource['xcentroid'] - x, isource['ycentroid'] - y)))<jupyter_output><empty_output><jupyter_text>Here we present two possible models that can fit the PSF distribution. A Gaussian and a Moffat profile:<jupyter_code>from astropy.modeling import models, fitting

model = 'moffat'

if model == 'gaussian':
    # Fit the data using a Gaussian
    fwhm_best_guess = 1.5
    model_init = models.Gaussian1D(amplitude=1.0, mean=0, stddev=fwhm_best_guess)
elif model == 'moffat':
    # Fit the data using a Moffat
    model_init = models.Moffat1D(amplitude=1.0, x_0=0, gamma=2., alpha=3.5)
else:
    raise Exception("Unknown model type: %s. Must be gaussian or moffat." % model)

fitter = fitting.SimplexLSQFitter()
fitted_model = fitter(model_init, pixel_distance, flux_counts)

print ("Fit value:",  fitter.fit_info['final_func_val'])
print ("SN:", isource['flux'] * sigma_detection)<jupyter_output><empty_output><jupyter_text>Once fitted the models, we need to convert from the parameters to the actual FWHM estimate.<jupyter_code># FHWM conversion
if model == 'gaussian':
    iFWHM = 2.355 * fitted_model.stddev * sdss_pixelscale.value
elif model == 'moffat':
    iFWHM = fitted_model.gamma * 2 * np.sqrt(2 ** (1. / fitted_model.alpha) - 1) * sdss_pixelscale.value
else:
    raise Exception("Unknown model type: %s. Must be gaussian or moffat." % model)

print ("FWHM estimated: ", iFWHM)<jupyter_output><empty_output><jupyter_text>We can finally plot and see how the model traces the pixel values traces our fitted model.<jupyter_code># Check fitting
if fitter.fit_info['final_func_val'] < 5.0:
    color = 'green'
else:
    color = 'red'
    
# Plot the data with the best-fit model
plt.figure()
plt.plot(pixel_distance, flux_counts, 'ko')
rx = np.linspace(0, int(max(pixel_distance)), int(max(pixel_distance)) * 10)
plt.plot(rx,
         fitted_model(rx),
         color=color,
         lw=3.0,
         label='SN: %.2f, Fit: %.2f, FWHM: %.2f"' % (isource['flux'] * sigma_detection,
                                                       fitter.fit_info['final_func_val'],
                                                       iFWHM))
plt.xlabel('Distance (pixels)')
plt.ylabel('Normalized flux (ADU)')
plt.title('%s profile fit' % model)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Background modelling As we have seen in the case with non-uniform background, the constant median can be insuficient. Here we produce a 2D model of the background that can be subtracted from the original image to improve the modelling of the stars close to a very large extended source (or when the backrgound is not flat for any other reason)<jupyter_code>from photutils import Background2D, SigmaClip, MedianBackground
sigma_clip = SigmaClip(sigma=3., iters=10)
bkg_estimator = MedianBackground()
bkg = Background2D(data=sdss_g_hdu_list[0].data, 
                   box_size=(50, 50), 
                   filter_size=(3, 3),
                   sigma_clip=sigma_clip, 
                   bkg_estimator=bkg_estimator)
my_python_ds9(bkg.background)<jupyter_output><empty_output><jupyter_text>Now let's go back to where the background was subtracted to verify the difference.  Aperture photometry We will work with the simple circular apertures. First we need to set the size we want and create the aperture objects:<jupyter_code>from photutils import CircularAperture

aperture_radius = 5.0

positions = (sources['xcentroid'], sources['ycentroid'])
apertures = CircularAperture(positions, r=aperture_radius)

my_python_ds9(sdss_g_hdu_list[0].data[0:400, 0:400])
apertures.plot(color='limegreen', lw=2)<jupyter_output><empty_output><jupyter_text>Global Background subtraction As we only want the flux from the sources we are interested in, we need to remove the contribution from the background. The simplest way is to remove the median value constant of the sigma clipped image that we calculated before to the entire array:<jupyter_code>from photutils import aperture_photometry

bkg_subtracted_image = sdss_g_hdu_list[0].data - median
phot_table_global_bkg = aperture_photometry(bkg_subtracted_image, apertures)
print (phot_table_global_bkg)<jupyter_output><empty_output><jupyter_text>Local Background subtraction We could also remove the 2D background model that we calculated before, but it is usually more precise to create an annulus around the source we are interested in, and estimate the background level there.<jupyter_code>from photutils import CircularAnnulus

r_in = 10
r_out = 15

annulus_apertures = CircularAnnulus(positions, r_in=r_in, r_out=r_out)

my_python_ds9(sdss_g_hdu_list[0].data[0:400, 0:400])
apertures.plot(color='limegreen', lw=2)
annulus_apertures.plot(color='purple', lw=2, alpha=1)<jupyter_output><empty_output><jupyter_text>The circular apertures and the annulus can be joined for a common photometry processing.<jupyter_code>apers = [apertures, annulus_apertures]
phot_table_local_bkg = aperture_photometry(sdss_g_hdu_list[0].data, apers)
print(phot_table_local_bkg)<jupyter_output><empty_output><jupyter_text>Now we use the aperture_sum_1 to estimate the level of background around the source. We need to know the area of the annulus for this estimation:<jupyter_code>bkg_mean = phot_table_local_bkg['aperture_sum_1'] / annulus_apertures.area()<jupyter_output><empty_output><jupyter_text>And finally we can remove the background estimation to all pixels in the aperture:<jupyter_code>bkg_sum = bkg_mean * apertures.area()
final_sum = phot_table_local_bkg['aperture_sum_0'] - bkg_sum
phot_table_local_bkg['residual_aperture_sum'] = final_sum
print(phot_table_local_bkg['residual_aperture_sum'])<jupyter_output><empty_output><jupyter_text>In this comparison we see that many sources have the same flux with both methods but a significant amount of sources (the ones in the galaxy halo) have significantly more flux in the global subtraction method, as the flux from M102 is added to the one of the stars.<jupyter_code>plt.scatter(phot_table_local_bkg['residual_aperture_sum'], phot_table_global_bkg['aperture_sum'])
plt.xlim(0,100)
plt.ylim(0,100)
plt.xlabel('Local background')
plt.ylabel('Global background')<jupyter_output><empty_output><jupyter_text>Let's make a plot to verify this! Astropy qtables work similarly to Pandas, so we can make iloc searches.  Exercise PhotUtils: Identify problematic sourcesLocate which sources differ their flux in more than 50% between the measurement with local and global background estimation. Plot the results on the image.TIP: Using the "pandas-like" tools can facilitate the selection<jupyter_code># %load -r 43-52 solutions/10_Astronomy.py
#phot_table_local_bkg.add_index('id')
very_different = phot_table_local_bkg[phot_table_local_bkg['residual_aperture_sum'] * 1.5
                                           < phot_table_global_bkg['aperture_sum']]

quite_similar = phot_table_local_bkg[phot_table_local_bkg['residual_aperture_sum'] * 1.5
                                           > phot_table_global_bkg['aperture_sum']]

my_python_ds9(sdss_g_hdu_list[0].data)
plt.scatter(very_different['xcenter'], very_different['ycenter'], c='red', alpha=0.5)
plt.scatter(quite_similar['xcenter'], quite_similar['ycenter'], c='limegreen', alpha=0.5)<jupyter_output><empty_output><jupyter_text>Accessing to annulus pixels By default when calling aperture_photometry, we only receive the sum of pixels and therefore we can only access to the mean of the pixel values inside an aperture. To properly measure the local background, we would need to get the median of the pixels in the annulus and preferably, perform a sigma clipping over the values.Recently photutils enabled masks so we can actually do that!<jupyter_code># Flat annulus - 11
# Source in annulus - 12

source_num = 12

x = phot_table_local_bkg[source_num]['xcenter'].value
y = phot_table_local_bkg[source_num]['ycenter'].value

stamp_radius = 20
my_python_ds9(sdss_g_hdu_list[0].data[int(y - stamp_radius):
                                      int(y + stamp_radius), 
                                      int(x - stamp_radius): 
                                      int(x + stamp_radius)])
plt.grid('off')<jupyter_output><empty_output><jupyter_text>We can access to the annulus mask:<jupyter_code># Create annulus mask
annulus_apertures = CircularAnnulus((x, y), r_in=r_in, r_out=r_out)
masks = annulus_apertures.to_mask(method='center')
m0 = masks[0]

plt.imshow(m0)
plt.grid('off')<jupyter_output><empty_output><jupyter_text>And apply it to the data:<jupyter_code># Cut to data
cutout_data = m0.apply(sdss_g_hdu_list[0].data)

plt.imshow(cutout_data)
plt.grid('off')<jupyter_output><empty_output><jupyter_text>Now that we have access to the pixels, we can perform the median and compare to the mean<jupyter_code>annulus_array = cutout_data[cutout_data != 0]

# Apply statistics to masked data
mean = np.mean(annulus_array)
median = np.median(annulus_array)
std = np.std(annulus_array)
print (mean, median, std)<jupyter_output><empty_output><jupyter_text>And even sigma-clip the sources that may appear on top of the background. This creates a numpy.masked array:<jupyter_code>from astropy.stats import sigma_clip

clip_annulus_array = sigma_clip(annulus_array, sigma=3, iters=2)

mean_clipped = np.ma.mean(clip_annulus_array)
median_clipped = np.ma.median(clip_annulus_array)
std_clipped = np.ma.std(clip_annulus_array)
print (mean_clipped, median_clipped, std_clipped)<jupyter_output><empty_output><jupyter_text>Exercise PhotUtils 2: Estimate aperture errorCompute what would be the measurement error for the last aperture (the one with sigma clipping)$$\sigma^2_{F} = \sum_i{\sigma^2_{i}+F/g}$$TIP: Gain is the exposure time when pixels are in e/s<jupyter_code># %load -r 56-61 solutions/10_Astronomy.py
g = float(sdss_g_hdu_list[0].header['EXPTIME'])
aperture_area = apertures.area()  # or also = math.pi * aperture_radius**2
F = phot_table_local_bkg['aperture_sum_1'][source_num] - (median_clipped * aperture_area)
sigma_F = np.sqrt((F / g) + (aperture_area * std_clipped ** 2))

print (phot_table_local_bkg['aperture_sum_1'][source_num],  sigma_F)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1021673362098894, 'fraction_numerical': 0.012101373803501464, 'mean_word_length': 5.35687263556116, 'pattern_counts': {'":': 0, '<': 122, '<?xml version=': 0, '>': 118, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '29', 'dir': 'jupyter-structured-clean-dedup', 'id': '290564', 'license': 'MIT', 'n_tokens_mistral': 6780, 'n_tokens_neox': 6391, 'n_words': 1776, 'path': 'notebooks/photutils.ipynb', 'repo_name': 'vilardellsalles/master-uab', 'size': '20163'}