Merge pull request #96 from bsipocz/MAINT_euclid4_astropyfication

bsipocz · web-flow · commit 1986cb24517c · 2025-04-15T19:15:39.000-07:00
MAINT: cleanup euclid4 and introduced a couple of astropy features
diff --git a/tutorials/euclid_access/4_Euclid_intro_PHZ_catalog.md b/tutorials/euclid_access/4_Euclid_intro_PHZ_catalog.md
@@ -51,23 +51,24 @@ If you have questions about this notebook, please contact the [IRSA helpdesk](ht
 
 ```{code-cell} ipython3
 # Uncomment the next line to install dependencies if needed.
-# !pip install matplotlib pandas 'astropy>=5.3' 'astroquery>=0.4.10' fsspec firefly_client
+# !pip install matplotlib 'astropy>=5.3' 'astroquery>=0.4.10' fsspec firefly_client
 ```
 
 ```{code-cell} ipython3
 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 Table
+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
+from astropy.visualization import ImageNormalize, PercentileInterval, AsinhStretch, LogStretch, quantity_support
 from astropy.wcs import WCS
 
 from firefly_client import FireflyClient
@@ -90,32 +91,30 @@ coord = SkyCoord(ra, dec, unit='deg', frame='icrs')
 
 This searches specifically in the euclid_DpdMerBksMosaic "collection" which is the MER images and catalogs.
 
-```{code-cell} ipython3
-image_table = Irsa.query_sia(pos=(coord, search_radius), collection='euclid_DpdMerBksMosaic')
-```
++++
 
 ```{note}
-This table lists all MER mosaic images available in this position. These mosaics include the Euclid VIS, Y, J, H images, as well as ground-based telescopes which have been put on the same pixel scale. For more information, see the [Euclid documentation at IPAC](https://euclid.caltech.edu/page/euclid-faq-tech/).
+This table lists all MER mosaic images available in this position. These mosaics include the Euclid VIS, Y, J, H images, as well as ground-based telescopes which have been put on the same pixel scale. For more information, see the [Euclid documentation at IPAC](https://euclid.caltech.edu/page/euclid-faq-tech/). We use the ``facility`` argument below to query for Euclid images only. 
 ```
 
 ```{code-cell} ipython3
-# Convert the table to pandas dataframe
-df_im_irsa=image_table.to_pandas()
+image_table = Irsa.query_sia(pos=(coord, search_radius), collection='euclid_DpdMerBksMosaic', facility='Euclid')
 ```
 
-```{code-cell} ipython3
-df_im_euclid=df_im_irsa[ (df_im_irsa['dataproduct_subtype']=='science') &  (df_im_irsa['facility_name']=='Euclid')]
+Note that there are various image types are returned as well, we filter out the science images from these:
 
-df_im_euclid.head()
+```{code-cell} ipython3
+science_images = image_table[image_table['dataproduct_subtype'] == 'science']
+science_images
 ```
 
 Choose the VIS image and pull the filename and tileID
 
 ```{code-cell} ipython3
-filename=df_im_euclid[df_im_euclid['energy_bandpassname']=='VIS']['access_url'].to_list()[0]
+filename = science_images[science_images['energy_bandpassname'] == 'VIS']['access_url'][0]
+tileID = science_images[science_images['energy_bandpassname'] == 'VIS']['obs_id'][0][:9]
 
-tileID=re.search(r'TILE\s*(\d{9})', filename).group(1)
-print('The MER tile ID for this object is :',tileID)
+print(f'The MER tile ID for this object is : {tileID}')
 ```
 
 ## 2. Download PHZ catalog from IRSA
@@ -172,10 +171,10 @@ The _unif fluxes are: "Unified flux recomputed after correction from galactic ex
 We specify the following conditions on our search:
 - We select just the galaxies where the flux is greater than zero, to ensure the appear in all four of the Euclid MER images.
 - Select only objects in a circle (search radius selected below) around our selected RA and Dec
-- phz_classification =2 means we select only galaxies
-- Using the phz_90_int1 and phz_90_int2, we select just the galaxies where the error on the photometric redshift is less than 20%
+- `phz_classification = 2` means we select only galaxies
+- Using the `phz_90_int1` and `phz_90_int2`, we select just the galaxies where the error on the photometric redshift is less than 20%
 - Select just the galaxies between a median redshift of 1.4 and 1.6
-- We search just a 3 arcminute box around an RA and Dec
+- We search just a 5 arcminute box around an RA and Dec
 
 +++
 
@@ -190,53 +189,47 @@ im_cutout= 5 * u.arcmin
 ra_cutout = 267.8
 dec_cutout =  66
 
-coords_cutout = SkyCoord(ra_cutout, dec_cutout, unit=(u.deg, u.deg), frame='icrs')
-##########################################################################
-im_cutout_deg=im_cutout.to(u.deg).value
-
-## Create ra and dec bounds for the box
-ra0=ra_cutout-im_cutout_deg/2
-ra1=ra_cutout+im_cutout_deg/2
-
-dec0=dec_cutout-im_cutout_deg/2
-dec1=dec_cutout+im_cutout_deg/2
+coords_cutout = SkyCoord(ra_cutout, dec_cutout, unit='deg', frame='icrs')
+size_cutout = im_cutout.to(u.deg).value
 ```
 
 ```{code-cell} ipython3
-adql = f"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 \
-FROM {table_mer} AS mer \
-JOIN {table_phz} as phz \
-ON mer.object_id = phz.object_id \
-WHERE 1 = CONTAINS(POINT('ICRS', mer.ra, mer.dec), CIRCLE('ICRS', {ra_cutout}, {dec_cutout}, {im_cutout_deg/2})) \
-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 \
-"
-adql
+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 = Irsa.query_tap(adql)
-
+result_galaxies = Irsa.query_tap(adql).to_table()
+result_galaxies[:5]
+```
 
-## Convert table to pandas dataframe
-df_g_irsa = result.to_table().to_pandas()
+```{warning}
+Note that we use `to_table` above rather than `to_qtable`. While astropy's `QTable` is more powerful than its `Table`, as it e.g. handles the column units properly, we cannot use it here due to a known bug; it mishandles the large integer numbers in the `object_id` column and recast them as float during which process some precision is being lost. 
 
-# Display first few rows
-df_g_irsa.head()
+Once the bug is fixed, we plan to update the code in this notebook and simplify some of the approaches below.
 ```
 
-```{code-cell} ipython3
-len(df_g_irsa)
-```
++++
 
 ## 3. Read in a cutout of the MER image from IRSA directly
 
++++
+
+Due to the large field of view of the MER mosaic, let's cut out a smaller section (5'x5') of the MER mosaic to inspect the image.
+
 ```{code-cell} ipython3
 ## Use fsspec to interact with the fits file without downloading the full file
 hdu = fits.open(filename, use_fsspec=True)
@@ -245,152 +238,136 @@ hdu = fits.open(filename, use_fsspec=True)
 header = hdu[0].header
 
 ## Read in the cutout of the image that you want
-cutout_data = Cutout2D(hdu[0].section, position=coords_cutout, size=im_cutout, wcs=WCS(hdu[0].header))
-
-## Define a new fits file based on this smaller cutout, with accurate WCS based on the cutout size
-new_hdu = fits.PrimaryHDU(data=cutout_data.data, header=header)
-new_hdu.header.update(cutout_data.wcs.to_header())
-```
-
-```{code-cell} ipython3
-im_mer_irsa = new_hdu.data
+cutout_image = Cutout2D(hdu[0].section, position=coords_cutout, size=im_cutout, wcs=WCS(header))
 ```
 
 ```{code-cell} ipython3
-im_mer_irsa.shape
+cutout_image.data.shape
 ```
 
 ```{code-cell} ipython3
-norm = ImageNormalize(im_mer_irsa, interval=PercentileInterval(99.9), stretch=AsinhStretch())
-plt.imshow(im_mer_irsa, cmap='gray', origin='lower', norm=norm)
+norm = ImageNormalize(cutout_image.data, interval=PercentileInterval(99.9), stretch=AsinhStretch())
+_ = plt.imshow(cutout_image.data, cmap='gray', origin='lower', norm=norm)
 ```
 
 ## 4. Overplot the catalog on the MER mosaic image
 
-```{code-cell} ipython3
-## Use the WCS package to extract the coordinates from the header of the image
-head_mer_irsa=new_hdu.header
-wcs_irsa=WCS(head_mer_irsa) # VIS
-```
++++
 
-```{code-cell} ipython3
-## Convert the catalog to match the pixels in the image
-xy_irsa = wcs_irsa.all_world2pix(df_g_irsa["ra"],df_g_irsa["dec"],0)
+```{tip}
+We can rely on astropy's WCSAxes framework for making plots of Astronomical data in Matplotlib. Please note the usage of `projection` and `transform` arguments in the code example below.
+
+For more info, please visit the [WCSAxes documentation](https://docs.astropy.org/en/stable/visualization/wcsaxes/index.html).
 ```
 
 ```{code-cell} ipython3
-df_g_irsa['x_pix']=xy_irsa[0]
-df_g_irsa['y_pix']=xy_irsa[1]
-```
+ax = plt.subplot(projection=cutout_image.wcs)
 
-Due to the large field of view of the MER mosaic, let's cut out a smaller section (3'x3')of the MER mosaic to inspect the image
+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')
+```
 
-Plot MER catalog sources on the Euclid VIS image 3'x3' cutout
+## 5. Pull the spectra on the top brightest source based on object ID
 
 ```{code-cell} ipython3
-plt.imshow(im_mer_irsa, cmap='gray', origin='lower', norm=ImageNormalize(im_mer_irsa, interval=PercentileInterval(99.9), stretch=LogStretch()))
-colorbar = plt.colorbar()
-plt.scatter(df_g_irsa['x_pix'], df_g_irsa['y_pix'], s=36, facecolors='none', edgecolors='red')
-
-plt.title('Galaxies between z = 1.4 and 1.6')
-plt.show()
+result_galaxies.sort(keys='flux_h_unif', reverse=True)
 ```
 
-Pull the spectra on the top brightest source based on object ID
-
 ```{code-cell} ipython3
-df_g_irsa_sort=df_g_irsa.sort_values(by='flux_h_unif',ascending=False)
+result_galaxies[:3]
 ```
 
+Let's pick one of these galaxies. Note that the table has been sorted above, we can use the same index here and below to access the data for this particular galaxy. 
+
 ```{code-cell} ipython3
-df_g_irsa_sort.iloc[0:3]
+index = 2
+
+obj_id = result_galaxies['object_id'][index]
+redshift = result_galaxies['phz_median'][index]
 ```
 
-```{code-cell} ipython3
-obj_id=df_g_irsa_sort['object_id'].iloc[2]
-redshift = df_g_irsa_sort['phz_median'].iloc[2]
+We will use TAP and an ASQL query to find the spectral data for this particular galaxy.
 
-## Pull the data on these objects
-adql_object = f"SELECT * \
-FROM {table_1dspectra} \
-WHERE objectid = {obj_id}"
+```{code-cell} ipython3
+adql_object = f"SELECT * FROM {table_1dspectra} WHERE objectid = {obj_id}"
 
 ## Pull the data on this particular galaxy
-result2 = Irsa.query_tap(adql_object)
-df2=result2.to_table().to_pandas()
-df2
+result_spectra = Irsa.query_tap(adql_object).to_table()
+result_spectra
 ```
 
-Pull out the file name from the ``result`` table:
+Pull out the file name from the ``result_spectra`` table:
 
 ```{code-cell} ipython3
-file_uri = urllib.parse.urljoin(Irsa.tap_url, result2['uri'][0])
+file_uri = urllib.parse.urljoin(Irsa.tap_url, result_spectra['uri'][0])
 file_uri
 ```
 
 ```{code-cell} ipython3
 with fits.open(file_uri) as hdul:
-    hdu = hdul[df2['hdu'].iloc[0]]
-    dat = Table.read(hdu, format='fits', hdu=1)
-    df_obj_irsa = dat.to_pandas()
+    spectrum = QTable.read(hdul[result_spectra['hdu'][0]], format='fits')
+    spectrum_header = hdul[result_spectra['hdu'][0]].header
 ```
 
 ### Now the data are read in, plot the spectrum
 
-Divide by 10000 to convert from Angstrom to micron
+```{tip}
+As we use astropy.visualization’s quantity_support, matplotlib automatically picks up the axis units from the quantitites we plot.
+```
 
 ```{code-cell} ipython3
-plt.plot(df_obj_irsa['WAVELENGTH']/10000., df_obj_irsa['SIGNAL'])
+quantity_support()
+```
+
+```{code-cell} ipython3
+plt.plot(spectrum['WAVELENGTH'].to(u.micron), spectrum['SIGNAL'])
 
-plt.xlabel('Wavelength (microns)')
-plt.ylabel('Flux (erg / (s cm2))')
 plt.xlim(1.25, 1.85)
-plt.ylim(-0.5,0.5)
-plt.title('Object ID is '+str(obj_id)+'with phz_median='+str(redshift))
+plt.ylim(-0.5, 0.5)
+_ = plt.title(f"Object {obj_id} with phz_median={redshift}")
 ```
 
-Let's cut out a very small patch of the MER image to see what this galaxy looks like
+Let's cut out a very small patch of the MER image to see what this galaxy looks like. Remember that we sorted the table above, so can reuse the same index to pick up the coordinates for the galaxy. Otherwise we could filter on the object ID.
 
 ```{code-cell} ipython3
-## How large do you want the image cutout to be?
-im_cutout= 2.0 * u.arcsec
-
-## Use the ra and dec of the galaxy
-ra = df_g_irsa[df_g_irsa['object_id']==obj_id]['ra'].iloc[0]
-dec =  df_g_irsa[df_g_irsa['object_id']==obj_id]['dec'].iloc[0]
+result_galaxies[index]
+```
 
-coords_cutout = SkyCoord(ra, dec, unit=(u.deg,u.deg), frame='icrs')
+```{code-cell} ipython3
+## How large do you want the image cutout to be?
+size_galaxy_cutout = 2.0 * u.arcsec
 ```
 
-Use ``fsspec`` to obtain a cutout of the fits file
+Use the `ra` and `dec` columns for the galaxy to create a `SkyCoord`.
 
 ```{code-cell} ipython3
-hdu = fits.open(filename, use_fsspec=True)
-
-header = hdu[0].header
+coords_galaxy = SkyCoord(result_galaxies['ra'][index], result_galaxies['dec'][index], unit='deg')
 ```
 
 ```{code-cell} ipython3
-## Read in the cutout of the image that you want
-cutout_data = Cutout2D(hdu[0].section, position=coords_cutout, size=im_cutout, wcs=WCS(hdu[0].header))
+coords_galaxy
 ```
 
+We haven't closed the image file above, so use `Cutout2D` again to cut out a section around the galaxy.
+
 ```{code-cell} ipython3
-new_hdu = fits.PrimaryHDU(data=cutout_data.data, header=header)
-new_hdu.header.update(cutout_data.wcs.to_header())
+cutout_galaxy = Cutout2D(hdu[0].section, position=coords_galaxy, size=size_galaxy_cutout, wcs=WCS(header))
 ```
 
+Plot to show the cutout on the galaxy
+
 ```{code-cell} ipython3
-## Plot a quick simple plot to show the cutout on the galaxy
+ax = plt.subplot(projection=cutout_galaxy.wcs)
 
-plt.imshow(new_hdu.data, cmap='gray', origin='lower',
-           norm=ImageNormalize(new_hdu.data, interval=PercentileInterval(99.9), stretch=AsinhStretch()))
-colorbar = plt.colorbar()
+ax.imshow(cutout_galaxy.data, cmap='gray', origin='lower', 
+          norm=ImageNormalize(cutout_galaxy.data, interval=PercentileInterval(99.9), stretch=AsinhStretch()))
 ```
 
-## 5. Load the image on Firefly to be able to interact with the data directly
+## 6. Load the image on Firefly to be able to interact with the data directly
 
 +++
 
@@ -423,7 +400,7 @@ fc.align_images(lock_match=True)
 
 ```{code-cell} ipython3
 csv_path = os.path.join(download_path, "mer_df.csv")
-df_g_irsa.to_csv(csv_path, index=False)
+result_galaxies.write(csv_path, format="csv")
 ```
 
 ### Upload the CSV table to Firefly and display as an overlay on the FITS image
@@ -439,6 +416,6 @@ fc.show_table(uploaded_table)
 
 **Author**: Tiffany Meshkat, Anahita Alavi, Anastasia Laity, Andreas Faisst, Brigitta Sipőcz, Dan Masters, Harry Teplitz, Jaladh Singhal, Shoubaneh Hemmati, Vandana Desai
 
-**Updated**: 2025-03-31
+**Updated**: 2025-04-10
 
 **Contact:** [the IRSA Helpdesk](https://irsa.ipac.caltech.edu/docs/help_desk.html) with questions or reporting problems.
