Euclid Q1: PHZ catalogs

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

import os
import re
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 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

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')

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 ############################
## How large do you want the image cutout to be?
im_cutout= 5 * u.arcmin

## What is 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/np.cos(coords_cutout.dec)}, {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]

## How large do you want the image cutout to be?
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()))

download_path = "data"
if os.path.exists(download_path):
    print("Output directory already created.")
else:
    print("Creating data directory.")
    os.mkdir(download_path)

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

fc.show_fits(url=filename)

fc.align_images(lock_match=True)

csv_path = os.path.join(download_path, "mer_df.csv")
result_galaxies.write(csv_path, format="csv")

uploaded_table = fc.upload_file(csv_path)
print(f"Uploaded Table URL: {uploaded_table}")

fc.show_table(uploaded_table)