Euclid Q1: PHZ Catalogs

# Uncomment the next line to install dependencies if needed.
# !pip install matplotlib 'astropy>=5.3' 'astroquery>=0.4.11' fsspec firefly_client

import os
import re
import io
import urllib

import numpy as np
import matplotlib.pyplot as plt

from astropy.coordinates import SkyCoord
from astropy.io import fits
from astropy.nddata import Cutout2D
from astropy.table import QTable
from astropy import units as u
from astropy.utils.data import conf, download_file
from astropy.visualization import ImageNormalize, PercentileInterval, AsinhStretch, LogStretch, quantity_support
from astropy.wcs import WCS

from firefly_client import FireflyClient
from astroquery.ipac.irsa import Irsa

# The Euclid spectrum files are large and the time it takes to read
# them can exceed astropy's default timeout limit. Increase it.
conf.remote_timeout = 120

ra = 268
dec = 66
search_radius= 10 * u.arcsec

coord = SkyCoord(ra, dec, unit='deg', frame='icrs')

image_table = Irsa.query_sia(pos=(coord, search_radius), collection='euclid_DpdMerBksMosaic', facility='Euclid')

science_images = image_table[image_table['dataproduct_subtype'] == 'science']
science_images

filename = science_images[science_images['energy_bandpassname'] == 'VIS']['access_url'][0]
tileID = science_images[science_images['energy_bandpassname'] == 'VIS']['obs_id'][0][:9]

print(f'The MER tile ID for this object is : {tileID}')

Irsa.list_catalogs(filter='euclid', include_metadata_tables=True)

table_mer = 'euclid_q1_mer_catalogue'
table_phz = 'euclid_q1_phz_photo_z'
table_1dspectra = 'euclid.objectid_spectrafile_association_q1'

columns_info = Irsa.list_columns(catalog=table_phz)
print(len(columns_info))

# Full list of columns and their description
columns_info

######################## User defined section ############################
# Set the image cutout size
im_cutout= 5 * u.arcmin

# Set the center of the cutout
ra_cutout = 267.8
dec_cutout =  66

coords_cutout = SkyCoord(ra_cutout, dec_cutout, unit='deg', frame='icrs')
size_cutout = im_cutout.to(u.deg).value

adql = ("SELECT DISTINCT mer.object_id, mer.ra, mer.dec, "
        "phz.flux_vis_unif, phz.flux_y_unif, phz.flux_j_unif, phz.flux_h_unif, "
        "phz.phz_classification, phz.phz_median, phz.phz_90_int1, phz.phz_90_int2 "
        f"FROM {table_mer} AS mer "
        f"JOIN {table_phz} as phz "
        "ON mer.object_id = phz.object_id "
        "WHERE 1 = CONTAINS(POINT('ICRS', mer.ra, mer.dec), "
                            f"BOX('ICRS', {ra_cutout}, {dec_cutout}, {size_cutout}, {size_cutout})) "
        "AND phz.flux_vis_unif> 0 "
        "AND  phz.flux_y_unif > 0 "
        "AND phz.flux_j_unif > 0 "
        "AND phz.flux_h_unif > 0 "
        "AND phz.phz_classification = 2 "
        "AND ((phz.phz_90_int2 - phz.phz_90_int1) / (1 + phz.phz_median)) < 0.20 "
        "AND phz.phz_median BETWEEN 1.4 AND 1.6")


# Use TAP with this ADQL string
result_galaxies = Irsa.query_tap(adql).to_table()
result_galaxies[:5]

# Use fsspec to interact with the fits file without downloading the full file
hdu = fits.open(filename, use_fsspec=True)

# Store the header
header = hdu[0].header

# Read in the cutout of the image that you want
cutout_image = Cutout2D(hdu[0].section, position=coords_cutout, size=im_cutout, wcs=WCS(header))

cutout_image.data.shape

norm = ImageNormalize(cutout_image.data, interval=PercentileInterval(99.9), stretch=AsinhStretch())
_ = plt.imshow(cutout_image.data, cmap='gray', origin='lower', norm=norm)

ax = plt.subplot(projection=cutout_image.wcs)

ax.imshow(cutout_image.data, cmap='gray', origin='lower',
          norm=ImageNormalize(cutout_image.data, interval=PercentileInterval(99.9), stretch=LogStretch()))
plt.scatter(result_galaxies['ra'], result_galaxies['dec'], s=36, facecolors='none', edgecolors='red',
            transform=ax.get_transform('world'))

_ = plt.title('Galaxies between z = 1.4 and 1.6')

result_galaxies.sort(keys='flux_h_unif', reverse=True)

result_galaxies[:3]

index = 2

obj_id = result_galaxies['object_id'][index]
redshift = result_galaxies['phz_median'][index]

adql_object = f"SELECT * FROM {table_1dspectra} WHERE objectid = {obj_id}"

# Pull the data on this particular galaxy
result_spectra = Irsa.query_tap(adql_object).to_table()
result_spectra

spectrum_path = f"https://irsa.ipac.caltech.edu/{result_spectra['path'][0]}"
spectrum_path

spectrum = QTable.read(spectrum_path)

quantity_support()

plt.plot(spectrum['WAVELENGTH'].to(u.micron), spectrum['SIGNAL'])

plt.xlim(1.25, 1.85)
plt.ylim(-0.5, 0.5)
_ = plt.title(f"Object {obj_id} with phz_median={redshift}")

result_galaxies[index]

# Set the image cutout size for the selected galaxy
size_galaxy_cutout = 2.0 * u.arcsec

coords_galaxy = SkyCoord(result_galaxies['ra'][index], result_galaxies['dec'][index], unit='deg')

coords_galaxy

cutout_galaxy = Cutout2D(hdu[0].section, position=coords_galaxy, size=size_galaxy_cutout, wcs=WCS(header))

ax = plt.subplot(projection=cutout_galaxy.wcs)

ax.imshow(cutout_galaxy.data, cmap='gray', origin='lower',
          norm=ImageNormalize(cutout_galaxy.data, interval=PercentileInterval(99.9), stretch=AsinhStretch()))

fc = FireflyClient.make_client('https://irsa.ipac.caltech.edu/irsaviewer')

fc.show_fits_image(filename)

tbl_stream = io.BytesIO()
result_galaxies.write(tbl_stream, format='votable')

fc.show_table(tbl_stream)