Clean up comments and formatting

spectroscopy/code_src/desi_functions.py

@@ -31,7 +31,6 @@ def DESIBOSS_get_spec(sample_table, search_radius_arcsec):
 
     # Set up client
     client = SparclClient()
-    # print(client.all_datasets) # print data sets
 
     # Initialize multi-index object:
     df_spec = MultiIndexDFObject()
@@ -40,20 +39,17 @@ def DESIBOSS_get_spec(sample_table, search_radius_arcsec):
 
         # Search
         data_releases = ['DESI-EDR', 'BOSS-DR16']
-        # data_releases = ['DESI-EDR','BOSS-DR16','SDSS-DR16'] # we want to use DR17 directly using SDSS query
 
         search_coords = stab["coord"]
         dra = (search_radius_arcsec*u.arcsec).to(u.degree)
         ddec = (search_radius_arcsec*u.arcsec).to(u.degree)
         out = ['sparcl_id', 'ra', 'dec', 'redshift', 'spectype', 'data_release', 'redshift_err']
         cons = {'spectype': ['GALAXY', 'STAR', 'QSO'],
                 'data_release': data_releases,
-                # 'redshift': [0.5, 0.9],
                 'ra': [search_coords.ra.deg-dra.value, search_coords.ra.deg+dra.value],
                 'dec': [search_coords.dec.deg-ddec.value, search_coords.dec.deg+ddec.value]
                 }
         found_I = client.find(outfields=out, constraints=cons, limit=20)  # search
-        # print(found_I)
         if len(found_I.records) == 0:
             continue
 

spectroscopy/code_src/herschel_functions.py

@@ -80,16 +80,22 @@ def Herschel_get_spec(sample_table, search_radius_arcsec, datadir,
         # each instrument and object..
 
         for instrument_name in ['PACS', 'SPIRE']:
-            querystring = "select observation_id from hsa.v_active_observation join hsa.instrument using (instrument_oid) where contains(point('ICRS', hsa.v_active_observation.ra, hsa.v_active_observation.dec), circle('ICRS', "+str(
-                search_coords.ra.deg)+", " + str(search_coords.dec.deg) + ", " + str(search_radius_arcsec) + "))=1 and hsa.instrument.instrument_name='"+str(instrument_name)+"'"
+            querystring = (
+                "select observation_id from hsa.v_active_observation join hsa.instrument using "
+                "(instrument_oid) where contains(point('ICRS', hsa.v_active_observation.ra, "
+                "hsa.v_active_observation.dec), circle('ICRS', " + str(search_coords.ra.deg)
+                + ", " + str(search_coords.dec.deg) + ", " + str(search_radius_arcsec)
+                + "))=1 and hsa.instrument.instrument_name='" + str(instrument_name) + "'")
             objectid_table = HSA.query_hsa_tap(querystring)
 
             # download_data only accepts one observation_id so we need to loop over each observation_id
             for tab_id in range(len(objectid_table)):
                 observation_id = str(objectid_table[tab_id]['observation_id'])
                 try:
                     HSA.download_data(observation_id=observation_id, retrieval_type='OBSERVATION',
-                                      instrument_name=instrument_name, product_level="LEVEL2, LEVEL_2_5, LEVEL_3", download_dir=datadir)
+                                      instrument_name=instrument_name,
+                                      product_level="LEVEL2, LEVEL_2_5, LEVEL_3",
+                                      download_dir=datadir)
 
                     # ok, now we have the tar files, need to read those into the right data structure
                     # first untar
@@ -146,7 +152,8 @@ def Herschel_get_spec(sample_table, search_radius_arcsec, datadir,
                                                              mission=["Herschel"],
                                                              instrument=[instrument_name],
                                                              filter=[df["band"][0]],
-                                                             )).set_index(["objectid", "label", "filter", "mission"])
+                                                             ))
+                                dfsingle = dfsingle.set_index(["objectid", "label", "filter", "mission"])
 
                                 df_spec.append(dfsingle)
 

spectroscopy/code_src/mast_functions.py

@@ -96,13 +96,11 @@ def JWST_get_spec_helper(sample_table, search_radius_arcsec, datadir, verbose,
 
         # Query results
         search_coords = stab["coord"]
-        query_results = Observations.query_criteria(coordinates=search_coords, radius=search_radius_arcsec * u.arcsec,
-                                                    dataproduct_type=["spectrum"], obs_collection=["JWST"], intentType="science", calib_level=[3, 4],
-                                                    # instrument_name=['NIRSPEC/MSA', 'NIRSPEC/SLIT', 'NIRCAM/GRISM', 'NIRISS/WFSS'],
-                                                    instrument_name=['NIRSPEC/MSA', 'NIRSPEC/SLIT'],
-                                                    # filters=['CLEAR;PRISM','F070LP;G140M'],
-                                                    dataRights=['PUBLIC']
-                                                    )
+        query_results = Observations.query_criteria(
+            coordinates=search_coords, radius=search_radius_arcsec * u.arcsec,
+            dataproduct_type=["spectrum"], obs_collection=["JWST"], intentType="science",
+            calib_level=[3, 4], instrument_name=['NIRSPEC/MSA', 'NIRSPEC/SLIT'],
+            dataRights=['PUBLIC'])
         print("Number of search results: {}".format(len(query_results)))
 
         if len(query_results) == 0:
@@ -113,25 +111,18 @@ def JWST_get_spec_helper(sample_table, search_radius_arcsec, datadir, verbose,
         data_products_list = Observations.get_product_list(query_results)
 
         # Filter
-        data_products_list_filter = Observations.filter_products(data_products_list,
-                                                                 productType=["SCIENCE"],
-                                                                 # filters=['CLEAR;PRISM','F070LP;G140M'],
-                                                                 # filters=['CLEAR;PRISM'],
-                                                                 extension="fits",
-                                                                 # only fully reduced or contributed
-                                                                 calib_level=[3, 4],
-                                                                 productSubGroupDescription=[
-                                                                     "X1D"],  # only 1D spectra
-                                                                 # only public data
-                                                                 dataRights=['PUBLIC']
-                                                                 )
+        data_products_list_filter = Observations.filter_products(
+            data_products_list, productType=["SCIENCE"], extension="fits",
+            calib_level=[3, 4],  # only fully reduced or contributed
+            productSubGroupDescription=["X1D"],  # only 1D spectra
+            dataRights=['PUBLIC'])  # only public data
         print("Number of files to download: {}".format(len(data_products_list_filter)))
 
         if len(data_products_list_filter) == 0:
             print("Nothing to download for source {}.".format(stab["label"]))
             continue
 
-        # Download (supress output)
+        # Download (suppress output)
         trap = io.StringIO()
         with redirect_stdout(trap):
             download_results = Observations.download_products(
@@ -152,7 +143,6 @@ def JWST_get_spec_helper(sample_table, search_radius_arcsec, datadir, verbose,
                                  data_products_list_filter["obsID"][jj])[0]
             tmp = query_results[idx_cross][keys]
             tab.add_row(list(tmp[0]) + [data_products_list_filter["productFilename"][jj]])
-        # print("number in table {}".format(len(tab)))
 
         # Create multi-index object
         for jj in range(len(tab)):
@@ -167,18 +157,13 @@ def JWST_get_spec_helper(sample_table, search_radius_arcsec, datadir, verbose,
             filepath = download_results["Local Path"][file_idx[0]]
             spec1d = Table.read(filepath, hdu=1)
 
-            # print(filepath)
-            # spec1d = Spectrum1D.read(filepath)
-
             dfsingle = pd.DataFrame(dict(
                 wave=[spec1d["WAVELENGTH"].data * spec1d["WAVELENGTH"].unit],
                 flux=[spec1d["FLUX"].data * spec1d["FLUX"].unit],
                 err=[spec1d["FLUX_ERROR"].data *
                      spec1d["FLUX_ERROR"].unit],
                 label=[stab["label"]],
                 objectid=[stab["objectid"]],
-                # objID=[tab["objID"][jj]], # REMOVE
-                # obsid=[tab["obsid"][jj]],
                 mission=[tab["obs_collection"][jj]],
                 instrument=[tab["instrument_name"][jj]],
                 filter=[tab["filters"][jj]],
@@ -233,7 +218,6 @@ def JWST_group_spectra(df, verbose, quickplot):
                 print("Processing {}: ".format(filt), end=" ")
 
             sel = np.where((tab["filter"] == filt) & (tab["label"] == obj))[0]
-            # tab_sel = tab.copy()[sel]
             tab_sel = tab.iloc[sel]
             if verbose:
                 print("Number of items: {}".format(len(sel)), end=" | ")
@@ -268,13 +252,12 @@ def JWST_group_spectra(df, verbose, quickplot):
                 u.erg / u.second / (u.centimeter**2) / u.angstrom)
 
             # Add to data frame
-            dfsingle = pd.DataFrame(dict(wave=[wave_grid.to(u.micrometer)], flux=[fluxes_stack_cgs], err=[np.repeat(0, len(fluxes_stack_cgs))],
-                                         label=[tab_sel["label"].iloc[0]],
-                                         objectid=[tab_sel["objectid"].iloc[0]],
-                                         mission=[tab_sel["mission"].iloc[0]],
-                                         instrument=[tab_sel["instrument"].iloc[0]],
-                                         filter=[tab_sel["filter"].iloc[0]],
-                                         )).set_index(["objectid", "label", "filter", "mission"])
+            dfsingle = pd.DataFrame(dict(
+                wave=[wave_grid.to(u.micrometer)], flux=[fluxes_stack_cgs],
+                err=[np.repeat(0, len(fluxes_stack_cgs))], label=[tab_sel["label"].iloc[0]],
+                objectid=[tab_sel["objectid"].iloc[0]], mission=[tab_sel["mission"].iloc[0]],
+                instrument=[tab_sel["instrument"].iloc[0]], filter=[tab_sel["filter"].iloc[0]]))
+            dfsingle = dfsingle.set_index(["objectid", "label", "filter", "mission"])
             df_spec.append(dfsingle)
 
             # Quick plot
@@ -328,9 +311,10 @@ def HST_get_spec(sample_table, search_radius_arcsec, datadir, verbose,
 
         # Query results
         search_coords = stab["coord"]
-        query_results = Observations.query_criteria(coordinates=search_coords, radius=search_radius_arcsec * u.arcsec,
-                                                    dataproduct_type=["spectrum"], obs_collection=["HST"], intentType="science", calib_level=[3, 4],
-                                                    )
+        query_results = Observations.query_criteria(
+            coordinates=search_coords, radius=search_radius_arcsec * u.arcsec,
+            dataproduct_type=["spectrum"], obs_collection=["HST"], intentType="science",
+            calib_level=[3, 4],)
         print("Number of search results: {}".format(len(query_results)))
 
         if len(query_results) == 0:
@@ -341,14 +325,10 @@ def HST_get_spec(sample_table, search_radius_arcsec, datadir, verbose,
         data_products_list = Observations.get_product_list(query_results)
 
         # Filter
-        data_products_list_filter = Observations.filter_products(data_products_list,
-                                                                 productType=["SCIENCE"],
-                                                                 extension="fits",
-                                                                 # only fully reduced or contributed
-                                                                 calib_level=[3, 4],
-                                                                 productSubGroupDescription=[
-                                                                     "SX1"]  # only 1D spectra
-                                                                 )
+        data_products_list_filter = Observations.filter_products(
+            data_products_list, productType=["SCIENCE"], extension="fits",
+            calib_level=[3, 4],  # only fully reduced or contributed
+            productSubGroupDescription=["SX1"])  # only 1D spectra
         print("Number of files to download: {}".format(len(data_products_list_filter)))
 
         if len(data_products_list_filter) == 0:
@@ -387,7 +367,6 @@ def HST_get_spec(sample_table, search_radius_arcsec, datadir, verbose,
 
             # open spectrum
             filepath = download_results["Local Path"][file_idx[0]]
-            # print(filepath)
             spec1d = Spectrum1D.read(filepath)
 
             # Note: this should be in erg/s/cm2/A and any wavelength unit.
@@ -396,8 +375,6 @@ def HST_get_spec(sample_table, search_radius_arcsec, datadir, verbose,
                     spec1d.uncertainty.array * spec1d.uncertainty.unit],
                 label=[stab["label"]],
                 objectid=[stab["objectid"]],
-                # objID=[tab["objID"][jj]],
-                # obsid=[tab["obsid"][jj]],
                 mission=[tab["obs_collection"][jj]],
                 instrument=[tab["instrument_name"][jj]],
                 filter=[tab["filters"][jj]],

spectroscopy/code_src/plot_functions.py

@@ -23,7 +23,6 @@
 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
@@ -112,9 +111,8 @@ def create_figures(df_spec, bin_factor, show_nbr_figures, save_output):
         else:
             LOGX = False
 
-        for ff, (filt, filt_df) in enumerate(singleobj_df.groupby('filter')):
+        for filt, filt_df in singleobj_df.groupby('filter'):
 
-            # print("{} entries for a object {} and filter {}".format(len(filt_df.flux), objectid , filt))
             for ii in range(len(filt_df.flux)):
 
                 # get data
@@ -128,10 +126,9 @@ def create_figures(df_spec, bin_factor, show_nbr_figures, save_output):
                 flux = flux[mask]
                 err = err[mask]
 
-                # ax1.plot(wave , flux , "-"  , label="{} ({})".format(filt, filt_df.reset_index().mission[ii]) )
                 wave_bin, flux_bin, _ = bin_spectra(wave, flux, bin_factor=bin_factor)
 
-                # do some more clearning (mainly to remove some very low values)
+                # do some more cleaning (mainly to remove some very low values)
                 selnotnan = np.where(~np.isnan(flux_bin))[0]
                 flux_bin = flux_bin[selnotnan]
                 wave_bin = wave_bin[selnotnan]
@@ -140,8 +137,10 @@ def create_figures(df_spec, bin_factor, show_nbr_figures, save_output):
                 wave_bin = wave_bin[~clip_mask]
                 flux_bin = flux_bin[~clip_mask]
 
-                ax1.step(wave_bin.to(u.micrometer), flux_bin.to(u.erg / u.second / (u.centimeter**2) / u.angstrom), "-",
-                         label="{} ({})".format(filt, filt_df.reset_index().mission[ii]), where="mid")
+                ax1.step(wave_bin.to(u.micrometer),
+                         flux_bin.to(u.erg / u.second / (u.centimeter**2) / u.angstrom), "-",
+                         label="{} ({})".format(filt, filt_df.reset_index().mission[ii]),
+                         where="mid")
 
         ax1.set_title(this_label)
         if LOGX:

spectroscopy/code_src/sample_selection.py

@@ -1,8 +1,4 @@
 from astropy.table import Table, join, join_skycoord, unique
-# from astroquery.ipac.ned import Ned
-# from astroquery.sdss import SDSS
-# from astroquery.simbad import Simbad
-# from astroquery.vizier import Vizier
 
 
 def clean_sample(coords_list, labels_list, precision, verbose=1):

spectroscopy/spectra_generator.md

@@ -6,9 +6,9 @@ jupytext:
     format_version: 0.13
     jupytext_version: 1.16.4
 kernelspec:
-  display_name: science_demo
+  display_name: notebook
   language: python
-  name: conda-env-science_demo-py
+  name: python3
 ---
 
 # Extract Multi-Wavelength Spectroscopy from Archival Data
@@ -81,8 +81,8 @@ The ones with an asterisk (*) are the challenging ones.
 • ...
 ## Runtime
 
-As of 2024 August, this notebook takes ~330s to run to completion on Fornax using the 'Astrophysics
-Default Image' and the 'Large' server with 16GB RAM/ 4CPU.
+As of 2024 August, this notebook takes about 3 minutes to run to completion on Fornax using
+Server Type: 'Standard - 8GB RAM/4 CPU' and Environment: 'Default Astrophysics' (image).
 
 ## Authors:
 Andreas Faisst, Jessica Krick, Shoubaneh Hemmati, Troy Raen, Brigitta Sipőcz, David Shupe