Run Image Pipeline and Create Source Catalog#
In this run, we are starting from the rate files which have already been calibrated with the detector1 pipeline. This will save time as we do not need to edit any of the steps being performed as part of detector1. Therefore, the first calibration that should be done as part of a WFSS run is to run the rate files direct images through the Image2 and Image3 steps of the JWST pipeline. This includes creating a source catalog, which most likely will need to be adjusted from the pipeline default values. Not having a good source catalog will result in non optimal extraction of sources in the dispersed, WFSS, images.
Use case: The default parameters for the pipeline do not extract the expected sources, so custom parameters need to be set to obtain new combined image and source catalog.
Data: JWST/NIRISS images and spectra from program 2079 observation 004. This should be stored in a single directory data, and can be downloaded from the previous notebook, 00_niriss_mast_query_data_setup.ipynb.
Tools: astropy, crds, glob, jdaviz, json, jwst, matplotlib, numpy, os, pandas, urllib, warnings, zipfile
Cross-instrument: NIRISS
Content
Imports & Data Setup
Default Imaging Pipeline Run
Image2
Image3
Inspecting Default Results
Custom Imaging Pipeline Run
Image3
Inspecting Custom Results
Refining the Source Catalog Further
Matching Source IDs Across Catalogs
Manually Editing the Source Catalog
Author: Rachel Plesha (rplesha@stsci.edu), Camilla Pacifici (cpacifici@stsci.edu), JWebbinar notebooks.
First Published: May 2024
Last tested: This notebook was last tested with JWST pipeline version 1.12.5 and the CRDS context jwst_1225.pmap.
Imports & Data Setup#
# Update the CRDS path to your local directory
%env CRDS_PATH=crds_cache
%env CRDS_SERVER_URL=https://jwst-crds.stsci.edu
import os
import glob
import json
import warnings
import urllib
import zipfile
import numpy as np
import pandas as pd
import astropy.units as u
from astropy.io import fits
from astropy.coordinates import SkyCoord
from astropy.table import Table
from jdaviz import Imviz
from matplotlib import pyplot as plt
%matplotlib widget
# %matplotlib inline
from jwst.pipeline import Image2Pipeline
from jwst.pipeline import Image3Pipeline
Check what version of the JWST pipeline you are using. To see what the latest version of the pipeline is available or install a previous version, check GitHub. Also verify what CRDS version you are using. CRDS documentation explains how to set a specific context to use in the JWST pipeline. If either of these values are different from the last tested note above there may be differences in this notebook.
import jwst
import crds
print('JWST Pipeliene Version:', jwst.__version__)
print('CRDS Context:', crds.get_context_name('jwst'))
The data directory, data_dir here should contain all of the association and rate files in a single, flat directory. default_run_image3 and custom_run_image3 are directories that we will use later for our calibrated data. They are separated so that we can compare the two outputs.
data_dir = 'data'
default_run_image3 = 'default_image3_calibrated' # where the results of the default image3 run will be saved (inside of data_dir)
custom_run_image3 = 'custom_image3_calibrated'# where the results of the custom image3 run will be saved (inside of data_dir)
The association files expect that 1) all of the data are in the same directory and 2) that you are performing the pipeline call also in that directory. Because of that, we need to change into the data directory to run the imaging pipelines.
# if you have not downloaded the data from notebook 00, run this cell. Otherwise, feel free to skip it.
# Download uncalibrated data from Box into the data directory:
boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/niriss_wfss_advanced/niriss_wfss_advanced_01_input.zip'
boxfile = os.path.basename(boxlink)
urllib.request.urlretrieve(boxlink, boxfile)
zf = zipfile.ZipFile(boxfile, 'r')
zf.extractall(path=data_dir)
# move the files downloaded from the box file into the top level data directory
box_download_dir = os.path.join(data_dir, boxfile.split('.zip')[0])
for filename in glob.glob(os.path.join(box_download_dir, '*')):
if '.csv' in filename:
# move to the current directory
os.rename(filename, os.path.basename(filename))
else:
# move to the data directory
os.rename(filename, os.path.join(data_dir, os.path.basename(filename)))
# remove unnecessary files now
os.remove(boxfile)
os.rmdir(box_download_dir)
# From the csv file we created earlier, find a list of all of the grism observations we will want to calibrate with spec2
listrate_file = 'list_ngdeep_rate.csv'
rate_df = pd.read_csv(listrate_file)
cwd = os.getcwd() # get the current working directory
if cwd != data_dir: # if you are not already in the location of the data, change into it
try:
os.chdir(data_dir)
except FileNotFoundError:
print(f"Not able to change into: {data_dir}.\nRemaining in: {cwd}")
pass
for temp_dir in [default_run_image3, custom_run_image3]:
if not os.path.exists(temp_dir):
os.mkdir(temp_dir)
Default Imaging Pipeline Run#
To start, run the default image2 and image3 steps of the pipeline on all direct images observed with the WFSS data.
Run Default Image2#
Image2 is run on the direct image rate files. While your program should have valid association files to download from MAST, if for any reason you need to make your own association file, see Creating Custom ASN Files.
Looking in a Level 2 Imaging Association File#
First, take a look inside the association (ASN) files to better understand everything that is contained in them.
For image2 association files, there should be one asn file for each dither position in an observing sequence which is set by the exposure strategy. In this case, that should match the number of direct images (FILTER=CLEAR) in rate_df because each direct image is at a unique dither position (XOFFSET, YOFFSET) within an observing sequence. For this program and observation, there is one direct image before a grism change with only one dither, another direct image with four dithers between the change in grisms, and a direct image at the end of a sequence with three dithers. This leads to a total of eight images per observing sequence, with five observing sequences in the observation using the blocking filters F115W -> F115W -> F150W -> F150W -> F200W.
image2_asns = glob.glob('*image2*asn*.json')
print(len(image2_asns), 'Image2 ASN files') # there should be 40 asn files for image2
# the number of association files should match the number of rate files
print(len(rate_df[rate_df['FILTER'] == 'CLEAR']), 'Direct Image rate files')
# look at one of the association files
asn_data = json.load(open(image2_asns[0]))
for key, data in asn_data.items():
print(f"{key} : {data}")
From this association, we can tell many things about the observation:
From
asn_typeandasn_rule, we can see that this is an image2 associationFrom
degraded_statuswe can see that there are no exposures to not be included in the calibration.From
constraints, we can see this is not a time series observation (TSO), the observation is part of program 2079, observed with NIRISS with the CLEAR (i.e. imaging for WFSS) and F150W filters.From
productswe can see there is only one exposure associated. This is typical for image2 where there is usually only one exposure per dither per observing sequence.
# in particular, take a closer look at the product filenames with the association file:
for product in asn_data['products']:
for key, value in product.items():
if key == 'members':
print(f"{key}:")
for member in value:
print(f" {member['expname']} {member['exptype']}")
else:
print(f"{key}: {value}")
Run image2#
The rate.fits products will be calibrated into cal.fits files. More information about the steps performed in the Image2 part of the pipeline can be found in the Image2 pipeline documentation.
In this case, we’re saving the outputs to the same directory we are running the pipeline in so that we can then use the output cal files to run the Image3 pipeline
for img2_asn in image2_asns:
# check if the calibrated file already exists
asn_data = json.load(open(img2_asn))
cal_file = f"{asn_data['products'][0]['name']}_cal.fits"
if os.path.exists(cal_file):
print(cal_file, 'cal file already exists.')
continue
# if not, calibrated with image2
img2 = Image2Pipeline.call(img2_asn, save_results=True)
Run Default Image3#
Looking in a Level 3 Association File#
The contents are quite similar to image2, but notice now that there are many more members that are associated together, and they use the individual pointing cal files from image2. Image3 resamples and combines images of the same blocking filter (PUPIL for NIRISS WFSS) from all dither and observing sequences to form a single image, which leads to fewer image3 association files.
image3_asns = glob.glob('*image3*asn*.json')
print(len(image3_asns), 'Image3 ASN files') # there should be 3 image3 association files
# the number of image3 association files should match the number of unique blocking filters used
uniq_filters = np.unique(rate_df[rate_df['FILTER'] == 'CLEAR']['PUPIL'])
print(f"{len(uniq_filters)} unique filters used: {uniq_filters}")
# look at one of the association files
image3_asn_data = json.load(open(image3_asns[0]))
for key, data in image3_asn_data.items():
print(f"{key} : {data}")
# in particular, take a closer look at the product filenames with the association file:
for product in image3_asn_data['products']:
for key, value in product.items():
if key == 'members':
print(f"{key}:")
for member in value:
print(f" {member['expname']} {member['exptype']}")
else:
print(f"{key}: {value}")
Run image3#
In Image3, the cal.fits individual pointing files will be calibrated into a single combined i2d.fits image. The Image3 step is also where the final source catalog is created, so we can change some input paramters to obtain a more refined output source catalog. This is done below in the Custom Imaging Pipeline Run section. However, we will first calibrate the data with the default parameters. More information about the steps performed in the Image3 part of the pipeline can be found in the Image3 pipeline documentation.
Note: Image3 can take a while to run
for img3_asn in image3_asns:
# check if the calibrated file already exists
asn_data = json.load(open(img3_asn))
cal_file = os.path.join(default_run_image3, f"{asn_data['products'][0]['name']}_i2d.fits")
if os.path.exists(cal_file):
print(cal_file, 'cal file already exists.')
continue
# if not, calibrated with image3
img3 = Image3Pipeline.call(img3_asn, save_results=True, output_dir=default_run_image3)
# remove unnecessary files to save disk space
for crf in glob.glob(os.path.join(default_run_image3, '*crf.fits')):
os.remove(crf)
Inspecting Default Results#
# These are all resuts from the Image3 pipeline
image3_i2d = np.sort(glob.glob(os.path.join(default_run_image3, '*i2d.fits'))) # combined image over multiple dithers/mosaic
image3_segm = np.sort(glob.glob(os.path.join(default_run_image3, '*segm.fits'))) # segmentation map that defines the extent of a source
image3_cat = np.sort(glob.glob(os.path.join(default_run_image3, '*cat.ecsv'))) # Source catalog that defines the RA/Dec of a source at a particular pixel
Matplotlib#
Matplotlib has limitations where ImViz might better suite your needs – especially if you like to look at things in WCS coordinates. For the notebook purposes, we are highlighting a few key areas using the matplotlib package instead.
Using the i2d combined image and the source catalog produced by Image3, we can visually inspect if we’re happy with where the pipeline found the sources. In the following figures, what has been defined as an extended source by the pipeline is shown in orange-red, and what has been defined as a point source by the pipeline is shown in grey. This definition affects the extraction box in the WFSS images.
fig = plt.figure(figsize=(10, 10))
cols = 2
rows = int(np.ceil(len(image3_i2d) / cols))
for plt_num, img in enumerate(image3_i2d):
# determine where the subplot should be
xpos = (plt_num % 40) % cols
ypos = ((plt_num % 40) // cols) # // to make it an int.
# make the subplot
ax = plt.subplot2grid((rows, cols), (ypos, xpos))
# plot the image
with fits.open(img) as hdu:
ax.imshow(hdu[1].data, vmin=0, vmax=0.3, origin='lower')
ax.set_title(f"obs{hdu[0].header['OBSERVTN']} {hdu[0].header['PUPIL']}")
# also plot the associated catalog
cat = Table.read(img.replace('i2d.fits', 'cat.ecsv'))
extended_sources = cat[cat['is_extended'] == 1] # 1 is True; i.e. is extended
point_sources = cat[cat['is_extended'] == 0] # 0 is False; i.e. is a point source
ax.scatter(extended_sources['xcentroid'], extended_sources['ycentroid'], s=20, facecolors='None', edgecolors='orangered', alpha=0.9)
ax.scatter(point_sources['xcentroid'], point_sources['ycentroid'], s=20, facecolors='None', edgecolors='dimgrey', alpha=0.9)
# Helps to make the axes not overlap ; you can also set this manually if this doesn't work
plt.tight_layout()
The segmentation maps are also a product of the Image3 pipeline, and they are used the help determine the source catalog. Let’s take a look at those to ensure we are happy with what it is defining as a source.
In the segmentation map, each blue blob should correspond to a physical target. There are cases where sources can be blended, in which case the parameters for making the semgentation map and source catalog should be changed. An example of this can be seen below in the observation 004 F200W filter image where two galaxies at ~(1600, 1300) have been blended into one source. This is discussed in more detail below in Custom Imaging Pipeline Run.
Note that because of the filter offset difference the field of views for the shown in the cutouts below differs for each filter.
# we will look at this multiple times, so let's define this as a function
def plot_image_and_segmentation_map(i2d_images, segm_images, xmin=1250, xmax=1750, ymin=1250, ymax=1750, cat_suffix='cat.ecsv'):
cols = 2
rows = len(i2d_images)
fig = plt.figure(figsize=(10, 10*(rows/2)))
for plt_num, img in enumerate(np.sort(np.concatenate([segm_images, i2d_images]))):
# determine where the subplot should be
xpos = (plt_num % 40) % cols
ypos = ((plt_num % 40) // cols) # // to make it an int.
# make the subplot
ax = plt.subplot2grid((rows, cols), (ypos, xpos))
if 'i2d' in img:
cat = Table.read(img.replace('i2d.fits', cat_suffix))
cmap = 'gist_gray'
else:
cmap = 'tab20c_r'
# plot the image
with fits.open(img) as hdu:
ax.imshow(hdu[1].data, vmin=0, vmax=0.3, origin='lower', cmap=cmap)
title = f"obs{hdu[0].header['OBSERVTN']} {hdu[0].header['PUPIL']}"
# also plot the associated catalog
extended_sources = cat[cat['is_extended'] == 1] # 1 is True; i.e. is extended
point_sources = cat[cat['is_extended'] == 0] # 0 is False; i.e. is a point source
for color, sources in zip(['darkred', 'black'], [extended_sources, point_sources]):
# plotting the sources
ax.scatter(sources['xcentroid'], sources['ycentroid'], s=20, facecolors='None', edgecolors=color, alpha=0.9)
# adding source labels
for i, source_num in enumerate(sources['label']):
ax.annotate(source_num,
(sources['xcentroid'][i]+0.5, sources['ycentroid'][i]+0.5),
fontsize=8,
color=color)
if 'i2d' in img:
ax.set_title(f"{title} combined image")
else:
ax.set_title(f"{title} segmentation map")
# zooming in on a smaller region
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
# Helps to make the axes not overlap ; you can also set this manually if this doesn't work
plt.tight_layout()
return fig
default_fig = plot_image_and_segmentation_map(image3_i2d, image3_segm)
ImViz#
Similarly to DS9, ImViz allows you to interactively view these images and the corresponding source catalog as well.
imviz = Imviz()
viewer = imviz.default_viewer
# plot each i2d image
catalogs = [] # for plotting the catalogs
labels = [] # for plotting the catalogs
for img in image3_i2d:
print(f'Plotting: {img}')
label = f"obs{fits.getval(img, 'OBSERVTN')} {fits.getval(img, 'PUPIL')}"
with warnings.catch_warnings():
warnings.simplefilter('ignore')
imviz.load_data(img, data_label=label)
# save info to plot the catalogs next
catalogs.append(img.replace('i2d.fits', 'cat.ecsv'))
labels.append(label)
# this aligns the image to use the WCS coordinates;
# the images need to be loaded first, but before adding markers
imviz.link_data(link_type='wcs')
# also plot the associated catalog
# this needs to be a separate loop due to linking in imviz when using sky coordinates
for catname, label in zip(catalogs, labels):
cat = Table.read(catname)
# format the table into the format imviz expects
sky_coords = Table({'coord': [SkyCoord(ra=cat['sky_centroid'].ra.degree,
dec=cat['sky_centroid'].dec.degree,
unit="deg")]})
viewer.marker = {'color': 'orange', 'alpha': 1, 'markersize': 20, 'fill': False}
viewer.add_markers(sky_coords, use_skycoord=True, marker_name=f"{label} catalog")
# This changes the stretch of all of the images
plotopt = imviz.plugins['Plot Options']
plotopt.select_all(viewers=True, layers=True)
plotopt.stretch_preset = '99.5%'
imviz.show()
Custom Imaging Pipeline Run#
Image3#
Try editing a few parameters and compare the outcomes to the default run above, at first for a single file.
When we call the image3 pipeline, we can add modifications to a specific step of the pipeline. In this case we’re going to edit the source_catalog and tweakreg steps of the pipeline. An explanation of the different parameters to tweak can be found in the further information below, while the default values are a combination of both the default pipeline values listed in there, and the parameter reference files that are supplied.
image3_asns = np.sort(glob.glob('*image3*asn*.json'))
test_asn = image3_asns[1]
# check if the calibrated file already exists
asn_data = json.load(open(test_asn))
i2d_file = os.path.join(custom_run_image3, f"{asn_data['products'][0]['name']}_i2d.fits")
if os.path.exists(i2d_file):
print(i2d_file, 'i2d file already exists.')
else:
# call the image3 pipeline in the same way as before, but add a few new modifications
cust_img3 = Image3Pipeline.call(test_asn,
steps={
'source_catalog': {'kernel_fwhm': 5.0,
'snr_threshold': 10.0,
'npixels': 50,
'deblend': True,
},
'tweakreg': {'snr_threshold': 20,
'abs_refcat': 'GAIADR3',
'save_catalogs': True,
'searchrad': 3.0,
'kernel_fwhm': 2.302,
'fitgeometry': 'shift',
},
},
save_results=True,
output_dir=custom_run_image3)
# remove unnecessary files to save disk space
for crf in glob.glob(os.path.join(custom_run_image3, '*crf.fits')):
os.remove(crf)
Inspecting Custom Results#
default_i2d = os.path.join(default_run_image3, os.path.basename(i2d_file))
compare_i2ds = [i2d_file, default_i2d]
compare_segm = [i2d_file.replace('i2d.fits', 'segm.fits'), default_i2d.replace('i2d.fits', 'segm.fits')]
compare_fig = plot_image_and_segmentation_map(compare_i2ds, compare_segm)
The cell below shows similar information using Imviz instead to visualize this.
imviz = Imviz()
viewer = imviz.default_viewer
for img, label in zip([i2d_file, os.path.join(default_run_image3, os.path.basename(i2d_file))], ['custom', 'default']):
print(f'Plotting: {img}')
title = f"{label} obs{fits.getval(img, 'OBSERVTN')} {fits.getval(img, 'PUPIL')}"
with warnings.catch_warnings():
warnings.simplefilter('ignore')
imviz.load_data(img, data_label=title)
# this aligns the image to use the WCS coordinates
linking = imviz.plugins['Links Control']
linking.link_type = 'WCS'
# also plot the associated catalog
cat = Table.read(img.replace('i2d.fits', 'cat.ecsv'))
# format the table into the format imviz expects
t_xy = Table({'x': cat['xcentroid'],
'y': cat['ycentroid']})
viewer.marker = {'color': 'orange', 'alpha': 1, 'markersize': 20, 'fill': False}
viewer.add_markers(t_xy, marker_name=f"{label} catalog")
# This changes the stretch of all of the images
plotopt = imviz.plugins['Plot Options']
plotopt.select_all(viewers=True, layers=True)
plotopt.stretch_preset = '99.5%'
imviz.show()
Calibrate the remaining images if you are happy with the above results
image3_asns = np.sort(glob.glob('*image3*asn*.json'))
for img3_asn in image3_asns:
# check if the calibrated file already exists
asn_data = json.load(open(img3_asn))
i2d_file = os.path.join(custom_run_image3, f"{asn_data['products'][0]['name']}_i2d.fits")
if os.path.exists(i2d_file):
print(i2d_file, 'cal file already exists.')
continue
# call the image3 pipeline in the same way as before, but add a few new modifications
cust_img3 = Image3Pipeline.call(img3_asn,
steps={
'source_catalog': {'kernel_fwhm': 5.0,
'snr_threshold': 10.0,
'npixels': 50,
'deblend': True,
},
'tweakreg': {'snr_threshold': 20,
'abs_refcat': 'GAIADR3',
'save_catalogs': True,
'searchrad': 3.0,
'kernel_fwhm': 2.302,
'fitgeometry': 'shift',
},
},
save_results=True,
output_dir=custom_run_image3)
# remove unnecessary files to save disk space
for crf in glob.glob(os.path.join(custom_run_image3, '*crf.fits')):
os.remove(crf)
# These are all resuts from the Image3 pipeline
cust_image3_i2d = np.sort(glob.glob(os.path.join(custom_run_image3, '*i2d.fits'))) # combined image over multiple dithers/mosaic
cust_image3_segm = np.sort(glob.glob(os.path.join(custom_run_image3, '*segm.fits'))) # segmentation map that defines the extent of a source
custom_fig = plot_image_and_segmentation_map(cust_image3_i2d, cust_image3_segm)
Refining the Source Catalog Further#
In the above cases, we have modified the results using the pipeline directly. It might be the case that perhaps you want to use someone else’s custom source catalog, or modify the source catalog even further from what was output by the pipeline. In these cases, we will then need to modify the spec2 ASN files to point to the new source catalog, which will be discussed in the spec2 notebook. Additionally, it can be useful to match all of the source IDs across the different catalogs. In this case, there are three different catalogs created by the pipeline that identify the same RA/Dec or X/Y pixel location as a different source ID, so we will edit those to have the same source ID values across the catalogs. These extra steps aren’t always necessary, but could be helpful in analyses of NIRISS WFSS data.
Matching Source IDs Across Catalogs#
In the above figures, you can see that the same source has multiple source IDs. Here we want to match all of the source IDs across all observations to be sure we are talking about the same source regardless of which filter or observation we look at. To do this, we use the astropy match_to_catalog_3d function and rebase the labels.
The first step is to decide on a base catalog to match all of the other catalogs to. Here we’ll use the the observation 004 F115W filter.
custom_cats = np.sort(glob.glob(os.path.join(custom_run_image3, '*niriss_clear-f????_cat.ecsv'))) # cat filename format from image3 results
print("All image3 catalogs:\n", custom_cats)
base_cat = Table.read(custom_cats[0])
# save out the base catalog with a new name to be consistent
base_cat_name = custom_cats[0].replace('cat.ecsv', 'source-match_cat.ecsv')
base_cat.write(base_cat_name, overwrite=True)
print("\nBase catalog:", base_cat_name)
Loop through the remaining catalogs to match the IDs based off of sky coordinates matching to within 1 arcsecond.
max_sep = 1 * u.arcsec # adjust if necessary
base_sky = base_cat['sky_centroid']
for to_match_cat in custom_cats[1:]:
# read in the catalog
other_cat = Table.read(to_match_cat)
other_sky = other_cat['sky_centroid']
# find the matching sources between the two catalogs based on sky coordinates
idx, d2d, d3d = base_sky.match_to_catalog_3d(other_sky)
sep_constraint = d2d < max_sep
base_matches = base_cat[sep_constraint]
other_matches = other_cat[idx[sep_constraint]]
# rebase the ID values to be the same
other_matches['label'] = base_matches['label']
# save out the new catalog
match_cat_name = to_match_cat.replace('cat.ecsv', 'source-match_cat.ecsv')
other_matches.write(match_cat_name, overwrite=True)
print('Saved:', match_cat_name)
Look at the new source label numbers. They should match!
new_cat_fig = plot_image_and_segmentation_map(cust_image3_i2d, cust_image3_segm, cat_suffix='source-match_cat.ecsv',
xmin=1500, xmax=2000, ymin=800, ymax=1300)
Manually Editing the Source Catalog#
Looking ahead to the WFSS extraction, it might be that we only really care about one or two sources at any particular moment. In this case, it’s helpful to pair down the source catalog to just that source to make it easier to look at the extracted 1-D spectrum in the WFSS data.
To start, look at the current custom source catalog for one of the filters to get an idea of what is contained in the catalogs.
# first, look at the current, custom source catalog for the F200W filter
cat_name = np.sort(glob.glob(os.path.join(custom_run_image3, '*source-match_cat.ecsv')))[2]
cat = Table.read(cat_name)
print(cat_name)
cat
There may be multiple ways to look for a source, so shown below are three options:
With a known RA/Dec of an object
A known x/y location of an object
With a source ID of an object. Note that we are using the rebased source catalogs here for the IDs.
This same concept can be extended to filter by any of the columns contained in the source catalog.
# with a known RA/Dec
desired_ra = 53.15437048946369
desired_dec = -27.771689847051736
c = SkyCoord(ra=desired_ra*u.degree, dec=desired_dec*u.degree)
nearest_id, distance_2d, distance_3d = c.match_to_catalog_sky(cat['sky_centroid'])
cat[['label', 'xcentroid', 'ycentroid', 'sky_centroid', 'is_extended', 'isophotal_abmag', 'isophotal_vegamag']][nearest_id]
# alternatively with a known X/Y pixel location in the F200W image (based on what you've defined cat to be)
known_x = 1880
known_y = 1100
nearest_pos = [np.sqrt((x-known_x)**2 + (y-known_y)**2) for x, y in zip(cat['xcentroid'], cat['ycentroid'])]
wh_nearest = np.where(np.array(nearest_pos) == min(nearest_pos))[0][0]
cat[['label', 'xcentroid', 'ycentroid', 'sky_centroid', 'is_extended', 'isophotal_abmag', 'isophotal_vegamag']][wh_nearest]
# alternatively with a known source number
# note that this number might be different for you depending on pipeline versions and parameters changed.
# source = 109 # if you want to specify the number
source = cat['label'][nearest_id] # to match the one chosen in this notebook
wh_source = np.where(np.array(cat['label'] == source))[0][0]
cat[['label', 'xcentroid', 'ycentroid', 'sky_centroid', 'is_extended', 'isophotal_abmag', 'isophotal_vegamag']][wh_source]
Using any of the three methods above, write out the new catalog.
with RA/Dec:
new_cat = Table(cat[nearest_id])with x/y:
new_cat = Table(cat[wh_nearest])with source ID:
new_cat = Table(cat[wh_source])
new_cat = Table(cat[wh_source]) # turn the row instance into a dataframe again
# save the new catalog with a unique name
new_cat_name = cat_name.replace('cat.ecsv', f'source{source}_cat.ecsv')
new_cat.write(new_cat_name, overwrite=True)
print('Saved:', new_cat_name)
Once we have an updated source catalog that we are content with, we can move on to the spec2 step of the pipeline. It likely will be necessary to come back to this step after running spec2. Let’s take a quick look at the source that we will be extracting in the following notebook with spec2.
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 6))
img = cust_image3_i2d[-1]
cat = Table.read(new_cat_name)
for ax in [ax1, ax2]:
# plot the image
with fits.open(img) as hdu:
ax.imshow(hdu[1].data, vmin=0, vmax=0.3, origin='lower')
title = f"obs{hdu[0].header['OBSERVTN']} {hdu[0].header['PUPIL']}"
# also plot the associated catalog
extended_sources = cat[cat['is_extended'] == 1] # 1 is True; i.e. is extended
point_sources = cat[cat['is_extended'] == 0] # 0 is False; i.e. is a point source
for color, sources in zip(['darkred', 'black'], [extended_sources, point_sources]):
# plotting the sources
ax.scatter(sources['xcentroid'], sources['ycentroid'], s=20, facecolors='None', edgecolors=color, alpha=0.9)
# adding source labels
for i, source_num in enumerate(sources['label']):
ax.annotate(source_num,
(sources['xcentroid'][i]+0.5, sources['ycentroid'][i]+0.5),
fontsize=8,
color=color)
fig.suptitle("Speicifc Source to Extract with Spec2")
# zooming in on a smaller region
ax2.set_xlim(known_x-50, known_x+50)
ax2.set_ylim(known_y-50, known_y+50)
fig.tight_layout()
