Exploring Star Clusters in the Euclid ERO Data#

Learning Goals#

By the end of this tutorial, you will be able to:

• Access the Euclid ERO images using astroquery

• Create cutouts on the large Euclid ERO images directly

• Extract sources on the Euclid image and run photometry tools

• Compare photometry to the Gaia catalog

• Visualize the Euclid ERO image and the Gaia catalog in Firefly

Introduction#

Euclid is a European Space Agency (ESA) space mission with NASA participation, to study the geometry and nature of the dark Universe. As part of its first observations, Euclid publicly released so-called Early Release Observations (ERO) targeting some press-worthy targets on the sky such as star clusters or local galaxies.

In this notebook, we will focus on the ERO data of NGC 6397, a globular cluster 78,000 light years away. (Note that there is another globular cluster, NGC 6254 that can also be used for this analysis - in fact the user can choose which one to use) The goal of this analysis is to extract the Euclid photometry of the stars belonging to the cluster and compare them to the photometry of Gaia. Due to the different pixel scales of the visible (VIS, \(0.1^{\prime\prime}\)) and near-IR (NISP, \(0.3^{\prime\prime}\)) wavelengths, we will first detect and extract the stars in the VIS filter and use their position to subsequently extract the photometry in the NISP (Y, J, H) filters (a method also referred to as forced photometry).

One challenge with Euclid data is their size due to the large sky coverage and small pixel size. This notebook demonstrates how to download only a cutout of the large Euclid ERO observation image (namely focussing only on the position of the globular cluster. Furthermore, this notebook demonstrates how to extract sources on a large astronomical image and how to measure their photometry across different pixel scales using a very simple implementation of the forced photometry method with positional priors. Finally, we also demonstrate how to visualize the Euclid images and catalog in Firefly, an open-source web-based UI library for astronomical data archive access and visualization developed at Caltech (https://github.com/Caltech-IPAC/firefly). Firefly is a convenient tool to visualize images similar to DS9 on your local computer, but it can run on a cloud-based science platform.

This notebook is written to be used in Fornax which is a cloud based computing platform using AWS. The advantage of this is that the user does not need to download actuall data, hence the analysis of large datasets is not limited by local computing resources. It also allows to access data across all archives fast and easy.

Data Volume#

The total data volume required for running this notebook is less than 20 MB.

Imports#

First, we import all necessary packages.

# Uncomment the next line to install dependencies if needed.
# !pip install tqdm numpy matplotlib astropy photutils>=2.0 astroquery>=0.4.10 sep>=1.4 firefly_client>=3.2

# General imports
import os
import numpy as np
from tqdm import tqdm

# Astroquery
from astroquery.ipac.irsa import Irsa
from astroquery.gaia import Gaia

# Astropy
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.io import fits
from astropy.nddata import Cutout2D
from astropy.wcs import WCS
from astropy.table import Table
from astropy.stats import sigma_clipped_stats

# Photometry tools
import sep
from photutils.detection import DAOStarFinder
from photutils.psf import PSFPhotometry, IterativePSFPhotometry, CircularGaussianSigmaPRF, make_psf_model_image
from photutils.background import LocalBackground, MMMBackground

# Firefly
from firefly_client import FireflyClient

# Matplotlib
import matplotlib.pyplot as plt
import matplotlib as mpl

/home/runner/work/irsa-tutorials/irsa-tutorials/.tox/py311-buildhtml/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Next, we define some parameters for Matplotlib plotting.

## Plotting stuff
mpl.rcParams['font.size'] = 14
mpl.rcParams['axes.labelpad'] = 7
mpl.rcParams['xtick.major.pad'] = 7
mpl.rcParams['ytick.major.pad'] = 7
mpl.rcParams['xtick.minor.visible'] = True
mpl.rcParams['ytick.minor.visible'] = True
mpl.rcParams['xtick.minor.top'] = True
mpl.rcParams['xtick.minor.bottom'] = True
mpl.rcParams['ytick.minor.left'] = True
mpl.rcParams['ytick.minor.right'] = True
mpl.rcParams['xtick.major.size'] = 5
mpl.rcParams['ytick.major.size'] = 5
mpl.rcParams['xtick.minor.size'] = 3
mpl.rcParams['ytick.minor.size'] = 3
mpl.rcParams['xtick.direction'] = 'in'
mpl.rcParams['ytick.direction'] = 'in'
#mpl.rc('text', usetex=True)
mpl.rc('font', family='serif')
mpl.rcParams['xtick.top'] = True
mpl.rcParams['ytick.right'] = True
mpl.rcParams['hatch.linewidth'] = 1

def_cols = plt.rcParams['axes.prop_cycle'].by_key()['color']

Setting up the Environment#

Next, we set up the environment. This includes

setting up an output data directory (will be created if it does not exist)
define the search radius around the target of interest to pull data from the Gaia catalog
define the cutout size that will be used to download a certain part of the Euclid ERO images

Finally, we also define the target of interest here. We can choose between two globular clusters, NGC 6254 and NGC6397, both covered by the Euclid ERO data.

Note that astropy units can be attached to the search_radius and cutout_size.

# create output directory
if os.path.exists("./data/"):
    print("Output directory already created.")
else:
    print("Creating data directory.")
    os.mkdir("./data/")

search_radius = 1.5 * u.arcmin # search radius
cutout_size = 1.5 * u.arcmin # cutout size

coord = SkyCoord.from_name('NGC 6397')

Output directory already created.

Search Euclid ERO Images#

Now, we search for the Euclid ERO images using the astroquery package. Note that the Euclid ERO images are no in the cloud currently, but we access them directly from IRSA using IRSA’s Simple Image Access (SIA) methods.

Note

The following only works for combined images (either extended or point source stacks). This would not work if there are multiple, let’s say, H-band images of Euclid at a given position. Therefore, no time domain studies here (which is anyway not one of the main goals of Euclid).

The IRSA SIA products can be listed via

Irsa.list_collections(servicetype='SIA')

Here we use the collection euclid_ero, containing the Euclid ERO images. We first create a SkyCoord object and then query the SIA.

image_tab = Irsa.query_sia(pos=(coord, search_radius), collection='euclid_ero')
print("Number of images available: {}".format(len(image_tab)))

Number of images available: 25

Sort the queried table by wavelength (column em_min). This allows us later when we plot the images to keep them sorted by wavelength (VIS, Y, J, H).

image_tab.sort('em_min')

Let’s inspec the table that was downloaded.

image_tab[0:3]

Table length=3

s_ra	s_dec	facility_name	instrument_name	dataproduct_subtype	calib_level	dataproduct_type	energy_bandpassname	energy_emband	obs_id	s_resolution	em_min	em_max	em_res_power	proposal_title	access_url	access_format	access_estsize	t_exptime	s_region	obs_collection	obs_intent	algorithm_name	facility_keywords	instrument_keywords	environment_photometric	proposal_id	proposal_pi	proposal_project	target_name	target_type	target_standard	target_moving	target_keywords	obs_release_date	s_xel1	s_xel2	s_pixel_scale	position_timedependent	t_min	t_max	t_resolution	t_xel	obs_publisher_did	s_fov	em_xel	pol_states	pol_xel	cloud_access	o_ucd	upload_row_id
deg	deg									arcsec	m	m					kbyte	s	deg																		arcsec		d	d	s			deg
float64	float64	object	object	object	int16	object	object	object	object	float64	float64	float64	float64	object	object	object	int64	float64	object	object	object	object	object	object	bool	object	object	object	object	object	bool	bool	object	object	int64	int64	float64	bool	float64	float64	float64	int64	object	float64	int64	object	int64	object	object	int64
265.1744234	-53.6588139	Euclid	VIS	weight	3	image	VIS	Optical	VIS-ERO-NGC6397	0.16	5.5e-07	9e-07	2.1	Euclid Early Release Observations	https://irsa.ipac.caltech.edu/data/Euclid/ERO/images/NGC6397/ERO-NGC6397/Euclid-VIS-ERO-NGC6397-Flattened.DR3.weight.fits	image/fits	5184009	--	POLYGON ICRS 266.0282156 -54.1557675 266.0082 -53.1559013 264.3406 -53.1559013 264.3205844 -54.1557675 266.0282156 -54.1557675	euclid_ero	science	coadd			--			Euclid	NGC 6397	object	False	False		2024-05-23 00:00:00	36000	36000	0.1000000000008	False	--	--	--	--	ivo://irsa.ipac/euclid_ero?VIS-ERO-NGC6397/Flattened	1.000000000008	--		--	{}		1
265.1744234	-53.6588139	Euclid	VIS	science	3	image	VIS	Optical	VIS-ERO-NGC6397	0.16	5.5e-07	9e-07	2.1	Euclid Early Release Observations	https://irsa.ipac.caltech.edu/data/Euclid/ERO/images/NGC6397/ERO-NGC6397/Euclid-VIS-ERO-NGC6397-Flattened.DR3.fits	image/fits	5184009	--	POLYGON ICRS 266.0282156 -54.1557675 266.0082 -53.1559013 264.3406 -53.1559013 264.3205844 -54.1557675 266.0282156 -54.1557675	euclid_ero	science	coadd			--			Euclid	NGC 6397	object	False	False		2024-05-23 00:00:00	36000	36000	0.1000000000008	False	--	--	--	--	ivo://irsa.ipac/euclid_ero?VIS-ERO-NGC6397/Flattened	1.000000000008	--		--	{}		1
265.1744234	-53.6588139	Euclid	VIS	auxiliary	3	image	VIS	Optical	VIS-ERO-NGC6397	0.16	5.5e-07	9e-07	2.1	Euclid Early Release Observations	https://irsa.ipac.caltech.edu/data/Euclid/ERO/images/NGC6397/ERO-NGC6397/Euclid-VIS-Mask-ERO-NGC6397-LSB.DR3.fits	image/fits	5184006	--	POLYGON ICRS 266.0282156 -54.1557675 266.0082 -53.1559013 264.3406 -53.1559013 264.3205844 -54.1557675 266.0282156 -54.1557675	euclid_ero	science	coadd			--			Euclid	NGC 6397	object	False	False		2024-05-23 00:00:00	36000	36000	0.1000000000008	False	--	--	--	--	ivo://irsa.ipac/euclid_ero?VIS-ERO-NGC6397/LSB	1.000000000008	--		--	{}		1

Next, we add additional restrictions to narrow down the images.

The Euclid ERO images come in two different flavors:

Flattened images: idealized for compact sources (for example stars)
LSB images: idealized for low surface brightness objects (for example galaxies)

Maybe counter-intuitively, we select the LSB images here by checking if the file name given in the URL (access_url column) contains that key word. We found that the Flattened images contain many masked pixels.

#sel_basic = np.where( ["Flattened" in tt["access_url"] for tt in image_tab] )[0]
sel_basic = np.where( ["LSB" in tt["access_url"] for tt in image_tab] )[0]
image_tab = image_tab[sel_basic]

Next, we want to check what filteres are available. This can be done by np.unique(), however, in that case it would sort the filters alphabetically. We want to keep the sorting based on the wavelength, therefore we opt for a more complicated way.

idxs = np.unique(image_tab["energy_bandpassname"], return_index=True)[1]
filters = [image_tab["energy_bandpassname"][idx] for idx in sorted(idxs)]
print("Filters: {}".format(filters))

Filters: ['VIS', 'Y', 'J', 'H']

We can now loop throught the filters and collect the images to create a handy summary table with all the data we have. This way we will have an easier time to later access the data.

# Create a dictionary with all the necessary information (science, weight, noise, mask)
summary_table = Table(names=["filter","products","facility_name","instrument_name"] , dtype=[str,str,str,str])
for filt in filters:
    sel = np.where(image_tab['energy_bandpassname'] == filt)[0]
    products = list( np.unique(image_tab["dataproduct_subtype"][sel].value) )
    if "science" in products: # sort so that science is the first in the list. This is the order we create the hdu extensions
        products.remove("science")
        products.insert(0,"science")
    else:
        print("WARNING: no science image found!")
    print(products)

    summary_table.add_row( [filt , ";".join(products), str(np.unique(image_tab["facility_name"][sel].value)[0]), str(np.unique(image_tab["instrument_name"][sel].value)[0])] )

['science', 'auxiliary']
['science', 'auxiliary']
['science', 'auxiliary']
['science', 'auxiliary']

Let’s check out the summary table that we have created. We see that we have all the 4 Euclid bands and what data products are available for each of them.

summary_table

Table length=4

filter	products	facility_name	instrument_name
str3	str17	str6	str4
VIS	science;auxiliary	Euclid	VIS
Y	science;auxiliary	Euclid	NISP
J	science;auxiliary	Euclid	NISP
H	science;auxiliary	Euclid	NISP

Create Cutout Images#

Now that we have a list of data products, we can create the cutouts. This is important as the full Euclid ERO images would be too large to run extraction and photometry software on them (they would simply fail due to memory issues).

For each image, we create a cutout around the target of interest, using the cutout_size defined above. The cutouts are then collected into HDUs. That way we can easily access the different data products. Note that we only use the science products as the ancillary products are not needed.

We save the HDU to disk as it will be later used when we visualize the Euclid ERO FITS images in Firefly.

Note

You will notice that Cutout2D can be applied to an URL. That way, it we do not need to download the full image to create a cutout. This is a useful trick to keep in mind when analyzing large images. This makes creating cutout images very fast.

%%time
for ii,filt in tqdm(enumerate(filters)):
    products = summary_table[summary_table["filter"] == filt]["products"][0].split(";")

    for product in products:
        sel = np.where( (image_tab["energy_bandpassname"] == filt) & (image_tab["dataproduct_subtype"] == product) )[0]

        with fits.open(image_tab['access_url'][sel[0]], use_fsspec=True) as hdul:
            tmp = Cutout2D(hdul[0].section, position=coord, size=cutout_size, wcs=WCS(hdul[0].header)) # create cutout


            if (product == "science") & (ii == 0): # if science image, then create a new hdu.
                hdu0 = fits.PrimaryHDU(data = tmp.data, header=hdul[0].header)
                hdu0.header.update(tmp.wcs.to_header()) # update header with WCS info
                hdu0.header["EXTNAME"] = "{}_{}".format(filt,product.upper())
                hdu0.header["PRODUCT"] = product.upper()
                hdu0.header["FILTER"] = filt.upper()
                hdulcutout = fits.HDUList([hdu0])
            elif (product == "science") & (ii > 0):
                hdu = fits.ImageHDU(data = tmp.data, header=hdul[0].header)
                hdu.header.update(tmp.wcs.to_header()) # update header with WCS info
                hdu.header["EXTNAME"] = "{}_{}".format(filt,product.upper())
                hdu.header["PRODUCT"] = product.upper()
                hdu.header["FILTER"] = filt.upper()
                hdulcutout.append(hdu)

## Save the HDUL cube:
hdulcutout.writeto("./data/euclid_images_test.fits", overwrite=True)

CPU times: user 1.57 s, sys: 326 ms, total: 1.9 s
Wall time: 26.1 s

0it [00:00, ?it/s]
1it [00:18, 18.24s/it]
2it [00:20,  8.83s/it]
3it [00:23,  6.29s/it]
4it [00:26,  4.71s/it]
4it [00:26,  6.51s/it]

Let’s look at the HDU that we created to check what information we have. You see that all filters are collected in different extensions. Also note the different dimensions of the FITS layers as the pixel scale of VIS and NISP are different.

hdulcutout.info()

Filename: (No file associated with this HDUList)
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  VIS_SCIENCE    1 PrimaryHDU     190   (900, 900)   float32   
  1  Y_SCIENCE     1 ImageHDU       191   (300, 300)   float32   
  2  J_SCIENCE     1 ImageHDU       191   (300, 300)   float32   
  3  H_SCIENCE     1 ImageHDU       191   (300, 300)   float32

We can now plot the image cutouts that we generated. The globular cluster is clearly visible.

nimages = len(filters) # number of images
ncols = int(4) # number of columns
nrows = int( nimages // ncols ) # number of rows

fig = plt.figure(figsize = (5*ncols,5*nrows) )
axs = [fig.add_subplot(nrows,ncols,ii+1) for ii in range(nimages)]

for ii,filt in enumerate(filters):
    img = hdulcutout["{}_SCIENCE".format(filt)].data
    axs[ii].imshow(img , origin="lower")
    axs[ii].text(0.05,0.05 , "{} ({}/{})".format(summary_table["facility_name"][ii],summary_table["instrument_name"][ii],filt) , fontsize=14 , color="white",
                 va="bottom", ha="left" , transform=axs[ii].transAxes)

plt.show()

../../_images/5124054e8573b0c569e89f846638cb597919dc69e838f9f83c1faaf7ab4efe4e.png

Extract Sources and Measure their Photometry on the VIS Image#

Now that we have the images in memory (and on disk - but we do not need them, yet), we can measure the fluxes of the individual stars. Our simple photometry pipeline has different parts:

First, we are using the Python package sep (similar to SExtractor) to extract the position of the sources. We do that on the VIS image, which provides the highest resolution.
Second, we use sep to perform aperture measurements of the photometry. We will use that to compare the obtained fluxes to our forced photometry method
Third, we apply a PSF fitting technique (using the photutils Python package) to improve the photometry measurement

We start by extracting the sources using sep. We first isolate the data that we want to look at (the VIS image only).

## Get Data (this will be replaced later)
img = hdulcutout["VIS_SCIENCE"].data
hdr = hdulcutout["VIS_SCIENCE"].header
img[img == 0] = np.nan

There are some NaN value on the image that we need to mask out. For this we create a mask image that we later feed to sep.

mask = np.isnan(img)

Next, we compute the background statistics what will be used by sep to extract the sources above a certain threshold.

mean, median, std = sigma_clipped_stats(img, sigma=3.0)
print(np.array((mean, median, std)))

[135.95872  99.62191 100.63356]

WARNING: Input data contains invalid values (NaNs or infs), which were automatically clipped. [astropy.stats.sigma_clipping]

Finally, we perform object detection with sep. There are also modules in photutils to do that, however, we found that sep works best here. We output the number of objects found on the image.

objects = sep.extract(img-median, thresh=1.2, err=std, minarea=3, mask=mask, deblend_cont=0.0002, deblend_nthresh=64 )
print("Number of sources extracted: ", len(objects))

Number of sources extracted:  2759

Next, we perform simple aperture photometry on the detected sources. Again, we use the sep package for this. We will use these aperture photometry later to compare to the PSF photometry.

flux, fluxerr, flag = sep.sum_circle(img-median, objects['x'], objects['y'],r=3.0, err=std, gain=1.0)

Now, we use the photutils Python package to perform PSF fitting. Here we assume a simple Gaussian with a FWHM given by psf_fwhm as PSF.

Note

We use a Gaussian PSF here for simplicity. The photometry can be improved by using a pixelated PSF measured directly from the Euclid images (for example by stacking stars).

psf_fwhm = 0.16 # PSF FWHM in arcsec
pixscale = 0.1 # VIS pixel scale in arcsec/px

init_params = Table([objects["x"],objects["y"]] , names=["x","y"]) # initial positions
psf_model = CircularGaussianSigmaPRF(flux=1, sigma=psf_fwhm/pixscale / 2.35)
psf_model.x_0.fixed = True
psf_model.y_0.fixed = True
psf_model.sigma.fixed = False
fit_shape = (5, 5)
psfphot = PSFPhotometry(psf_model,
                        fit_shape,
                        finder = DAOStarFinder(fwhm=0.1, threshold=3.*std, exclude_border=True), # not needed because we are using fixed initial positions.
                        aperture_radius = 4,
                        progress_bar = True)

After initiating the PSFPhotometry object, we can run the PSF photometry measurement. This can take a while (typically between 1 and 2 minutes).

phot = psfphot(img-median, error=None, mask=mask, init_params=init_params)

Fit source/group:   0%|          | 0/2759 [00:00<?, ?it/s]
Fit source/group:   1%|          | 25/2759 [00:00<00:11, 248.53it/s]
Fit source/group:   2%|▏         | 50/2759 [00:00<00:11, 237.64it/s]
Fit source/group:   3%|▎         | 76/2759 [00:00<00:10, 246.56it/s]
Fit source/group:   4%|▎         | 102/2759 [00:00<00:10, 248.15it/s]
Fit source/group:   5%|▍         | 129/2759 [00:00<00:10, 253.31it/s]
Fit source/group:   6%|▌         | 155/2759 [00:00<00:10, 246.56it/s]
Fit source/group:   7%|▋         | 180/2759 [00:00<00:10, 243.56it/s]
Fit source/group:   7%|▋         | 205/2759 [00:00<00:10, 244.52it/s]
Fit source/group:   8%|▊         | 232/2759 [00:00<00:10, 249.69it/s]
Fit source/group:   9%|▉         | 258/2759 [00:01<00:09, 251.66it/s]
Fit source/group:  10%|█         | 284/2759 [00:01<00:09, 248.59it/s]
Fit source/group:  11%|█         | 310/2759 [00:01<00:09, 251.85it/s]
Fit source/group:  12%|█▏        | 336/2759 [00:01<00:09, 252.89it/s]
Fit source/group:  13%|█▎        | 362/2759 [00:01<00:09, 254.23it/s]
Fit source/group:  14%|█▍        | 388/2759 [00:01<00:10, 232.48it/s]
Fit source/group:  15%|█▍        | 413/2759 [00:01<00:09, 236.04it/s]
Fit source/group:  16%|█▌        | 437/2759 [00:01<00:11, 194.89it/s]
Fit source/group:  17%|█▋        | 462/2759 [00:01<00:11, 207.49it/s]
Fit source/group:  18%|█▊        | 488/2759 [00:02<00:10, 220.12it/s]
Fit source/group:  19%|█▊        | 513/2759 [00:02<00:09, 226.87it/s]
Fit source/group:  19%|█▉        | 537/2759 [00:02<00:09, 226.48it/s]
Fit source/group:  20%|██        | 562/2759 [00:02<00:09, 230.93it/s]
Fit source/group:  21%|██▏       | 587/2759 [00:02<00:09, 235.32it/s]
Fit source/group:  22%|██▏       | 611/2759 [00:02<00:09, 228.01it/s]
Fit source/group:  23%|██▎       | 635/2759 [00:02<00:09, 225.73it/s]
Fit source/group:  24%|██▍       | 661/2759 [00:02<00:09, 232.72it/s]
Fit source/group:  25%|██▍       | 685/2759 [00:02<00:09, 218.79it/s]
Fit source/group:  26%|██▌       | 710/2759 [00:03<00:09, 225.52it/s]
Fit source/group:  27%|██▋       | 733/2759 [00:03<00:08, 226.32it/s]
Fit source/group:  27%|██▋       | 757/2759 [00:03<00:08, 229.92it/s]
Fit source/group:  28%|██▊       | 782/2759 [00:03<00:08, 233.50it/s]
Fit source/group:  29%|██▉       | 807/2759 [00:03<00:08, 238.00it/s]
Fit source/group:  30%|███       | 831/2759 [00:03<00:08, 236.29it/s]
Fit source/group:  31%|███       | 856/2759 [00:03<00:07, 239.83it/s]
Fit source/group:  32%|███▏      | 882/2759 [00:03<00:07, 243.17it/s]
Fit source/group:  33%|███▎      | 907/2759 [00:03<00:07, 243.29it/s]
Fit source/group:  34%|███▍      | 933/2759 [00:03<00:07, 246.36it/s]
Fit source/group:  35%|███▍      | 958/2759 [00:04<00:07, 227.98it/s]
Fit source/group:  36%|███▌      | 982/2759 [00:04<00:07, 230.33it/s]
Fit source/group:  36%|███▋      | 1006/2759 [00:04<00:07, 226.35it/s]
Fit source/group:  37%|███▋      | 1029/2759 [00:04<00:07, 225.98it/s]
Fit source/group:  38%|███▊      | 1055/2759 [00:04<00:07, 233.30it/s]
Fit source/group:  39%|███▉      | 1079/2759 [00:04<00:07, 234.32it/s]
Fit source/group:  40%|████      | 1104/2759 [00:04<00:06, 236.82it/s]
Fit source/group:  41%|████      | 1130/2759 [00:04<00:06, 242.37it/s]
Fit source/group:  42%|████▏     | 1156/2759 [00:04<00:06, 247.20it/s]
Fit source/group:  43%|████▎     | 1182/2759 [00:05<00:06, 248.36it/s]
Fit source/group:  44%|████▍     | 1208/2759 [00:05<00:06, 249.01it/s]
Fit source/group:  45%|████▍     | 1233/2759 [00:05<00:06, 244.67it/s]
Fit source/group:  46%|████▌     | 1260/2759 [00:05<00:06, 249.60it/s]
Fit source/group:  47%|████▋     | 1285/2759 [00:05<00:06, 243.50it/s]
Fit source/group:  47%|████▋     | 1310/2759 [00:05<00:06, 234.60it/s]
Fit source/group:  48%|████▊     | 1334/2759 [00:05<00:06, 234.31it/s]
Fit source/group:  49%|████▉     | 1359/2759 [00:05<00:05, 235.62it/s]
Fit source/group:  50%|█████     | 1383/2759 [00:05<00:06, 228.32it/s]
Fit source/group:  51%|█████     | 1407/2759 [00:05<00:06, 216.46it/s]
Fit source/group:  52%|█████▏    | 1432/2759 [00:06<00:05, 224.38it/s]
Fit source/group:  53%|█████▎    | 1458/2759 [00:06<00:05, 232.98it/s]
Fit source/group:  54%|█████▍    | 1483/2759 [00:06<00:05, 237.78it/s]
Fit source/group:  55%|█████▍    | 1507/2759 [00:06<00:05, 238.01it/s]
Fit source/group:  55%|█████▌    | 1531/2759 [00:06<00:05, 219.35it/s]
Fit source/group:  56%|█████▋    | 1555/2759 [00:06<00:05, 224.85it/s]
Fit source/group:  57%|█████▋    | 1582/2759 [00:06<00:04, 236.18it/s]
Fit source/group:  58%|█████▊    | 1607/2759 [00:06<00:04, 238.65it/s]
Fit source/group:  59%|█████▉    | 1632/2759 [00:06<00:04, 235.24it/s]
Fit source/group:  60%|██████    | 1656/2759 [00:07<00:04, 227.23it/s]
Fit source/group:  61%|██████    | 1681/2759 [00:07<00:04, 232.86it/s]
Fit source/group:  62%|██████▏   | 1707/2759 [00:07<00:04, 239.27it/s]
Fit source/group:  63%|██████▎   | 1733/2759 [00:07<00:04, 244.93it/s]
Fit source/group:  64%|██████▎   | 1758/2759 [00:07<00:04, 233.25it/s]
Fit source/group:  65%|██████▍   | 1782/2759 [00:07<00:04, 226.15it/s]
Fit source/group:  66%|██████▌   | 1808/2759 [00:07<00:04, 233.84it/s]
Fit source/group:  66%|██████▋   | 1832/2759 [00:07<00:03, 233.05it/s]
Fit source/group:  67%|██████▋   | 1860/2759 [00:07<00:03, 246.35it/s]
Fit source/group:  68%|██████▊   | 1885/2759 [00:08<00:03, 245.22it/s]
Fit source/group:  69%|██████▉   | 1911/2759 [00:08<00:03, 246.39it/s]
Fit source/group:  70%|███████   | 1937/2759 [00:08<00:03, 248.26it/s]
Fit source/group:  71%|███████   | 1962/2759 [00:08<00:03, 244.61it/s]
Fit source/group:  72%|███████▏  | 1987/2759 [00:08<00:03, 225.85it/s]
Fit source/group:  73%|███████▎  | 2012/2759 [00:08<00:03, 231.50it/s]
Fit source/group:  74%|███████▍  | 2039/2759 [00:08<00:02, 241.26it/s]
Fit source/group:  75%|███████▍  | 2064/2759 [00:08<00:02, 235.59it/s]
Fit source/group:  76%|███████▌  | 2088/2759 [00:08<00:03, 216.17it/s]
Fit source/group:  76%|███████▋  | 2110/2759 [00:09<00:03, 211.76it/s]
Fit source/group:  77%|███████▋  | 2134/2759 [00:09<00:02, 219.06it/s]
Fit source/group:  78%|███████▊  | 2161/2759 [00:09<00:02, 231.44it/s]
Fit source/group:  79%|███████▉  | 2188/2759 [00:09<00:02, 240.77it/s]
Fit source/group:  80%|████████  | 2214/2759 [00:09<00:02, 244.46it/s]
Fit source/group:  81%|████████  | 2241/2759 [00:09<00:02, 249.34it/s]
Fit source/group:  82%|████████▏ | 2269/2759 [00:09<00:01, 257.87it/s]
Fit source/group:  83%|████████▎ | 2295/2759 [00:09<00:01, 255.78it/s]
Fit source/group:  84%|████████▍ | 2322/2759 [00:09<00:01, 259.89it/s]
Fit source/group:  85%|████████▌ | 2349/2759 [00:09<00:01, 257.20it/s]
Fit source/group:  86%|████████▌ | 2377/2759 [00:10<00:01, 262.45it/s]
Fit source/group:  87%|████████▋ | 2404/2759 [00:10<00:01, 260.46it/s]
Fit source/group:  88%|████████▊ | 2431/2759 [00:10<00:01, 235.14it/s]
Fit source/group:  89%|████████▉ | 2456/2759 [00:10<00:01, 230.68it/s]
Fit source/group:  90%|████████▉ | 2482/2759 [00:10<00:01, 238.18it/s]
Fit source/group:  91%|█████████ | 2507/2759 [00:10<00:01, 193.96it/s]
Fit source/group:  92%|█████████▏| 2528/2759 [00:10<00:01, 169.25it/s]
Fit source/group:  93%|█████████▎| 2553/2759 [00:10<00:01, 185.89it/s]
Fit source/group:  93%|█████████▎| 2574/2759 [00:11<00:01, 169.30it/s]
Fit source/group:  94%|█████████▍| 2593/2759 [00:11<00:01, 147.19it/s]
Fit source/group:  95%|█████████▍| 2609/2759 [00:11<00:01, 137.79it/s]
Fit source/group:  95%|█████████▌| 2624/2759 [00:11<00:00, 140.44it/s]
Fit source/group:  96%|█████████▌| 2639/2759 [00:11<00:00, 132.28it/s]
Fit source/group:  96%|█████████▌| 2654/2759 [00:11<00:00, 125.75it/s]
Fit source/group:  97%|█████████▋| 2669/2759 [00:11<00:00, 131.26it/s]
Fit source/group:  97%|█████████▋| 2685/2759 [00:12<00:00, 134.45it/s]
Fit source/group:  98%|█████████▊| 2699/2759 [00:12<00:00, 102.39it/s]
Fit source/group:  98%|█████████▊| 2716/2759 [00:12<00:00, 116.35it/s]
Fit source/group:  99%|█████████▉| 2729/2759 [00:12<00:00, 117.03it/s]
Fit source/group: 100%|█████████▉| 2754/2759 [00:12<00:00, 150.06it/s]
Fit source/group: 100%|██████████| 2759/2759 [00:12<00:00, 219.54it/s]
WARNING: One or more fit(s) may not have converged. Please check the "flags" column in the output table. [photutils.psf.photometry]

Once this is done, we can create a residual image.

resimage = psfphot.make_residual_image(data = img-median, psf_shape = (9, 9))

Add model sources:   0%|          | 0/2759 [00:00<?, ?it/s]
Add model sources:  19%|█▊        | 516/2759 [00:00<00:00, 5151.16it/s]
Add model sources:  38%|███▊      | 1042/2759 [00:00<00:00, 5210.13it/s]
Add model sources:  58%|█████▊    | 1596/2759 [00:00<00:00, 5356.99it/s]
Add model sources:  78%|███████▊  | 2143/2759 [00:00<00:00, 5400.76it/s]
Add model sources:  97%|█████████▋| 2684/2759 [00:00<00:00, 5137.40it/s]
Add model sources: 100%|██████████| 2759/2759 [00:00<00:00, 5190.41it/s]

We now want to add the best-fit coordinates (R.A. and Decl.) to the VIS photometry catalog. For this, we have to convert the image coordinates into sky coordinates using the WCS information. We will need these coordinates because we want to use them as positional priors for the photometry measurement on the NISP images.

## Add coordinates to catalog
wcs1 = WCS(hdr) # VIS
radec = wcs1.all_pix2world(phot["x_fit"],phot["y_fit"],0)
phot["ra_fit"] = radec[0]
phot["dec_fit"] = radec[1]

Finally, we plot the VIS image and the residual with the extracted sources overlaid.

fig = plt.figure(figsize=(20,10))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.imshow(np.log10(img), cmap="Greys", origin="lower")
ax1.plot(phot["x_fit"], phot["y_fit"] , "o", markersize=8 , markeredgecolor="red", fillstyle="none")

ax2.imshow(resimage,vmin=-5*std, vmax=5*std, cmap="RdBu", origin="lower")
ax2.plot(phot["x_fit"], phot["y_fit"] , "o", markersize=8 , markeredgecolor="red", fillstyle="none")

plt.show()

../../_images/5310e7d4e4d3ddae2b8fd76a93c625e9a91fd785f6b5d4a6913c3aa0c8a499d8.png

As an additional check, we can compare the aperture photometry and the PSF photometry. We should find a good agreement between those two measurement methods. However, note that the PSF photometry should do a better job in deblending the fluxes of sources that are close by.

x = objects["flux"]
y = phot["flux_fit"]

fig = plt.figure(figsize=(6,5))
ax1 = fig.add_subplot(1,1,1)

ax1.plot(x , y , "o", markersize=2, alpha=0.5)
minlim = np.nanmin(np.concatenate((x,y)))
maxlim = np.nanmax(np.concatenate((x,y)))

ax1.fill_between(np.asarray([minlim,maxlim]),np.asarray([minlim,maxlim])/1.5,np.asarray([minlim,maxlim])*1.5, color="gray", alpha=0.2, linewidth=0)
ax1.fill_between(np.asarray([minlim,maxlim]),np.asarray([minlim,maxlim])/1.2,np.asarray([minlim,maxlim])*1.2, color="gray", alpha=0.4, linewidth=0)
ax1.plot(np.asarray([minlim,maxlim]),np.asarray([minlim,maxlim]), ":", color="gray")

ax1.set_xlabel("Aperture Photometry [flux]")
ax1.set_ylabel("PSF forced-photometry [flux]")
ax1.set_xscale('log')
ax1.set_yscale('log')
plt.show()

../../_images/17aa75c20aa9ce7d5bcf9073811354441aa9c0803b1b8dfa46eef8ca6bfdfbfe.png

Measure the Photometry on the NISP Images#

We now have the photometry and the position of sources on the VIS image. We can now proceed with similar steps on the NISP images. Because the NISP PSF and pixel scale are larger that those of the VIS images, we utilize the advantage of position prior-based forced photometry. For this, we use the positions of the VIS measurements and perform PSF fitting on the NISP image using these priors.

The steps below are almost identical to the ones applied to the VIS images.

We first isolate the data, which is in this case the NISP H-band filter. Note that this is an arbitrary choice and you should be encouraged to try other filters as well!

img2 = hdulcutout["H_SCIENCE"].data
hdr2 = hdulcutout["H_SCIENCE"].header
img2[img2 == 0] = np.nan

We again need to create a mask that will be fed to the PSF fitting module.

mask2 = np.isnan(img2)

… and we also get some statistics on the sky background.

mean2, median2, std2 = sigma_clipped_stats(img2, sigma=3.0)
print(np.array((mean2, median2, std2)))

[2626.302  1889.8518 2043.9843]

Now, we need to obtain the (x,y) image coordinates on the NISP image that correspond to the extracted sources on the VIS image. We use the WCS information from the NISP image for this case and apply it to the sky coordinates obtained on the VIS image.

wcs = WCS(hdr) # VIS
wcs2 = WCS(hdr2) # NISP
radec = wcs.all_pix2world(objects["x"],objects["y"], 0)
xy = wcs2.all_world2pix(radec[0],radec[1],0)

Having all this set up, we can now again perform the PSF photometry using PSFPhotometry() from the photutils package. This again can take a while, typically 1 minute.

psf_fwhm = 0.3 # arcsec
pixscale = 0.3 # arcsec/px

init_params = Table([xy[0],xy[1]] , names=["x","y"]) # initial positions
psf_model = CircularGaussianSigmaPRF(flux=1, sigma=psf_fwhm/pixscale / 2.35)
psf_model.x_0.fixed = True
psf_model.y_0.fixed = True
psf_model.sigma.fixed = False
fit_shape = (3, 3)
psfphot2 = PSFPhotometry(psf_model,
                        fit_shape,
                        finder = DAOStarFinder(fwhm=0.1, threshold=3.*std2, exclude_border=True), # not needed because we are using fixed initial positions.
                        aperture_radius = 4,
                        progress_bar = True)
phot2 = psfphot2(img2-median2, error=None, mask=mask2, init_params=init_params)
resimage2 = psfphot2.make_residual_image(data = img2-median2, psf_shape = (3, 3))

Fit source/group:   0%|          | 0/2759 [00:00<?, ?it/s]
Fit source/group:   0%|          | 9/2759 [00:00<00:31, 87.80it/s]
Fit source/group:   1%|          | 24/2759 [00:00<00:22, 123.43it/s]
Fit source/group:   1%|▏         | 41/2759 [00:00<00:19, 140.45it/s]
Fit source/group:   2%|▏         | 56/2759 [00:00<00:20, 129.45it/s]
Fit source/group:   3%|▎         | 75/2759 [00:00<00:18, 141.98it/s]
Fit source/group:   3%|▎         | 90/2759 [00:00<00:22, 120.10it/s]
Fit source/group:   4%|▍         | 104/2759 [00:00<00:21, 125.24it/s]
Fit source/group:   4%|▍         | 117/2759 [00:00<00:22, 117.87it/s]
Fit source/group:   5%|▍         | 133/2759 [00:01<00:20, 128.14it/s]
Fit source/group:   5%|▌         | 147/2759 [00:01<00:23, 110.75it/s]
Fit source/group:   6%|▌         | 159/2759 [00:01<00:23, 108.65it/s]
Fit source/group:   6%|▌         | 171/2759 [00:01<00:25, 100.84it/s]
Fit source/group:   7%|▋         | 192/2759 [00:01<00:21, 119.27it/s]
Fit source/group:   8%|▊         | 208/2759 [00:01<00:19, 128.77it/s]
Fit source/group:   8%|▊         | 229/2759 [00:01<00:17, 143.19it/s]
Fit source/group:   9%|▉         | 245/2759 [00:01<00:17, 145.44it/s]
Fit source/group:   9%|▉         | 260/2759 [00:02<00:17, 146.06it/s]
Fit source/group:  10%|▉         | 275/2759 [00:02<00:18, 137.65it/s]
Fit source/group:  11%|█         | 291/2759 [00:02<00:17, 143.00it/s]
Fit source/group:  11%|█         | 306/2759 [00:02<00:18, 131.63it/s]
Fit source/group:  12%|█▏        | 320/2759 [00:02<00:23, 101.67it/s]
Fit source/group:  12%|█▏        | 342/2759 [00:02<00:18, 127.91it/s]
Fit source/group:  13%|█▎        | 361/2759 [00:02<00:17, 139.18it/s]
Fit source/group:  14%|█▎        | 377/2759 [00:02<00:18, 131.51it/s]
Fit source/group:  14%|█▍        | 392/2759 [00:03<00:17, 135.42it/s]
Fit source/group:  15%|█▍        | 408/2759 [00:03<00:16, 140.24it/s]
Fit source/group:  16%|█▌        | 429/2759 [00:03<00:14, 158.84it/s]
Fit source/group:  16%|█▌        | 446/2759 [00:03<00:14, 158.70it/s]
Fit source/group:  17%|█▋        | 463/2759 [00:03<00:15, 146.25it/s]
Fit source/group:  17%|█▋        | 479/2759 [00:03<00:16, 135.42it/s]
Fit source/group:  18%|█▊        | 500/2759 [00:03<00:14, 153.87it/s]
Fit source/group:  19%|█▊        | 517/2759 [00:03<00:15, 143.18it/s]
Fit source/group:  19%|█▉        | 532/2759 [00:04<00:16, 133.79it/s]
Fit source/group:  20%|█▉        | 546/2759 [00:04<00:17, 123.84it/s]
Fit source/group:  20%|██        | 560/2759 [00:04<00:18, 120.10it/s]
Fit source/group:  21%|██        | 576/2759 [00:04<00:18, 120.38it/s]
Fit source/group:  21%|██▏       | 591/2759 [00:04<00:17, 127.32it/s]
Fit source/group:  22%|██▏       | 604/2759 [00:04<00:18, 117.94it/s]
Fit source/group:  22%|██▏       | 617/2759 [00:04<00:21, 97.71it/s]
Fit source/group:  23%|██▎       | 634/2759 [00:04<00:19, 110.58it/s]
Fit source/group:  24%|██▍       | 658/2759 [00:05<00:15, 131.41it/s]
Fit source/group:  24%|██▍       | 674/2759 [00:05<00:15, 138.25it/s]
Fit source/group:  25%|██▍       | 689/2759 [00:05<00:16, 128.66it/s]
Fit source/group:  25%|██▌       | 703/2759 [00:05<00:19, 107.15it/s]
Fit source/group:  26%|██▌       | 715/2759 [00:05<00:19, 105.50it/s]
Fit source/group:  26%|██▋       | 727/2759 [00:05<00:19, 104.65it/s]
Fit source/group:  27%|██▋       | 746/2759 [00:05<00:16, 122.55it/s]
Fit source/group:  28%|██▊       | 759/2759 [00:06<00:18, 108.55it/s]
Fit source/group:  28%|██▊       | 778/2759 [00:06<00:15, 123.89it/s]
Fit source/group:  29%|██▉       | 796/2759 [00:06<00:15, 130.69it/s]
Fit source/group:  30%|██▉       | 818/2759 [00:06<00:13, 141.25it/s]
Fit source/group:  30%|███       | 833/2759 [00:06<00:14, 131.39it/s]
Fit source/group:  31%|███       | 849/2759 [00:06<00:14, 127.67it/s]
Fit source/group:  31%|███▏      | 865/2759 [00:06<00:14, 135.15it/s]
Fit source/group:  32%|███▏      | 879/2759 [00:06<00:14, 129.57it/s]
Fit source/group:  32%|███▏      | 895/2759 [00:07<00:13, 135.88it/s]
Fit source/group:  33%|███▎      | 910/2759 [00:07<00:14, 129.82it/s]
Fit source/group:  34%|███▎      | 926/2759 [00:07<00:13, 137.60it/s]
Fit source/group:  34%|███▍      | 945/2759 [00:07<00:12, 145.79it/s]
Fit source/group:  35%|███▍      | 960/2759 [00:07<00:13, 133.90it/s]
Fit source/group:  35%|███▌      | 974/2759 [00:07<00:14, 122.45it/s]
Fit source/group:  36%|███▌      | 987/2759 [00:07<00:16, 108.21it/s]
Fit source/group:  36%|███▌      | 999/2759 [00:07<00:17, 101.51it/s]
Fit source/group:  37%|███▋      | 1021/2759 [00:08<00:13, 129.62it/s]
Fit source/group:  38%|███▊      | 1038/2759 [00:08<00:12, 138.26it/s]
Fit source/group:  38%|███▊      | 1053/2759 [00:08<00:14, 118.65it/s]
Fit source/group:  39%|███▉      | 1072/2759 [00:08<00:13, 128.68it/s]
Fit source/group:  39%|███▉      | 1088/2759 [00:08<00:13, 125.52it/s]
Fit source/group:  40%|████      | 1109/2759 [00:08<00:12, 135.73it/s]
Fit source/group:  41%|████      | 1123/2759 [00:08<00:12, 126.17it/s]
Fit source/group:  41%|████      | 1136/2759 [00:08<00:13, 120.55it/s]
Fit source/group:  42%|████▏     | 1149/2759 [00:09<00:14, 107.92it/s]
Fit source/group:  42%|████▏     | 1165/2759 [00:09<00:13, 119.38it/s]
Fit source/group:  43%|████▎     | 1183/2759 [00:09<00:12, 130.86it/s]
Fit source/group:  43%|████▎     | 1197/2759 [00:09<00:12, 128.82it/s]
Fit source/group:  44%|████▍     | 1211/2759 [00:09<00:13, 111.27it/s]
Fit source/group:  44%|████▍     | 1223/2759 [00:09<00:18, 84.74it/s]
Fit source/group:  45%|████▍     | 1233/2759 [00:09<00:17, 86.41it/s]
Fit source/group:  45%|████▌     | 1246/2759 [00:10<00:16, 89.79it/s]
Fit source/group:  46%|████▌     | 1267/2759 [00:10<00:12, 117.31it/s]
Fit source/group:  47%|████▋     | 1285/2759 [00:10<00:11, 130.16it/s]
Fit source/group:  47%|████▋     | 1300/2759 [00:10<00:10, 135.03it/s]
Fit source/group:  48%|████▊     | 1315/2759 [00:10<00:11, 123.65it/s]
Fit source/group:  48%|████▊     | 1329/2759 [00:10<00:12, 114.74it/s]
Fit source/group:  49%|████▊     | 1342/2759 [00:10<00:12, 111.96it/s]
Fit source/group:  49%|████▉     | 1354/2759 [00:10<00:13, 101.61it/s]
Fit source/group:  50%|████▉     | 1366/2759 [00:11<00:13, 103.60it/s]
Fit source/group:  50%|████▉     | 1377/2759 [00:11<00:13, 101.55it/s]
Fit source/group:  50%|█████     | 1388/2759 [00:11<00:13, 99.39it/s]
Fit source/group:  51%|█████     | 1399/2759 [00:11<00:13, 98.01it/s]
Fit source/group:  51%|█████     | 1409/2759 [00:11<00:14, 92.65it/s]
Fit source/group:  51%|█████▏    | 1419/2759 [00:11<00:14, 92.30it/s]
Fit source/group:  52%|█████▏    | 1442/2759 [00:11<00:10, 128.53it/s]
Fit source/group:  53%|█████▎    | 1456/2759 [00:11<00:09, 131.40it/s]
Fit source/group:  54%|█████▎    | 1481/2759 [00:11<00:07, 163.66it/s]
Fit source/group:  54%|█████▍    | 1499/2759 [00:12<00:07, 166.52it/s]
Fit source/group:  55%|█████▍    | 1516/2759 [00:12<00:08, 145.25it/s]
Fit source/group:  56%|█████▌    | 1532/2759 [00:12<00:10, 117.75it/s]
Fit source/group:  56%|█████▌    | 1548/2759 [00:12<00:10, 118.46it/s]
Fit source/group:  57%|█████▋    | 1565/2759 [00:12<00:09, 127.52it/s]
Fit source/group:  57%|█████▋    | 1579/2759 [00:12<00:09, 128.51it/s]
Fit source/group:  58%|█████▊    | 1599/2759 [00:12<00:07, 146.68it/s]
Fit source/group:  59%|█████▊    | 1615/2759 [00:13<00:09, 117.07it/s]
Fit source/group:  59%|█████▉    | 1629/2759 [00:13<00:09, 122.05it/s]
Fit source/group:  60%|█████▉    | 1643/2759 [00:13<00:10, 106.18it/s]
Fit source/group:  60%|██████    | 1657/2759 [00:13<00:09, 113.27it/s]
Fit source/group:  61%|██████    | 1670/2759 [00:13<00:09, 111.39it/s]
Fit source/group:  61%|██████    | 1684/2759 [00:13<00:09, 117.59it/s]
Fit source/group:  62%|██████▏   | 1699/2759 [00:13<00:08, 125.24it/s]
Fit source/group:  62%|██████▏   | 1713/2759 [00:13<00:08, 128.09it/s]
Fit source/group:  63%|██████▎   | 1728/2759 [00:14<00:08, 123.79it/s]
Fit source/group:  63%|██████▎   | 1741/2759 [00:14<00:08, 115.78it/s]
Fit source/group:  64%|██████▎   | 1753/2759 [00:14<00:08, 112.67it/s]
Fit source/group:  64%|██████▍   | 1765/2759 [00:14<00:09, 108.19it/s]
Fit source/group:  64%|██████▍   | 1776/2759 [00:14<00:09, 105.23it/s]
Fit source/group:  65%|██████▍   | 1787/2759 [00:14<00:10, 93.78it/s]
Fit source/group:  65%|██████▌   | 1803/2759 [00:14<00:08, 110.09it/s]
Fit source/group:  66%|██████▌   | 1821/2759 [00:14<00:07, 127.17it/s]
Fit source/group:  67%|██████▋   | 1835/2759 [00:14<00:08, 113.16it/s]
Fit source/group:  67%|██████▋   | 1852/2759 [00:15<00:07, 123.82it/s]
Fit source/group:  68%|██████▊   | 1865/2759 [00:15<00:07, 117.39it/s]
Fit source/group:  68%|██████▊   | 1878/2759 [00:15<00:07, 114.04it/s]
Fit source/group:  69%|██████▊   | 1895/2759 [00:15<00:06, 124.23it/s]
Fit source/group:  69%|██████▉   | 1910/2759 [00:15<00:07, 119.56it/s]
Fit source/group:  70%|██████▉   | 1923/2759 [00:15<00:07, 108.32it/s]
Fit source/group:  70%|███████   | 1938/2759 [00:15<00:06, 117.90it/s]
Fit source/group:  71%|███████   | 1951/2759 [00:16<00:08, 99.37it/s]
Fit source/group:  72%|███████▏  | 1974/2759 [00:16<00:06, 128.70it/s]
Fit source/group:  72%|███████▏  | 1989/2759 [00:16<00:07, 108.73it/s]
Fit source/group:  73%|███████▎  | 2002/2759 [00:16<00:07, 107.08it/s]
Fit source/group:  73%|███████▎  | 2015/2759 [00:16<00:06, 112.32it/s]
Fit source/group:  74%|███████▎  | 2028/2759 [00:16<00:06, 107.23it/s]
Fit source/group:  74%|███████▍  | 2045/2759 [00:16<00:06, 118.80it/s]
Fit source/group:  75%|███████▍  | 2058/2759 [00:16<00:06, 114.94it/s]
Fit source/group:  75%|███████▌  | 2074/2759 [00:17<00:05, 119.36it/s]
Fit source/group:  76%|███████▌  | 2087/2759 [00:17<00:05, 116.05it/s]
Fit source/group:  76%|███████▌  | 2099/2759 [00:17<00:05, 111.41it/s]
Fit source/group:  77%|███████▋  | 2117/2759 [00:17<00:05, 126.59it/s]
Fit source/group:  77%|███████▋  | 2134/2759 [00:17<00:04, 137.38it/s]
Fit source/group:  78%|███████▊  | 2149/2759 [00:17<00:04, 129.09it/s]
Fit source/group:  78%|███████▊  | 2163/2759 [00:17<00:05, 114.37it/s]
Fit source/group:  79%|███████▉  | 2177/2759 [00:17<00:04, 120.65it/s]
Fit source/group:  80%|███████▉  | 2194/2759 [00:18<00:04, 132.48it/s]
Fit source/group:  80%|████████  | 2209/2759 [00:18<00:04, 135.62it/s]
Fit source/group:  81%|████████  | 2223/2759 [00:18<00:04, 107.24it/s]
Fit source/group:  81%|████████  | 2235/2759 [00:18<00:05, 104.67it/s]
Fit source/group:  81%|████████▏ | 2247/2759 [00:18<00:05, 102.27it/s]
Fit source/group:  82%|████████▏ | 2264/2759 [00:18<00:04, 118.06it/s]
Fit source/group:  83%|████████▎ | 2277/2759 [00:18<00:04, 119.33it/s]
Fit source/group:  83%|████████▎ | 2290/2759 [00:18<00:04, 116.61it/s]
Fit source/group:  84%|████████▍ | 2311/2759 [00:18<00:03, 140.58it/s]
Fit source/group:  84%|████████▍ | 2326/2759 [00:19<00:03, 140.49it/s]
Fit source/group:  85%|████████▍ | 2341/2759 [00:19<00:02, 139.82it/s]
Fit source/group:  85%|████████▌ | 2356/2759 [00:19<00:03, 111.50it/s]
Fit source/group:  86%|████████▌ | 2369/2759 [00:19<00:03, 108.15it/s]
Fit source/group:  86%|████████▋ | 2381/2759 [00:19<00:03, 106.48it/s]
Fit source/group:  87%|████████▋ | 2393/2759 [00:19<00:03, 105.45it/s]
Fit source/group:  87%|████████▋ | 2404/2759 [00:19<00:03, 100.84it/s]
Fit source/group:  88%|████████▊ | 2420/2759 [00:19<00:02, 115.87it/s]
Fit source/group:  88%|████████▊ | 2433/2759 [00:20<00:03, 96.18it/s]
Fit source/group:  89%|████████▊ | 2448/2759 [00:20<00:02, 107.70it/s]
Fit source/group:  89%|████████▉ | 2460/2759 [00:20<00:03, 97.98it/s]
Fit source/group:  90%|████████▉ | 2473/2759 [00:20<00:02, 103.14it/s]
Fit source/group:  90%|█████████ | 2486/2759 [00:20<00:02, 103.19it/s]
Fit source/group:  91%|█████████ | 2497/2759 [00:20<00:02, 92.60it/s]
Fit source/group:  91%|█████████▏| 2518/2759 [00:20<00:02, 113.60it/s]
Fit source/group:  92%|█████████▏| 2530/2759 [00:21<00:02, 96.49it/s]
Fit source/group:  92%|█████████▏| 2541/2759 [00:21<00:02, 95.36it/s]
Fit source/group:  93%|█████████▎| 2557/2759 [00:21<00:01, 109.85it/s]
Fit source/group:  93%|█████████▎| 2569/2759 [00:21<00:01, 108.50it/s]
Fit source/group:  94%|█████████▍| 2590/2759 [00:21<00:01, 125.83it/s]
Fit source/group:  95%|█████████▍| 2614/2759 [00:21<00:01, 142.10it/s]
Fit source/group:  96%|█████████▌| 2635/2759 [00:21<00:00, 158.95it/s]
Fit source/group:  96%|█████████▌| 2652/2759 [00:21<00:00, 156.35it/s]
Fit source/group:  97%|█████████▋| 2668/2759 [00:22<00:00, 156.35it/s]
Fit source/group:  97%|█████████▋| 2687/2759 [00:22<00:00, 160.08it/s]
Fit source/group:  98%|█████████▊| 2708/2759 [00:22<00:00, 162.97it/s]
Fit source/group:  99%|█████████▉| 2725/2759 [00:22<00:00, 119.01it/s]
Fit source/group:  99%|█████████▉| 2739/2759 [00:22<00:00, 123.06it/s]
Fit source/group: 100%|██████████| 2759/2759 [00:22<00:00, 121.44it/s]
WARNING: One or more fit(s) may not have converged. Please check the "flags" column in the output table. [photutils.psf.photometry]
Add model sources:   0%|          | 0/2759 [00:00<?, ?it/s]
Add model sources:  19%|█▉        | 518/2759 [00:00<00:00, 5179.68it/s]
Add model sources:  38%|███▊      | 1042/2759 [00:00<00:00, 5210.66it/s]
Add model sources:  57%|█████▋    | 1570/2759 [00:00<00:00, 5238.80it/s]
Add model sources:  77%|███████▋  | 2114/2759 [00:00<00:00, 5316.65it/s]
Add model sources:  97%|█████████▋| 2664/2759 [00:00<00:00, 5379.84it/s]
Add model sources: 100%|██████████| 2759/2759 [00:00<00:00, 5321.02it/s]

Again, we convert the pixel coordinates to sky coordinates and add them to the catalog.

wcs2 = WCS(hdr2) # NISP
radec = wcs2.all_pix2world(phot2["x_fit"],phot2["y_fit"],0)
phot2["ra_fit"] = radec[0]
phot2["dec_fit"] = radec[1]

Finally, we create the same figure as above, showing the NISP image and the residual with the (VIS-extracted) sources overlaid.

fig = plt.figure(figsize=(20,10))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.imshow(np.log10(img2), cmap="Greys", origin="lower")
ax1.plot(phot2["x_fit"], phot2["y_fit"] , "o", markersize=8 , markeredgecolor="red", fillstyle="none")

ax2.imshow(resimage2,vmin=-20*std2, vmax=20*std2, cmap="RdBu", origin="lower")
ax2.plot(phot2["x_fit"], phot2["y_fit"] , "o", markersize=8 , markeredgecolor="red", fillstyle="none")

plt.show()

../../_images/58bf79dea59d97bda9da6ffc6764773aa294520934408b955d2be2dca2b26432.png

Load Gaia Catalog#

We now load the Gaia sources at the location of the globular clusters. The goal is to compare the photometry of Gaia to the one derived above for the Euclid VIS and NISP images. This is scientifically useful, for example we can compute the colors of the stars in the Gaia optical bands and the Euclid near-IR bands. To search for Gaia sources, we use astroquery again.

We first have to elimiate the row limit for the Gaia query by setting

Gaia.ROW_LIMIT = -1

Next, we request the Gaia catalog around the position of the globular cluster. We use the same size as the cutout size.

gaia_objects = Gaia.query_object_async(coordinate=coord, radius = cutout_size/2)
print("Number of Gaia stars found: {}".format(len(gaia_objects)))

INFO: Query finished. [astroquery.utils.tap.core]
Number of Gaia stars found: 1202

We then convert the sky coordinates of the Gaia stars to (x,y) image coordinates for VIS and NISP images using the corresponding WCS. This makes it more easy to plot the Gaia sources later on the images.

wcs = WCS(hdr) # VIS
wcs2 = WCS(hdr2) # NISP
xy = wcs.all_world2pix(gaia_objects["ra"],gaia_objects["dec"],0)
xy2 = wcs2.all_world2pix(gaia_objects["ra"],gaia_objects["dec"],0)

gaia_objects["x_vis"] = xy[0]
gaia_objects["y_vis"] = xy[1]
gaia_objects["x_nisp"] = xy2[0]
gaia_objects["y_nisp"] = xy2[1]

We save the Gaia table to disk as we will later use it for the visualization in Firefly.

gaia_objects.write("./data/gaiatable.csv", format="csv", overwrite=True)

Now we can overlay the Gaia sources on the VIS and NISP images (here the x/y coordinates become handy).

fig = plt.figure(figsize=(20,10))
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.imshow(np.log10(img), cmap="Greys", origin="lower")
ax1.plot(gaia_objects["x_vis"], gaia_objects["y_vis"] , "o", markersize=8 , markeredgecolor="red", fillstyle="none")
ax1.set_title("VIS")

ax2.imshow(np.log10(img2), cmap="Greys", origin="lower")
ax2.plot(gaia_objects["x_nisp"], gaia_objects["y_nisp"] , "o", markersize=8 , markeredgecolor="red", fillstyle="none")
ax2.set_title("NISP")

plt.show()

../../_images/536c91360b21ace8c38e6f405642d77bb159028d73b12f3927a3c404148d640c.png

Match the Gaia Catalog to the VIS and NISP Catalogs#

Now, we match the Gaia source positions to the extracted sources in the VIS and NISP images.

We first define which Gaia columns to copy to the matched catalog as well as the matching distance.

gaia_keys = ["source_id", "phot_g_mean_mag", "phot_bp_mean_mag", "phot_rp_mean_mag","ra","dec","pmra","pmdec"]
matching_distance = 0.6*u.arcsecond

First match to the VIS image. We use the astropy SkyCoord() function for matching in sky coordinates.

c = SkyCoord(ra=phot["ra_fit"]*u.degree, dec=phot["dec_fit"]*u.degree )
catalog = SkyCoord(ra=gaia_objects["ra"].data*u.degree, dec=gaia_objects["dec"].data*u.degree)
idx, d2d, d3d = c.match_to_catalog_sky(catalog)

sel_matched = np.where(d2d.to(u.arcsecond) < (matching_distance))[0]
print("Gaia Sources matched to VIS: {}".format( len(sel_matched) ) )
phot["gaia_distance"] = d2d.to(u.arcsecond)

for gaia_key in gaia_keys:
    phot["gaia_{}".format(gaia_key)] = 0.0
    phot["gaia_{}".format(gaia_key)][sel_matched] = gaia_objects[gaia_key][idx[sel_matched]]

Gaia Sources matched to VIS: 1223

And then we add the NISP sources. Note that we do not have to perform matching here because by design the VIS and NISP sources are matched (spatial prior forced photometry).

phot2["gaia_distance"] = d2d.to(u.arcsecond)

for gaia_key in gaia_keys:
    phot2["gaia_{}".format(gaia_key)] = 0.0
    phot2["gaia_{}".format(gaia_key)][sel_matched] = gaia_objects[gaia_key][idx[sel_matched]]

Once matched, we can now compare the Gaia and Euclid/NISP magnitudes of the stars.

# Data
x = phot["gaia_phot_rp_mean_mag"]
y = -2.5*np.log10(phot["flux_fit"]) + hdr["ZP_STACK"]

# selection
sel_good = np.where(phot["gaia_source_id"] > 0)[0]
x = x[sel_good]
y = y[sel_good]

fig = plt.figure(figsize=(6,5))
ax1 = fig.add_subplot(1,1,1)

ax1.plot(x , y , "o", markersize=2)
minlim = np.nanmin(np.concatenate((x,y)))
maxlim = np.nanmax(np.concatenate((x,y)))

ax1.fill_between(np.asarray([minlim,maxlim]),np.asarray([minlim,maxlim])/1.5,np.asarray([minlim,maxlim])*1.5, color="gray", alpha=0.2, linewidth=0)
ax1.fill_between(np.asarray([minlim,maxlim]),np.asarray([minlim,maxlim])/1.2,np.asarray([minlim,maxlim])*1.2, color="gray", alpha=0.4, linewidth=0)
ax1.plot(np.asarray([minlim,maxlim]),np.asarray([minlim,maxlim]), ":", color="gray")

ax1.set_xlabel("Gaia Rp [mag]")
ax1.set_ylabel("I$_E$ [mag]")
plt.show()

../../_images/04b5036471dc462ded5483fef537f5c4c3f046b489c3d6258080a044fe534edb.png

Visualization with Firefly#

At the end of this Notebook, we demonstrate how we can visualize the images and catalogs created above in Firefly.

We start by initializing the Firefly client. The following line will open a new Firefly GUI in a separate tab inside the Jupyter Notebook environment. The user can drag the tab onto the currently open tab to create a “split tab”. This the user to see the code and images side-by-side.

# Uncomment when using within Jupyter Lab with jupyter_firefly_extensions installed
# fc = FireflyClient.make_lab_client()

# Uncomment for contexts other than the above
fc = FireflyClient.make_client(url="https://irsa.ipac.caltech.edu/irsaviewer")

Open your web browser to this link

In order to display in image or catalog in Firefly, it needs to be uploaded to the Firefly server. We do this here using the upload_file() function. We first upload the FITS image that we created above.

fval = fc.upload_file('./data/euclid_images_test.fits')

Once the image is uploaded we can use the show_fits() function to display it. Note that our FITS image has multiple extensions (VIS, and NISP bands). We can open them separately in new Firefly tabs by looping over the HDUs and specifying the plot ID by the extension’s name.

for hh,hdu in enumerate(hdulcutout):
    fc.show_fits(fval, MultiImageIdx=hh, plot_id=hdu.header["EXTNAME"] )

We can lock the WCS between the images (allowing the user to pan and zoom the images simultaneously) by running:

fc.align_images(lock_match=True)

{'success': True}

In the same way, we can upload a table, in this case our Gaia table. We again use upload_file() but in this case we use show_table() to show it in Firefly.

tval = fc.upload_file('./data/gaiatable.csv')
fc.show_table(tval, tbl_id = "gaiatable")

{'success': True}

Now, check out the Firefly GUI. You can zoom the images, click on sources, filter the table, display different selection, and much more!

About this Notebook#

Author: Andreas Faisst (IPAC Scientist)

Updated: 2025-03-17

Contact: the IRSA Helpdesk with questions or reporting problems.

Exploring Star Clusters in the Euclid ERO Data

Contents

Exploring Star Clusters in the Euclid ERO Data#

Learning Goals#

Introduction#

Data Volume#

Imports#

Setting up the Environment#

Search Euclid ERO Images#

Create Cutout Images#

Extract Sources and Measure their Photometry on the VIS Image#

Measure the Photometry on the NISP Images#

Load Gaia Catalog#

Match the Gaia Catalog to the VIS and NISP Catalogs#

Visualization with Firefly#

About this Notebook#