diff --git a/spectroscopy/code_src/desi_functions.py b/spectroscopy/code_src/desi_functions.py index 94939ddd..f946a6b9 100644 --- a/spectroscopy/code_src/desi_functions.py +++ b/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,7 +39,6 @@ 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) @@ -48,12 +46,10 @@ def DESIBOSS_get_spec(sample_table, search_radius_arcsec): 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 diff --git a/spectroscopy/code_src/herschel_functions.py b/spectroscopy/code_src/herschel_functions.py index 1ead4f2a..aec4e679 100644 --- a/spectroscopy/code_src/herschel_functions.py +++ b/spectroscopy/code_src/herschel_functions.py @@ -80,8 +80,12 @@ 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 @@ -89,7 +93,9 @@ def Herschel_get_spec(sample_table, search_radius_arcsec, datadir, 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) diff --git a/spectroscopy/code_src/mast_functions.py b/spectroscopy/code_src/mast_functions.py index 664008f4..94aa7839 100644 --- a/spectroscopy/code_src/mast_functions.py +++ b/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,9 +157,6 @@ 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], @@ -177,8 +164,6 @@ def JWST_get_spec_helper(sample_table, search_radius_arcsec, datadir, verbose, 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]], diff --git a/spectroscopy/code_src/plot_functions.py b/spectroscopy/code_src/plot_functions.py index fad19dcb..e9876839 100644 --- a/spectroscopy/code_src/plot_functions.py +++ b/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: diff --git a/spectroscopy/code_src/sample_selection.py b/spectroscopy/code_src/sample_selection.py index ee14ddea..873003d9 100644 --- a/spectroscopy/code_src/sample_selection.py +++ b/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): diff --git a/spectroscopy/spectra_generator.md b/spectroscopy/spectra_generator.md index dc4cfc2f..b5507e0a 100644 --- a/spectroscopy/spectra_generator.md +++ b/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