Aperture Photometry JWST for Image Analyst

Exploring JWST image raw data from space

November 08, 2022 · 55 mins read

Aperture Photometry JWST for Image Analyst

Aperture photometry is the measurement of light which falls inside a particular aperture; usually, we mean a circular aperture of some fixed size. This project is dedicated to performing aperture photometry on identified sources of protostellar outflows in the NGC 3324 region of the Carina Nebula.

The formula for calculating the radius of the apertures is defined as,

png

Where $\lambda$ is the central wavelength and D is the diameter of the telescope (6.5m). This is twice the diffraction limit of each individual filter in arcseconds.

# use the entire screen width for the notebook
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:95% !important; }</style>"))

# make interactive plots
%matplotlib notebook

# Imports
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import glob # for filepath handling
import sys
#
import warnings
warnings.filterwarnings('ignore')
#
# Astropy:
from astropy.stats import SigmaClip
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.io import fits
from astropy.wcs import WCS
from astropy.visualization import simple_norm
from photutils.aperture import CircularAperture, CircularAnnulus, ApertureStats, aperture_photometry
import reproject as rpj # needed for reproject module

sys.path.append('../code')
from reproject_combine import reproject, make_rgb
from continuum_sub import ContinuumSubtract
import outflow_lookup as ol

path_to_file = '/mnt/c/Users/mukhamad-a/OneDrive - Playground/Data Science/outflow_sources_NGC3324.xlsx' # Set file path of the excel file
dict_ = ol.get_dictionary(path_to_file)

def sort_fnames(file_names):
    """ Function for sorting string lists of file names. Follows JWST naming conventions for naming .fits files."""
    # First get list of integers
    int_ind = []
    for f in file_names:
        num = f[-14:-10]

        if num[0] == 'f': # Accounts for MIRI f1130 band and up
            num = int(num[1:])
        else:
            pass

        int_ind.append(int(num))
    
    # Now sort the list of int returning indices that would sort the array
    sort_indices = np.argsort(int_ind)
    file_arr = np.array(file_names)
    sorted_fnames = file_arr[sort_indices] # re-arrange fnames
    
    return sorted_fnames # returns array of sorted fnames

# Get image directory and get files for imaging and plot generation
img_dir = '/mnt/d/st_images/Carina_level3/'
files = glob.glob(img_dir+'*')[3:] # Gets all JWST images including MIRI, excluding folders
files = sort_fnames(files) # sort them in ascending wavelength number
print('Files found:')
for f in files:
    print(f)

Files found:
/mnt/d/st_images/Carina_level3/jw02731-o001_t017_nircam_clear-f090w_i2d.fits
/mnt/d/st_images/Carina_level3/jw02731-o001_t017_nircam_clear-f187n_i2d.fits
/mnt/d/st_images/Carina_level3/jw02731-o001_t017_nircam_clear-f200w_i2d.fits
/mnt/d/st_images/Carina_level3/jw02731-o001_t017_nircam_clear-f335m_i2d.fits
/mnt/d/st_images/Carina_level3/jw02731-o001_t017_nircam_clear-f444w_i2d.fits
/mnt/d/st_images/Carina_level3/jw02731-o001_t017_nircam_f444w-f470n_i2d.fits
/mnt/d/st_images/Carina_level3/jw02731-o002_t017_miri_f770w_i2d.fits
/mnt/d/st_images/Carina_level3/jw02731-o002_t017_miri_f1130w_i2d.fits
/mnt/d/st_images/Carina_level3/jw02731-o002_t017_miri_f1280w_i2d.fits
/mnt/d/st_images/Carina_level3/jw02731-o002_t017_miri_f1800w_i2d.fits

# Set target header/data for reprojection
with fits.open(files[5]) as hdu:
    tar_header = hdu[1].header
    tar_data = hdu[1].data

target_wcs = WCS(tar_header)

# Reproject all image files onto target image (F444W-F470N)
reproj = reproject(files, target_wcs, tar_data)

We will now do a similar process to what was done to get the pixel values for aperture photometry. We will crop the reprojected images and then do the photometry on the products. Using sky coordinates would be optimal, but this batch of MIRI images from JWST had incorrect WCS header information and so the reprojected images are slightly offset from the NIRCam images. So we must result to using pixel values rather than coordinates for aperture setting. Also, we will work with individual sources rather than setting apertures on the full image.

source-1 (PJ001)

png

focus = 'PJ001' # set desired focus
coords = ol.get_coordinates(path_to_file, focus) # Get coordinates of the source from the excel file
pixel_x, pixel_y = target_wcs.world_to_pixel(coords) # get pixel coordinates of the outflow
pixel_x, pixel_y = int(pixel_x), int(pixel_y)
region = ol.pixel_region(pixel_x, pixel_y, custom_bounds=(200,200,300,300))

# Crop images onto desired region
cropped = ol.crop_images(reproj, region)

factor = ((0.063 * u.arcsec) ** 2) / u.pixel # pixel scale (area per pixel)
data = u.Quantity(cropped[-1], unit= u.MJy / u.sr) # F1800W MIRI image as data
converted_data = (data * factor).to(u.Jy / u.pixel) # convert from units of MJy/sr to Jy/pix

We will place an aperture on the theorised host star on the MIRI image along with an annulus and plot it below to confirm we have a correct radius set.

Other possible source of PJ001

The exact source of PJ001 is not clear. So we will also perform astro-photometry on the nearby star, another possible source.

pos = [(225, 257), (212, 237)] # set central positions of apertures
factor = 1 / 0.063 # plate scale = 0.063"/pixel
radius = 1.9 * (1.22) * (18e-6 / 6.5) * (180 / np.pi) * 3600 * factor # radius defined as twice the diffraction limit
apertures = CircularAperture(pos, r=radius)
annulus = CircularAnnulus(pos, r_in=23, r_out=30)

fig, axs = plt.subplots(1,1)
fig.set_size_inches(20, 16)

axs.imshow(cropped[-1], origin='lower', norm = simple_norm(cropped[-1], stretch='log', log_a=100, invalid=0, min_cut=170), cmap='hot')
axs.set_title(f'MIRI image of {focus}')

ap_patch = apertures.plot(color='w', lw=2, label='Photometry aperture')
ann_patch = annulus.plot(color='k', lw=2, label='Background annulus')
plt.legend(loc='best')
plt.tight_layout()

<IPython.core.display.Javascript object>


# Get stats on background, so we can subtract it. (Local Background Subtraction)
sigclip = SigmaClip(sigma=1, maxiters=10) # sigma clip due to bright star in annulus
bkg_stats = ApertureStats(converted_data.value, annulus, sigma_clip=sigclip)
aper_stats = ApertureStats(converted_data.value, apertures)
aper_stats_bkgsub = ApertureStats(converted_data.value, apertures, local_bkg=bkg_stats.median) # Take the background median away from data
bkg_median = bkg_stats.median
print('Median of apertures:')
print(bkg_median)

Median of apertures:
[1.77942881e-05 2.20291722e-05]

We now have aperture statistics both with and without background. Now a table of properties is made for the apertures above.

# Print a table with full list of properties
t = aper_stats_bkgsub.to_table()
for col in t.colnames:
    if col == 'sky_centroid':
        pass
    else:
        t[col].info.format = '%.8g'
t['id'][:] = [1, 1]
t.show_in_notebook()

QTable length=2

idxidxcentroidycentroidsky_centroidsumsum_errsum_aper_areacenter_aper_areaminmaxmeanmedianmodestdmad_stdvarbiweight_locationbiweight_midvariancefwhmsemimajor_sigmasemiminor_sigmaorientationeccentricity
pix2pix2pixpixpixdeg
01214.8561258.02183None0.061886952nan1387.599713735.8347791e-070.000953280774.3991127e-051.61833e-05-3.9432353e-050.000112809261.567327e-051.272593e-081.5640835e-052.5287746e-1019.3846169.81913116.2540099-14.8797230.77092915
11207.9717244.04664None0.033468442nan1387.59971373-1.9842662e-060.000676857732.3656491e-051.5493832e-05-8.3148628e-074.3875713e-051.8064624e-051.9250782e-091.6341404e-053.2103485e-1021.90449610.1528788.364976735.8765190.56673232
mask = apertures.to_mask()
values_1 = []
for i, _ in enumerate(apertures): 
    val = mask[i].get_values(converted_data.value) # get all values for the aperture if needed.
    values_1.append(val)

# Write array values to .txt file
fname = 'miri_aper_values_source1.txt'
arr = np.array(values_1[0])
out_file = open(fname, 'w')
for i in range(len(arr)):
    print(arr[i], file=out_file)

source-6 (PJ006)

png

focus = 'PJ006' # set desired focus
coords = ol.get_coordinates(path_to_file, focus) # Get coordinates of the source from the excel file
pixel_x, pixel_y = target_wcs.world_to_pixel(coords) # get pixel coordinates of the outflow
pixel_x, pixel_y = int(pixel_x), int(pixel_y)
region = ol.pixel_region(pixel_x, pixel_y, custom_bounds=(100,160,100,120))

# Crop images onto desired region
cropped = ol.crop_images(reproj, region)

factor = ((0.063 * u.arcsec) ** 2) / u.pixel # pixel scale (area per pixel)
data = u.Quantity(cropped[7], unit= u.MJy / u.sr) # F1130W MIRI image as data
converted_data = (data * factor).to(u.Jy / u.pixel) # convert from units of MJy/sr to Jy/pix

pos = (100, 97.5) # set central positions of apertures
factor = 1 / 0.063 # plate scale = 0.063"/pixel
radius = 2 * (1.22) * (11.3e-6 / 6.5) * (180 / np.pi) * 3600 * factor # radius defined as twice the diffraction limit
aperture = CircularAperture(pos, r=radius)
annulus = CircularAnnulus(pos, r_in=20, r_out=25)

fig, axs = plt.subplots(1,1)
fig.set_size_inches(20, 16)

axs.imshow(cropped[7], origin='lower', norm = simple_norm(cropped[7], stretch='log', log_a=50, invalid=0, min_cut=0), cmap='hot')
axs.set_title(f'MIRI image of {focus}')

ap_patch = aperture.plot(color='w', lw=2, label='Photometry aperture')
ann_patch = annulus.plot(color='k', lw=2, label='Background annulus')
plt.legend(loc='best')
plt.tight_layout()

<IPython.core.display.Javascript object>


# Get stats on background, so we can subtract it. (Local Background Subtraction)
bkg_stats = ApertureStats(converted_data.value, annulus)
aper_stats = ApertureStats(converted_data.value, aperture)
aper_stats_bkgsub = ApertureStats(converted_data.value, aperture, local_bkg=bkg_stats.median) # Take the background median away from data
bkg_median = bkg_stats.median
print('Median of apertures:')
print(bkg_median)

Median of apertures:
6.088006454125928e-06

# Print a table with full list of properties
t = aper_stats_bkgsub.to_table()
for col in t.colnames:
    if col == 'sky_centroid':
        pass
    else:
        t[col].info.format = '%.8g'
t['id'][:] = [6]
t.show_in_notebook()

QTable length=1

idxidxcentroidycentroidsky_centroidsumsum_errsum_aper_areacenter_aper_areaminmaxmeanmedianmodestdmad_stdvarbiweight_locationbiweight_midvariancefwhmsemimajor_sigmasemiminor_sigmaorientationeccentricity
pix2pix2pixpixpixdeg
06100.3681797.623202None0.0091899611nan605.93895608-1.0394475e-060.000158286851.5122451e-054.4280449e-06-1.6960767e-052.9661511e-052.8760425e-068.7980521e-103.9122153e-067.3454668e-128.99545833.94072363.695374716.1511080.34733763
mask = aperture.to_mask()
values_6 = mask.get_values(converted_data.value) # get all values for the aperture if needed.

source-7 (PJ007)

png

focus = 'PJ007' # set desired focus
coords = ol.get_coordinates(path_to_file, focus) # Get coordinates of the source from the excel file
pixel_x, pixel_y = target_wcs.world_to_pixel(coords) # get pixel coordinates of the outflow
pixel_x, pixel_y = int(pixel_x), int(pixel_y)
region = ol.pixel_region(pixel_x, pixel_y, custom_bounds=(100,100,100,110))

# Crop images onto desired region
cropped = ol.crop_images(reproj, region)

factor = ((0.063 * u.arcsec) ** 2) / u.pixel # pixel scale (area per pixel)
data = u.Quantity(cropped[6], unit= u.MJy / u.sr) # F770W MIRI image as data
converted_data = (data * factor).to(u.Jy / u.pixel) # convert from units of MJy/sr to Jy/pix

pos = (101, 97.5) # set central positions of apertures

factor = 1 / 0.063 # plate scale = 0.063"/pixel
radius = 2 * (1.22) * (7.7e-6 / 6.5) * (180 / np.pi) * 3600 * factor # radius defined as twice the diffraction limit

aperture = CircularAperture(pos, r=radius)
annulus = CircularAnnulus(pos, r_in=20, r_out=25)

fig, axs = plt.subplots(1,1)
fig.set_size_inches(20, 16)

axs.imshow(cropped[6], origin='lower', norm = simple_norm(cropped[6], stretch='log', invalid=0, min_cut=0), cmap='hot')
axs.set_title(f'MIRI image of {focus}')

ap_patch = aperture.plot(color='w', lw=2, label='Photometry aperture')
ann_patch = annulus.plot(color='k', lw=2, label='Background annulus')
plt.legend(loc='best')
plt.tight_layout()

<IPython.core.display.Javascript object>


# Get stats on background, so we can subtract it. (Local Background Subtraction)
bkg_stats = ApertureStats(converted_data.value, annulus)
aper_stats = ApertureStats(converted_data.value, aperture)
aper_stats_bkgsub = ApertureStats(converted_data.value, aperture, local_bkg=bkg_stats.median) # Take the background median away from data
bkg_median = bkg_stats.median
print('Median of apertures:')
print(bkg_median)

Median of apertures:
4.106795917915792e-06

# Print a table with full list of properties
t = aper_stats_bkgsub.to_table()
for col in t.colnames:
    if col == 'sky_centroid':
        pass
    else:
        t[col].info.format = '%.8g'
t['id'][:] = [7]
t.show_in_notebook()

QTable length=1

idxidxcentroidycentroidsky_centroidsumsum_errsum_aper_areacenter_aper_areaminmaxmeanmedianmodestdmad_stdvarbiweight_locationbiweight_midvariancefwhmsemimajor_sigmasemiminor_sigmaorientationeccentricity
pix2pix2pixpixpixdeg
0799.80214596.912526None0.00056846865nan281.35422278-9.3192261e-072.1047617e-052.0462707e-065.7888722e-07-2.3558797e-064.1924192e-061.128552e-061.7576379e-114.4809382e-071.2308537e-124.65199952.24912021.657357878.2565290.67601091
mask = aperture.to_mask()
values_7 = mask.get_values(converted_data.value) # get all values for the aperture if needed.

source-11

png

focus = 'PJ011' # set desired focus
coords = ol.get_coordinates(path_to_file, focus) # Get coordinates of the source from the excel file
pixel_x, pixel_y = target_wcs.world_to_pixel(coords) # get pixel coordinates of the outflow
pixel_x, pixel_y = int(pixel_x), int(pixel_y)
region = ol.pixel_region(pixel_x, pixel_y, custom_bounds=(20,40,100,100))

# Crop images onto desired region
cropped = ol.crop_images(reproj, region)

factor = ((0.063 * u.arcsec) ** 2) / u.pixel # pixel scale (area per pixel)
data = u.Quantity(cropped[5], unit= u.MJy / u.sr) # F444w-F470N NIRCam image as data
converted_data = (data * factor).to(u.Jy / u.pixel) # convert from units of MJy/sr to Jy/pix

pos = (20, 101) # set central positions of apertures
factor = 1 / 0.063 # plate scale = 0.063"/pixel
radius = 2 * (1.22) * (4.7e-6 / 6.5) * (180 / np.pi) * 3600 * factor # radius defined as twice the diffraction limit
aperture = CircularAperture(pos, r=7)
annulus = CircularAnnulus(pos, r_in=20, r_out=25)

fig, axs = plt.subplots(1,1)
fig.set_size_inches(20, 16)

axs.imshow(cropped[5], origin='lower', norm = simple_norm(cropped[5], stretch='log', log_a=100, invalid=0, min_cut=0), cmap='hot')
axs.set_title(f'NIRCam image of {focus}')

ap_patch = aperture.plot(color='w', lw=2, label='Photometry aperture')
ann_patch = annulus.plot(color='k', lw=2, label='Background annulus')
plt.legend(loc='best')
plt.tight_layout()

<IPython.core.display.Javascript object>


# Get stats on background, so we can subtract it. (Local Background Subtraction)
bkg_stats = ApertureStats(converted_data.value, annulus)
aper_stats = ApertureStats(converted_data.value, aperture)
aper_stats_bkgsub = ApertureStats(converted_data.value, aperture, local_bkg=bkg_stats.median) # Take the background median away from data
bkg_median = bkg_stats.median
print('Median of apertures:')
print(bkg_median)

Median of apertures:
3.5729745461913834e-07

# Print a table with full list of properties
t = aper_stats_bkgsub.to_table()
for col in t.colnames:
    if col == 'sky_centroid':
        pass
    else:
        t[col].info.format = '%.8g'
t['id'][:] = [11]
t.show_in_notebook()

QTable length=1

idxidxcentroidycentroidsky_centroidsumsum_errsum_aper_areacenter_aper_areaminmaxmeanmedianmodestdmad_stdvarbiweight_locationbiweight_midvariancefwhmsemimajor_sigmasemiminor_sigmaorientationeccentricity
pix2pix2pixpixpixdeg
01120.107608100.23149None4.7459365e-05nan153.938041455.1410348e-082.2184838e-063.222854e-071.58484e-07-1.691188e-074.1859386e-071.1339579e-071.7522082e-131.5561315e-071.4231125e-145.66755052.63161132.158673277.8630460.57195397
mask = aperture.to_mask()
values_11 = mask.get_values(converted_data.value) # get all values for the aperture if needed.

source-12 (PJ012)

png

We will use a NIRCam image to perform photometry on the star as the theorised source is unresolved in the MIRI image.

focus = 'PJ012' # set desired focus
coords = ol.get_coordinates(path_to_file, focus) # Get coordinates of the source from the excel file
pixel_x, pixel_y = target_wcs.world_to_pixel(coords) # get pixel coordinates of the outflow
pixel_x, pixel_y = int(pixel_x), int(pixel_y)
region = ol.pixel_region(pixel_x, pixel_y, custom_bounds=(200,130,70,200))

# Crop images onto desired region
cropped = ol.crop_images(reproj, region)

factor = ((0.063 * u.arcsec) ** 2) / u.pixel # pixel scale (area per pixel)
data = u.Quantity(cropped[4], unit= u.MJy / u.sr) # F444w NIRCam image as data
converted_data = (data * factor).to(u.Jy / u.pixel) # convert from units of MJy/sr to Jy/pix

pos = (200.7, 69.5) # set central positions of apertures
factor = 1 / 0.063 # plate scale = 0.063"/pixel
radius = 2 * (1.22) * (4.44e-6 / 6.5) * (180 / np.pi) * 3600 * factor # radius defined as twice the diffraction limit
aperture = CircularAperture(pos, r=radius)
annulus = CircularAnnulus(pos, r_in=10, r_out=15)

fig, axs = plt.subplots(1,1)
fig.set_size_inches(20, 16)

axs.imshow(cropped[4], origin='lower', norm = simple_norm(cropped[4], stretch='log', log_a=100, invalid=0, min_cut=0), cmap='hot')
axs.set_title(f'NIRCam image of {focus}')

ap_patch = aperture.plot(color='w', lw=2, label='Photometry aperture')
ann_patch = annulus.plot(color='k', lw=2, label='Background annulus')
plt.legend(loc='best')
plt.tight_layout()

<IPython.core.display.Javascript object>


# Get stats on background, so we can subtract it. (Local Background Subtraction)
bkg_stats = ApertureStats(converted_data.value, annulus)
aper_stats = ApertureStats(converted_data.value, aperture)
aper_stats_bkgsub = ApertureStats(converted_data.value, aperture, local_bkg=bkg_stats.median) # Take the background median away from data
bkg_median = bkg_stats.median
print('Median of apertures:')
print(bkg_median)

Median of apertures:
4.1188862690205246e-07

# Print a table with full list of properties
t = aper_stats_bkgsub.to_table()
for col in t.colnames:
    if col == 'sky_centroid':
        pass
    else:
        t[col].info.format = '%.8g'
t['id'][:] = [12]
t.show_in_notebook()

QTable length=1

idxidxcentroidycentroidsky_centroidsumsum_errsum_aper_areacenter_aper_areaminmaxmeanmedianmodestdmad_stdvarbiweight_locationbiweight_midvariancefwhmsemimajor_sigmasemiminor_sigmaorientationeccentricity
pix2pix2pixpixpixdeg
012200.1574168.435335None9.528288e-05nan93.54873590-2.6521996e-094.7061433e-061.0252579e-065.3024994e-07-4.5976592e-071.1482339e-066.0760922e-071.3184411e-125.410148e-076.2015472e-134.75768812.31389231.676295857.0380750.68932893
mask = aperture.to_mask()
values_12 = mask.get_values(converted_data.value) # get all values for the aperture if needed.

source-19 (PJ019)

png

We will once again use a NIRCam image to perform photometry on the star for similar reasons listed above.

focus = 'PJ019' # set desired focus
coords = ol.get_coordinates(path_to_file, focus) # Get coordinates of the source from the excel file
pixel_x, pixel_y = target_wcs.world_to_pixel(coords) # get pixel coordinates of the outflow
pixel_x, pixel_y = int(pixel_x), int(pixel_y)
region = ol.pixel_region(pixel_x, pixel_y, custom_bounds=(180, 180, 230, 120))

# Crop images onto desired region
cropped = ol.crop_images(reproj, region)

factor = ((0.063 * u.arcsec) ** 2) / u.pixel # pixel scale (area per pixel)
data = u.Quantity(cropped[5], unit= u.MJy / u.sr) # F444w-F470N NIRCam image as data
converted_data = (data * factor).to(u.Jy / u.pixel) # convert from units of MJy/sr to Jy/pix

pos = (181, 230) # set central positions of apertures
factor = 1 / 0.063 # plate scale = 0.063"/pixel
radius = 2 * (1.22) * (4.7e-6 / 6.5) * (180 / np.pi) * 3600 * factor # radius defined as twice the diffraction limit
aperture = CircularAperture(pos, r=radius)
annulus = CircularAnnulus(pos, r_in=15, r_out=20)

fig, axs = plt.subplots(1,1)
fig.set_size_inches(20, 16)

axs.imshow(cropped[5], origin='lower', norm = simple_norm(cropped[5], stretch='log', log_a=100, invalid=0, min_cut=0), cmap='hot')
axs.set_title(f'NIRCam image of {focus}')

ap_patch = aperture.plot(color='w', lw=2, label='Photometry aperture')
ann_patch = annulus.plot(color='k', lw=2, label='Background annulus')
plt.legend(loc='best')
plt.tight_layout()

<IPython.core.display.Javascript object>


# Get stats on background, so we can subtract it. (Local Background Subtraction)
bkg_stats = ApertureStats(converted_data.value, annulus)
aper_stats = ApertureStats(converted_data.value, aperture)
aper_stats_bkgsub = ApertureStats(converted_data.value, aperture, local_bkg=bkg_stats.median) # Take the background median away from data
bkg_median = bkg_stats.median
print('Median of apertures:')
print(bkg_median)

Median of apertures:
5.91098810624582e-07

# Print a table with full list of properties
t = aper_stats_bkgsub.to_table()
for col in t.colnames:
    if col == 'sky_centroid':
        pass
    else:
        t[col].info.format = '%.8g'
t['id'][:] = [19]
t.show_in_notebook()

QTable length=1

idxidxcentroidycentroidsky_centroidsumsum_errsum_aper_areacenter_aper_areaminmaxmeanmedianmodestdmad_stdvarbiweight_locationbiweight_midvariancefwhmsemimajor_sigmasemiminor_sigmaorientationeccentricity
pix2pix2pixpixpixdeg
019180.98546229.95316None0.0007261677nan104.825681016.5229762e-075.7776108e-057.1617053e-062.8176777e-06-5.8703776e-061.1049918e-051.6160435e-061.2210068e-102.5580423e-061.7973869e-124.19803411.82721871.7371248-28.4756830.3101322
mask = aperture.to_mask()
values_19 = mask.get_values(converted_data.value) # get all values for the aperture if needed.

source-21 (PJ021)

focus = 'PJ021' # set desired focus
coords = ol.get_coordinates(path_to_file, focus) # Get coordinates of the source from the excel file
pixel_x, pixel_y = target_wcs.world_to_pixel(coords) # get pixel coordinates of the outflow
pixel_x, pixel_y = int(pixel_x), int(pixel_y)
region = ol.pixel_region(pixel_x, pixel_y, custom_bounds=(50, 50, 80, 50))

# Crop images onto desired region
cropped = ol.crop_images(reproj, region)

factor = ((0.063 * u.arcsec) ** 2) / u.pixel # pixel scale (area per pixel)
data = u.Quantity(cropped[5], unit= u.MJy / u.sr) # F444w-F470N NIRCam image as data
converted_data = (data * factor).to(u.Jy / u.pixel) # convert from units of MJy/sr to Jy/pix

pos = (50, 80) # set central positions of apertures
factor = 1 / 0.063 # plate scale = 0.063"/pixel
radius = 2 * (1.22) * (4.7e-6 / 6.5) * (180 / np.pi) * 3600 * factor # radius defined as twice the diffraction limit
aperture = CircularAperture(pos, r=radius)
annulus = CircularAnnulus(pos, r_in=15, r_out=20)

fig, axs = plt.subplots(1,1)
fig.set_size_inches(20, 16)

axs.imshow(cropped[5], origin='lower', norm = simple_norm(cropped[5], stretch='log', log_a=100, invalid=0, min_cut=0), cmap='hot')
axs.set_title(f'NIRCam image of {focus}')

ap_patch = aperture.plot(color='w', lw=2, label='Photometry aperture')
ann_patch = annulus.plot(color='k', lw=2, label='Background annulus')
plt.legend(loc='best')
plt.tight_layout()

<IPython.core.display.Javascript object>


# Get stats on background, so we can subtract it. (Local Background Subtraction)
bkg_stats = ApertureStats(converted_data.value, annulus)
aper_stats = ApertureStats(converted_data.value, aperture)
aper_stats_bkgsub = ApertureStats(converted_data.value, aperture, local_bkg=bkg_stats.median) # Take the background median away from data
bkg_median = bkg_stats.median
print('Median of apertures:')
print(bkg_median)

Median of apertures:
5.63654412687679e-07

# Print a table with full list of properties
t = aper_stats_bkgsub.to_table()
for col in t.colnames:
    if col == 'sky_centroid':
        pass
    else:
        t[col].info.format = '%.8g'
t['id'][:] = [21]
t.show_in_notebook()

QTable length=1

idxidxcentroidycentroidsky_centroidsumsum_errsum_aper_areacenter_aper_areaminmaxmeanmedianmodestdmad_stdvarbiweight_locationbiweight_midvariancefwhmsemimajor_sigmasemiminor_sigmaorientationeccentricity
pix2pix2pixpixpixdeg
02150.10501480.078556None0.0010612818nan104.825681019.6577569e-078.1531731e-051.0482058e-053.517837e-06-1.0410605e-051.6501911e-052.4533412e-062.7231306e-103.1331458e-064.7657673e-124.03061251.79551931.6234417.058780.42719015
mask = aperture.to_mask()
values_21 = mask.get_values(converted_data.value) # get all values for the aperture if needed.
← Previous Post Next Post →