diff --git a/spectroscopy/code_src/desi_functions.py b/spectroscopy/code_src/desi_functions.py index 7d8c7b7f..a3abe0e8 100644 --- a/spectroscopy/code_src/desi_functions.py +++ b/spectroscopy/code_src/desi_functions.py @@ -5,6 +5,7 @@ import astropy.units as u from astropy.coordinates import SkyCoord from astropy.table import Table +from astropy import nddata import pandas as pd @@ -47,7 +48,7 @@ def DESIBOSS_get_spec(sample_table, search_radius_arcsec): ## Search data_releases = ['DESI-EDR','BOSS-DR16'] - #data_releases = ['DESI-EDR','BOSS-DR16','SDSS-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) diff --git a/spectroscopy/code_src/mast_functions.py b/spectroscopy/code_src/mast_functions.py index 05782e65..500f303d 100644 --- a/spectroscopy/code_src/mast_functions.py +++ b/spectroscopy/code_src/mast_functions.py @@ -1,8 +1,10 @@ -import os +import os, sys, io import numpy as np +from contextlib import redirect_stdout import astropy.units as u +import astropy.constants as const from astropy.coordinates import SkyCoord from astropy.table import Table @@ -14,7 +16,257 @@ from specutils import Spectrum1D -def HST_get_spec(sample_table, search_radius_arcsec, datadir): +import matplotlib.pyplot as plt + + +def JWST_get_spec(sample_table, search_radius_arcsec, datadir, verbose): + ''' + Retrieves HST spectra for a list of sources and groups/stacks them. + This main function runs two sub-functions: + - JWST_get_spec_helper() which searches, downloads, retrieves the spectra + - JWST_group_spectra() which groups and stacks the spectra + + Parameters + ---------- + sample_table : `~astropy.table.Table` + Table with the coordinates and journal reference labels of the sources + search_radius_arcsec : `float` + Search radius in arcseconds. + datadir : `str` + Data directory where to store the data. Each function will create a + separate data directory (for example "[datadir]/HST/" for HST data). + verbose : `bool` + Verbosity level. Set to True for extra talking. + + Returns + ------- + df_spec : MultiIndexDFObject + The main data structure to store all spectra + + ''' + + ## Get the spectra + print("Searching and Downloading Spectra... ") + df_jwst_all = JWST_get_spec_helper(sample_table, search_radius_arcsec, datadir, verbose) + print("done") + + ## Group + print("Grouping Spectra... ") + df_jwst_group = JWST_group_spectra(df_jwst_all, verbose=verbose, quickplot=False) + print("done") + + return(df_jwst_group) + + +def JWST_get_spec_helper(sample_table, search_radius_arcsec, datadir, verbose): + ''' + Retrieves HST spectra for a list of sources. + + Parameters + ---------- + sample_table : `~astropy.table.Table` + Table with the coordinates and journal reference labels of the sources + search_radius_arcsec : `float` + Search radius in arcseconds. + datadir : `str` + Data directory where to store the data. Each function will create a + separate data directory (for example "[datadir]/HST/" for HST data). + verbose : `bool` + Verbosity level. Set to True for extra talking. + + Returns + ------- + df_spec : MultiIndexDFObject + The main data structure to store all spectra + + ''' + + ## Create directory + if not os.path.exists(datadir): + os.mkdir(datadir) + this_data_dir = os.path.join(datadir , "JWST/") + + + ## Initialize multi-index object: + df_spec = MultiIndexDFObject() + + for stab in sample_table: + + print("Processing source {}".format(stab["label"])) + + ## 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'] + ) + print("Number of search results: {}".format(len(query_results))) + + if len(query_results) > 0: # found some spectra + + + ## Retrieve spectra + 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", + 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: + + ## Download (supress output) + trap = io.StringIO() + with redirect_stdout(trap): + download_results = Observations.download_products(data_products_list_filter, download_dir=this_data_dir) + if verbose: print(trap.getvalue()) + + + ## Create table + # NOTE: `download_results` has NOT the same order as `data_products_list_filter`. We therefore + # have to "manually" get the product file names here and then use those to open the files. + keys = ["filters","obs_collection","instrument_name","calib_level","t_obs_release","proposal_id","obsid","objID","distance"] + tab = Table(names=keys + ["productFilename"] , dtype=[str,str,str,int,float,int,int,int,float]+[str]) + for jj in range(len(data_products_list_filter)): + idx_cross = np.where(query_results["obsid"] == 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)): + + # find correct path name: + # Note that `download_results` does NOT have same order as `tab`!! + file_idx = np.where( [ tab["productFilename"][jj] in download_results["Local Path"][iii] for iii in range(len(download_results))] )[0] + + # open spectrum + # Note that specutils returns the wrong units. Use Table.read() instead. + 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.spectral_axis] , flux=[spec1d.flux], err=[np.repeat(0,len(spec1d.flux))],, + 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]], + )).set_index(["objectid", "label", "filter", "mission"]) + df_spec.append(dfsingle) + + else: + print("Nothing to download for source {}.".format(stab["label"])) + else: + print("Source {} could not be found".format(stab["label"])) + + + return(df_spec) + +def JWST_group_spectra(df, verbose, quickplot): + ''' + Groups the JWST spectra and removes entries that have no spectra. Stacks + spectra that are similar and creates new DF. + + Parameters + ---------- + df : MultiIndexDFObject + Raw JWST multi-index object (output from `JWST_get_spec()`). + verbose : bool + Flag for verbosity: `True` or `False` + quickplot : bool + If `True`, quick plots are made for each spectral group. + + Returns + ------- + df_cons : MultiIndexDFObject + Consolidated and grouped data structure storing the spectra. + + ''' + + ## Initialize multi-index object: + df_spec = MultiIndexDFObject() + + ## Create data table from DF. + tab = df.data.reset_index() + + ## Get objects + objects_unique = np.unique(tab["label"]) + + for obj in objects_unique: + print("Grouping object {}".format(obj)) + + ## Get filters + filters_unique = np.unique(tab["filter"]) + if verbose: print("List of filters in data frame: {}".format( " | ".join(filters_unique)) ) + + for filt in filters_unique: + if verbose: 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=" | ") + + ## get good ones + sum_flux = np.asarray([np.nansum(tab_sel.iloc[iii]["flux"]).value for iii in range(len(tab_sel)) ]) + idx_good = np.where(sum_flux > 0)[0] + if verbose: print("Number of good spectra: {}".format(len(idx_good))) + + if len(idx_good) > 0: + ## Create wavelength grid + wave_grid = tab_sel.iloc[idx_good[0]]["wave"] # NEED TO BE MADE BETTER LATER + + ## Interpolate fluxes + fluxes_int = np.asarray([ np.interp(wave_grid , tab_sel.iloc[idx]["wave"] , tab_sel.iloc[idx]["flux"]) for idx in idx_good ]) + fluxes_units = [tab_sel.iloc[idx]["flux"].unit for idx in idx_good] + fluxes_stack = np.nanmedian(fluxes_int,axis=0) + if verbose: print("Units of fluxes for each spectrum: {}".format( ",".join([str(tt) for tt in fluxes_units]) ) ) + + ## Unit conversion to erg/s/cm2/A (note fluxes are nominally in Jy. So have to do the step with dividing by lam^2) + fluxes_stack_cgs = (fluxes_stack * fluxes_units[0]).to(u.erg / u.second / (u.centimeter**2) / u.hertz) * (const.c.to(u.angstrom/u.second)) / (wave_grid.to(u.angstrom)**2) + fluxes_stack_cgs = fluxes_stack_cgs.to(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"]) + df_spec.append(dfsingle) + + ## Quick plot + if quickplot: + tmp = np.percentile(fluxes_stack, q=(1,50,99) ) + plt.plot(wave_grid , fluxes_stack) + plt.ylim(tmp[0],tmp[2]) + plt.xlabel(r"Observed Wavelength [$\mu$m]") + plt.ylabel(r"Flux [Jy]") + plt.show() + + + return(df_spec) + +def HST_get_spec(sample_table, search_radius_arcsec, datadir, verbose): ''' Retrieves HST spectra for a list of sources. @@ -27,6 +279,8 @@ def HST_get_spec(sample_table, search_radius_arcsec, datadir): datadir : `str` Data directory where to store the data. Each function will create a separate data directory (for example "[datadir]/HST/" for HST data). + verbose : `bool` + Verbosity level. Set to True for extra talking. Returns ------- @@ -36,6 +290,8 @@ def HST_get_spec(sample_table, search_radius_arcsec, datadir): ''' ## Create directory + if not os.path.exists(datadir): + os.mkdir(datadir) this_data_dir = os.path.join(datadir , "HST/") @@ -71,25 +327,36 @@ def HST_get_spec(sample_table, search_radius_arcsec, datadir): if len(data_products_list_filter) > 0: ## Download - download_results = Observations.download_products(data_products_list_filter, download_dir=this_data_dir) + trap = io.StringIO() + with redirect_stdout(trap): + download_results = Observations.download_products(data_products_list_filter, download_dir=this_data_dir) + if verbose: print(trap.getvalue()) - ## Create table + # NOTE: `download_results` has NOT the same order as `data_products_list_filter`. We therefore + # have to "manually" get the product file names here and then use those to open the files. keys = ["filters","obs_collection","instrument_name","calib_level","t_obs_release","proposal_id","obsid","objID","distance"] - tab = Table(names=keys , dtype=[str,str,str,int,float,int,int,int,float]) - for jj in range(len(download_results)): - tmp = query_results[query_results["obsid"] == data_products_list_filter["obsID"][jj]][keys] - tab.add_row( list(tmp[0]) ) + tab = Table(names=keys + ["productFilename"] , dtype=[str,str,str,int,float,int,int,int,float]+[str]) + for jj in range(len(data_products_list_filter)): + idx_cross = np.where(query_results["obsid"] == 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]] ) ## Create multi-index object for jj in range(len(tab)): - # open spectrum - filepath = download_results[jj]["Local Path"] - print(filepath) - spec1d = Spectrum1D.read(filepath) + # find correct path name: + # Note that `download_results` does NOT have same order as `tab`!! + file_idx = np.where( [ tab["productFilename"][jj] in download_results["Local Path"][iii] for iii in range(len(download_results))] )[0] - dfsingle = pd.DataFrame(dict(wave=[spec1d.spectral_axis] , flux=[spec1d.flux], err=[np.repeat(0,len(spec1d.flux))], + # 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. + dfsingle = pd.DataFrame(dict(#wave=[spec1d.spectral_axis] , flux=[spec1d.flux], err=[np.repeat(0,len(spec1d.flux))], + wave=[spec1d.spectral_axis] , flux=[spec1d.flux], err=[spec1d.uncertainty.array * spec1d.uncertainty.unit], label=[stab["label"]], objectid=[stab["objectid"]], #objID=[tab["objID"][jj]], diff --git a/spectroscopy/code_src/plot_functions.py b/spectroscopy/code_src/plot_functions.py index cf71389d..972fa5bb 100644 --- a/spectroscopy/code_src/plot_functions.py +++ b/spectroscopy/code_src/plot_functions.py @@ -3,8 +3,10 @@ import numpy as np import astropy.units as u +import astropy.constants as const from astropy.coordinates import SkyCoord from astropy.table import Table +from astropy.stats import sigma_clip import matplotlib.pyplot as plt import matplotlib as mpl @@ -63,10 +65,26 @@ def bin_spectra(wave,flux, bin_factor): dlam = np.nanmedian(np.diff(wave.value)) * bin_factor + ## This is the faster way, however, it outputs warnings because there are empty slices. + #lam_bin = np.arange(np.nanmin(wave.value)+dlam/2, np.nanmax(wave.value)+dlam/2 + dlam , dlam) + #flux_bin = np.asarray( [ np.nanmedian(flux[(wave.value >= (ll-dlam/2)) & (wave.value < (ll+dlam/2) )].value) for ll in lam_bin ] ) + + ## This way is a bit slower but we can avoid empty slices. lam_bin = np.arange(np.nanmin(wave.value)+dlam/2, np.nanmax(wave.value)+dlam/2 + dlam , dlam) - flux_bin = np.asarray( [ np.nanmedian(flux[(wave.value >= (ll-dlam/2)) & (wave.value < (ll+dlam/2) )].value) for ll in lam_bin ] ) + lam_bins = [] + flux_bins = [] + for ll in lam_bin: + sel_tmp = np.where( (wave.value >= (ll-dlam/2)) & (wave.value < (ll+dlam/2) ) )[0] + if len(sel_tmp) > 0: + flux_bins.append( np.nanmedian(flux[sel_tmp].value) ) + lam_bins.append(ll) + + ## Add the units back! + lam_bins = np.asarray(lam_bins) * wave[0].unit + flux_bins = np.asarray(flux_bins) * flux[0].unit + dlam = dlam * wave[0].unit - return(lam_bin , flux_bin, dlam) + return(lam_bins , flux_bins, dlam) def create_figures(df_spec, bin_factor, show_nbr_figures , save_output): @@ -89,7 +107,6 @@ def create_figures(df_spec, bin_factor, show_nbr_figures , save_output): Whether to save the lightcurve figures. If saved, they will be in the "output" directory. ''' - for cc, (objectid, singleobj_df) in enumerate(df_spec.data.groupby('objectid')): @@ -107,20 +124,32 @@ def create_figures(df_spec, bin_factor, show_nbr_figures , save_output): for ff, (filt,filt_df) in enumerate(singleobj_df.groupby('filter')): - print("{} entries for a object {} and filter {}".format(len(filt_df.flux), objectid , filt)) + #print("{} entries for a object {} and filter {}".format(len(filt_df.flux), objectid , filt)) for ii in range(len(filt_df.flux)): - - wave = filt_df.reset_index().wave[ii] + + # get data + wave = filt_df.reset_index().wave[ii].to(u.micrometer) flux = filt_df.reset_index().flux[ii] err = filt_df.reset_index().err[ii] - mask = np.where(np.isfinite(err))[0] + + # do masking to remove value that are not finite + mask = np.where((np.isfinite(err)) & (flux > 0))[0] wave = wave[mask] flux = flux[mask] err = err[mask] - - #ax1.plot(wave/1e4 , flux , "-" , label="{} ({})".format(filt, filt_df.reset_index().mission[0])) + + #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) - ax1.step(wave_bin/1e4 , flux_bin , "-" , label="{} ({})".format(filt, filt_df.reset_index().mission[ii]), where="mid") + + # do some more clearning (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] + clip_mask = sigma_clip(flux_bin , sigma_lower = 3, cenfunc=np.nanmedian , sigma_upper = 5, stdfunc = np.nanstd).mask + 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") # all in um and.. ax1.set_title(this_label) if LOGX: ax1.set_xscale("log") diff --git a/spectroscopy/code_src/sample_selection.py b/spectroscopy/code_src/sample_selection.py index 1a328e12..ad6a049a 100644 --- a/spectroscopy/code_src/sample_selection.py +++ b/spectroscopy/code_src/sample_selection.py @@ -8,7 +8,7 @@ #from astroquery.vizier import Vizier -def clean_sample(coords_list, labels_list, verbose=1): +def clean_sample(coords_list, labels_list, precision, verbose=1): """Makes a unique sample of skycoords and labels with no repeats. Attaches an object ID to the coords. Parameters @@ -17,6 +17,8 @@ def clean_sample(coords_list, labels_list, verbose=1): list of Astropy SkyCoords derived from literature sources labels_list : list List of the first author name and publication year for tracking the sources + precision : float (astropy units) + Precision of matching/removing doubles. For example 0.5*u.arcsecond. verbose : int, optional Print out the length of the sample after applying this function @@ -31,7 +33,7 @@ def clean_sample(coords_list, labels_list, verbose=1): # now join the table with itself within a defined radius. # We keep one set of original column names to avoid later need for renaming tjoin = join(sample_table, sample_table, keys='coord', - join_funcs={'coord': join_skycoord(0.005 * u.deg)}, + join_funcs={'coord': join_skycoord(precision)}, uniq_col_name='{col_name}{table_name}', table_names=['', '_2']) # this join will return 4 entries for each redundant coordinate: diff --git a/spectroscopy/code_src/sdss_functions.py b/spectroscopy/code_src/sdss_functions.py index c5aea391..ffffe54c 100644 --- a/spectroscopy/code_src/sdss_functions.py +++ b/spectroscopy/code_src/sdss_functions.py @@ -14,7 +14,7 @@ from specutils import Spectrum1D -def SDSS_get_spec(sample_table, search_radius_arcsec): +def SDSS_get_spec(sample_table, search_radius_arcsec, data_release): ''' Retrieves SDSS spectra for a list of sources. Note that no data will be directly downloaded. All will be saved in cache. @@ -25,6 +25,8 @@ def SDSS_get_spec(sample_table, search_radius_arcsec): Table with the coordinates and journal reference labels of the sources search_radius_arcsec : `float` Search radius in arcseconds. + data_release : `int` + SDSS data release (e.g., 17 or 18) Returns ------- @@ -42,10 +44,10 @@ def SDSS_get_spec(sample_table, search_radius_arcsec): ## Get Spectra for SDSS search_coords = stab["coord"] - xid = SDSS.query_region(search_coords, radius=search_radius_arcsec * u.arcsec, spectro=True) + xid = SDSS.query_region(search_coords, radius=search_radius_arcsec * u.arcsec, spectro=True, data_release=data_release) if str(type(xid)) != "": - sp = SDSS.get_spectra(matches=xid, show_progress=True) + sp = SDSS.get_spectra(matches=xid, show_progress=True , data_release=data_release) ## Get data wave = 10**sp[0]["COADD"].data.loglam * u.angstrom # only one entry because we only search for one xid at a time. Could change that? diff --git a/spectroscopy/spectra_generator.md b/spectroscopy/spectra_generator.md index 785e952e..81db5990 100644 --- a/spectroscopy/spectra_generator.md +++ b/spectroscopy/spectra_generator.md @@ -5,11 +5,11 @@ jupyter: extension: .md format_name: markdown format_version: '1.3' - jupytext_version: 1.16.0 + jupytext_version: 1.15.2 kernelspec: - display_name: science_demo + display_name: Python 3 (ipykernel) language: python - name: conda-env-science_demo-py + name: python3 --- @@ -45,8 +45,8 @@ The notebook may focus on the COSMOS field for now, which has a large overlap of | IRSA | Herschel* | Some spectra, need to check reduction stage | | | | IRSA | Euclid | Spectra hosted at IRSA in FY25 -> preparation for ingestion | | Will use mock spectra with correct format for testing | | IRSA | SPHEREx | Spectra/cubes will be hosted at IRSA, first release in FY25 -> preparation for ingestion | | Will use mock spectra with correct format for testing | -| MAST | HST* | Slitless spectra would need reduction and extraction. There are some reduced slit spectra from COS in the Hubble Archive | `astroquery.mast`? | Implemented using `astroquery.mast` | -| MAST | JWST* | Reduced slit MSA spectra that can be queried | `astroquery.mast`? | Should be straight forward using `astroquery.mast` | +| MAST | HST* | Slitless spectra would need reduction and extraction. There are some reduced slit spectra from COS in the Hubble Archive | `astroquery.mast` | Implemented using `astroquery.mast` | +| MAST | JWST* | Reduced slit MSA and Slit spectra that can be queried | `astroquery.mast` | Implemented using `astroquery.mast` | | SDSS | SDSS optical| Optical spectra that are reduced | [Sky Server](https://skyserver.sdss.org/dr18/SearchTools) or `astroquery.sdss` (preferred) | Implemented using `astroquery.sdss`. | | DESI | DESI* | Optical spectra | [DESI public data release](https://data.desi.lbl.gov/public/) | Implemented with `SPARCL` library | | BOSS | BOSS* | Optical spectra | [BOSS webpage (part of SDSS)](https://www.sdss4.org/surveys/boss/) | Implemented with `SPARCL` library together with DESI | @@ -85,7 +85,6 @@ Andreas Faisst, Jessica Krick, Shoubaneh Hemmati, Troy Raen, Brigitta Sipőcz, D - ### Datasets that were considered but didn't end up being used: #### IRTF: - https://irsa.ipac.caltech.edu/Missions/irtf.html \ @@ -93,7 +92,6 @@ Andreas Faisst, Jessica Krick, Shoubaneh Hemmati, Troy Raen, Brigitta Sipőcz, D - large library of stellar spectra \ - Not included here because the data are not currently available in an easily accessible, searchable format - ```python # Ensure all dependencies are installed @@ -103,7 +101,7 @@ Andreas Faisst, Jessica Krick, Shoubaneh Hemmati, Troy Raen, Brigitta Sipőcz, D ```python ## IMPORTS -import sys +import sys, os import numpy as np import matplotlib.pyplot as plt @@ -120,7 +118,7 @@ from sample_selection import clean_sample from desi_functions import DESIBOSS_get_spec from spitzer_functions import SpitzerIRS_get_spec from sdss_functions import SDSS_get_spec -from mast_functions import HST_get_spec +from mast_functions import HST_get_spec, JWST_get_spec from keck_functions import KeckDEIMOS_get_spec from plot_functions import create_figures ``` @@ -152,22 +150,27 @@ coords.append(SkyCoord( 150.1024475 , 2.2815559, unit=u.deg )) labels.append("COSMOS2") coords.append(SkyCoord("{} {}".format("150.000" , "+2.00"), unit=(u.deg, u.deg) )) -labels.append("None") +labels.append("COSMOS3") +coords.append(SkyCoord("{} {}".format("+53.15508" , "-27.80178"), unit=(u.deg, u.deg) )) +labels.append("JADESGS-z7-01-QU") +coords.append(SkyCoord("{} {}".format("+53.15398", "-27.80095"), unit=(u.deg, u.deg) )) +labels.append("TestJWST") -sample_table = clean_sample(coords, labels , verbose=1) +sample_table = clean_sample(coords, labels, precision=2.0* u.arcsecond , verbose=1) -print("Number of sources in sample table: {}".format(len(sample_table))) ``` ### 1.2 Write out your sample to disk -At this point you may wish to write out your sample to disk and reuse that in future work sessions, instead of creating it from scratch again. +At this point you may wish to write out your sample to disk and reuse that in future work sessions, instead of creating it from scratch again. Note that we first check if the `data` directory exists and if not, we will create one. For the format of the save file, we would suggest to choose from various formats that fully support astropy objects(eg., SkyCoord). One example that works is Enhanced Character-Separated Values or ['ecsv'](https://docs.astropy.org/en/stable/io/ascii/ecsv.html) ```python +if not os.path.exists("./data"): + os.mkdir("./data") sample_table.write('data/input_sample.ecsv', format='ascii.ecsv', overwrite = True) ``` @@ -204,33 +207,41 @@ This archive includes spectra taken by ```python %%time - ## Get Keck Spectra (COSMOS only) df_spec_DEIMOS = KeckDEIMOS_get_spec(sample_table = sample_table, search_radius_arcsec=1) df_spec.append(df_spec_DEIMOS) +``` +```python +%%time ## Get Spitzer IRS Spectra df_spec_IRS = SpitzerIRS_get_spec(sample_table, search_radius_arcsec=1 , COMBINESPEC=False) df_spec.append(df_spec_IRS) - ``` ### 2.2 MAST Archive This archive includes spectra taken by - • HST + • HST (including slit spectroscopy) - • JWST (not implemented, yet) + • JWST (including MSA and slit spectroscopy) ```python %%time ## Get Spectra for HST -df_spec_HST = HST_get_spec(sample_table , search_radius_arcsec = 0.5, datadir = "./data/") +df_spec_HST = HST_get_spec(sample_table , search_radius_arcsec = 0.5, datadir = "./data/", verbose = False) df_spec.append(df_spec_HST) ``` +```python +%%time +## Get Spectra for JWST +df_jwst = JWST_get_spec(sample_table , search_radius_arcsec = 0.5, datadir = "./data/", verbose = False) +df_spec.append(df_jwst) +``` + ### 2.3 SDSS Archive This includes SDSS spectra. @@ -238,13 +249,14 @@ This includes SDSS spectra. ```python %%time ## Get SDSS Spectra -df_spec_SDSS = SDSS_get_spec(sample_table , search_radius_arcsec=5) +df_spec_SDSS = SDSS_get_spec(sample_table , search_radius_arcsec=5, data_release=17) df_spec.append(df_spec_SDSS) ``` ### 2.4 DESI Archive -This includes DESI spectra. Here, we use the `SPARCL` query. +This includes DESI spectra. Here, we use the `SPARCL` query. Note that this can also be used +for SDSS searches, however, according to the SPARCL webpage, only up to DR16 is included. Therefore, we will not include SDSS DR16 here (this is treated in the SDSS search above). ```python %%time @@ -253,10 +265,6 @@ df_spec_DESIBOSS = DESIBOSS_get_spec(sample_table, search_radius_arcsec=5) df_spec.append(df_spec_DESIBOSS) ``` -```python -df_spec.data.sort_values(by = ['objectid']) -``` - ## 3. Make plots of luminosity as a function of time We show flux in mJy as a function of time for all available bands for each object. `show_nbr_figures` controls how many plots are actually generated and returned to the screen. If you choose to save the plots with `save_output`, they will be put in the output directory and labelled by sample number. @@ -265,12 +273,12 @@ We show flux in mJy as a function of time for all available bands for each objec ```python ### Plotting #### create_figures(df_spec = df_spec, - bin_factor=10, + bin_factor=5, show_nbr_figures = 10, save_output = False, ) ``` -```python + -``` +