From 7e18abd4f7d4fd287e13846166f721af817aad76 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Sun, 12 Nov 2023 20:39:32 +0100 Subject: [PATCH] Retrieving 2MASS magnitudes from VizieR instead of Simbad, included the Gaia XP spectrum when available --- .gitignore | 2 + Makefile | 2 + calistar/calistar.py | 304 +++++++++++++++++++++++++++++-------------- docs/tutorial.ipynb | 30 ++++- requirements.txt | 1 + 5 files changed, 239 insertions(+), 100 deletions(-) diff --git a/.gitignore b/.gitignore index 074a492..7db57ff 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,8 @@ .DS_Store docs/_build docs/.ipynb_checkpoints/ +docs/6843672087120107264_gaiaxp_0.jpg +docs/6843672087120107264_gaiaxp_spec.dat build/ dist/ calistar.egg-info/ diff --git a/Makefile b/Makefile index 93dc230..7279686 100644 --- a/Makefile +++ b/Makefile @@ -37,6 +37,8 @@ clean: rm -f coverage.xml rm -f docs/calib_* rm -f docs/target_* + rm -f docs/6843672087120107264_gaiaxp_0.jpg + rm -f docs/6843672087120107264_gaiaxp_spec.dat rm -f calib_* rm -f target_* rm -rf .pytest_cache/ diff --git a/calistar/calistar.py b/calistar/calistar.py index 7f2632a..4842184 100644 --- a/calistar/calistar.py +++ b/calistar/calistar.py @@ -21,6 +21,7 @@ from astroquery.gaia import Gaia from astroquery.simbad import Simbad from astroquery.vizier import Vizier +from gaiaxpy import calibrate, plot_spectra from rich import print as rprint from rich.progress import track from typeguard import typechecked @@ -411,18 +412,59 @@ def target_star( if "has_rvs" in gaia_result.columns: print(f"RVS spectrum = {gaia_result['has_rvs'][0]}") - # Add spectral type and 2MASS JHKs magnitudes to the Simbad output + # Gaia XP spectrum + + if "has_xp_continuous" in gaia_result.columns and gaia_result['has_xp_continuous'][0]: + # Sampling adopted from the Gaia XP documentation + df_cal, sampling = calibrate( + input_object=[f"{self.gaia_source}"], + sampling=np.geomspace(330, 1049.9999999999, 361), + truncation=False, + with_correlation=False, + output_path="./", + # output_file=f"{self.gaia_source}_gaiaxp", + output_format=None, + save_file=False, + username=None, + password=None, + ) + + xp_plot = f"{self.gaia_source}_gaiaxp" + print(f"\nStoring Gaia XP plot: {xp_plot}_0.jpg") + + plot_spectra( + spectra=df_cal, + sampling=sampling, + multi=False, + show_plot=False, + output_path="./", + output_file=xp_plot, + format=None, + legend=True, + ) + + xp_wavel = sampling*1e-3 # (nm) -> (um) + xp_flux = 1e3*df_cal["flux"].to_numpy()[0] # (W m-2 nm-1) -> (W m-2 um-1) + xp_error = 1e3*df_cal["flux_error"].to_numpy()[0] + + header = "Wavelength (um) - Flux (W m-2 um-1) - Uncertainty (W m-2 um-1)" + xp_file = f"{self.gaia_source}_gaiaxp_spec.dat" + xp_spec = np.column_stack([xp_wavel, xp_flux, xp_error]) + np.savetxt(xp_file, xp_spec, header=header) + print(f"Storing Gaia XP spectrum: {xp_file}") + + # Add spectral type to the Simbad output # print(Simbad.list_votable_fields()) Simbad.add_votable_fields( "sptype", - "flux(J)", - "flux(H)", - "flux(K)", - "flux_error(J)", - "flux_error(H)", - "flux_error(K)", + # "flux(J)", + # "flux(H)", + # "flux(K)", + # "flux_error(J)", + # "flux_error(H)", + # "flux_error(K)", ) # Simbad query for selected Gaia source ID @@ -439,104 +481,152 @@ def target_star( print(f"Spectral type = {simbad_result['SP_TYPE']}") - print( - f"\n2MASS J mag = {simbad_result['FLUX_J']:.3f} " - f"+/- {simbad_result['FLUX_ERROR_J']:.3f}" - ) - - print( - f"2MASS H mag = {simbad_result['FLUX_H']:.3f} " - f"+/- {simbad_result['FLUX_ERROR_H']:.3f}" - ) - - print( - f"2MASS Ks mag = {simbad_result['FLUX_K']:.3f} " - f"+/- {simbad_result['FLUX_ERROR_K']:.3f}" - ) + # print( + # f"\n2MASS J mag = {simbad_result['FLUX_J']:.3f} " + # f"+/- {simbad_result['FLUX_ERROR_J']:.3f}" + # ) + # + # print( + # f"2MASS H mag = {simbad_result['FLUX_H']:.3f} " + # f"+/- {simbad_result['FLUX_ERROR_H']:.3f}" + # ) + # + # print( + # f"2MASS Ks mag = {simbad_result['FLUX_K']:.3f} " + # f"+/- {simbad_result['FLUX_ERROR_K']:.3f}" + # ) target_dict["Simbad ID"] = simbad_result["MAIN_ID"] target_dict["SpT"] = simbad_result["SP_TYPE"] - target_dict["2MASS/2MASS.J"] = ( - float(simbad_result["FLUX_J"]), - float(simbad_result["FLUX_ERROR_J"]), - ) - - target_dict["2MASS/2MASS.H"] = ( - float(simbad_result["FLUX_H"]), - float(simbad_result["FLUX_ERROR_H"]), - ) - - target_dict["2MASS/2MASS.Ks"] = ( - float(simbad_result["FLUX_K"]), - float(simbad_result["FLUX_ERROR_K"]), - ) + # target_dict["2MASS/2MASS.J"] = ( + # float(simbad_result["FLUX_J"]), + # float(simbad_result["FLUX_ERROR_J"]), + # ) + # + # target_dict["2MASS/2MASS.H"] = ( + # float(simbad_result["FLUX_H"]), + # float(simbad_result["FLUX_ERROR_H"]), + # ) + # + # target_dict["2MASS/2MASS.Ks"] = ( + # float(simbad_result["FLUX_K"]), + # float(simbad_result["FLUX_ERROR_K"]), + # ) # VizieR query for the selected Gaia source ID # Sort the result by distance from the queried object print("\n-> Querying VizieR...\n") - vizier_obj = Vizier(columns=["*", "+_r"], catalog="II/328/allwise") + vizier_obj = Vizier( + columns=["*", "+_r"], catalog=["II/246/out", "II/328/allwise"] + ) radius = u.Quantity(1.0 * u.arcmin) vizier_result = vizier_obj.query_object( f"GAIA {self.gaia_release} {self.gaia_source}", radius=radius ) - vizier_result = vizier_result["II/328/allwise"] - # print(f"Found {len(vizier_result)} object(s) in VizieR. Selecting the first object...") - vizier_result = vizier_result[0] + # 2MASS data from VizieR + + vizier_2mass = vizier_result["II/246/out"] + vizier_2mass = vizier_2mass[0] - print(f"ALLWISE source ID = {vizier_result['AllWISE']}") + print(f"2MASS source ID = {vizier_2mass['_2MASS']}") + + print( + f"Separation between Gaia and 2MASS source = " + f"{1e3*vizier_2mass['_r']:.1f} mas" + ) + + print( + f"\n2MASS J mag = {vizier_2mass['Jmag']:.3f} " + f"+/- {vizier_2mass['e_Jmag']:.3f}" + ) + + print( + f"2MASS H mag = {vizier_2mass['Hmag']:.3f} " + f"+/- {vizier_2mass['e_Hmag']:.3f}" + ) + + print( + f"2MASS Ks mag = {vizier_2mass['Kmag']:.3f} " + f"+/- {vizier_2mass['e_Kmag']:.3f}" + ) + + target_dict["2MASS/2MASS.J"] = ( + float(vizier_2mass["Jmag"]), + float(vizier_2mass["e_Jmag"]), + ) + target_dict["2MASS/2MASS.H"] = ( + float(vizier_2mass["Hmag"]), + float(vizier_2mass["e_Hmag"]), + ) + target_dict["2MASS/2MASS.Ks"] = ( + float(vizier_2mass["Kmag"]), + float(vizier_2mass["e_Kmag"]), + ) + + # WISE data from VizieR + + vizier_wise = vizier_result["II/328/allwise"] + vizier_wise = vizier_wise[0] + + print(f"\nALLWISE source ID = {vizier_wise['AllWISE']}") print( f"Separation between Gaia and ALLWISE source = " - f"{1e3*vizier_result['_r']:.1f} mas" + f"{1e3*vizier_wise['_r']:.1f} mas" ) print( - f"\nWISE W1 mag = {vizier_result['W1mag']:.3f} " - f"+/- {vizier_result['e_W1mag']:.3f}" + f"\nWISE W1 mag = {vizier_wise['W1mag']:.3f} " + f"+/- {vizier_wise['e_W1mag']:.3f}" ) print( - f"WISE W2 mag = {vizier_result['W2mag']:.3f} " - f"+/- {vizier_result['e_W2mag']:.3f}" + f"WISE W2 mag = {vizier_wise['W2mag']:.3f} " + f"+/- {vizier_wise['e_W2mag']:.3f}" ) print( - f"WISE W3 mag = {vizier_result['W3mag']:.3f} " - f"+/- {vizier_result['e_W3mag']:.3f}" + f"WISE W3 mag = {vizier_wise['W3mag']:.3f} " + f"+/- {vizier_wise['e_W3mag']:.3f}" ) print( - f"WISE W4 mag = {vizier_result['W4mag']:.3f} " - f"+/- {vizier_result['e_W4mag']:.3f}" + f"WISE W4 mag = {vizier_wise['W4mag']:.3f} " + f"+/- {vizier_wise['e_W4mag']:.3f}" ) target_dict["WISE/WISE.W1"] = ( - float(vizier_result["W1mag"]), - float(vizier_result["e_W4mag"]), + float(vizier_wise["W1mag"]), + float(vizier_wise["e_W4mag"]), ) + target_dict["WISE/WISE.W2"] = ( - float(vizier_result["W2mag"]), - float(vizier_result["e_W4mag"]), + float(vizier_wise["W2mag"]), + float(vizier_wise["e_W4mag"]), ) + target_dict["WISE/WISE.W3"] = ( - float(vizier_result["W3mag"]), - float(vizier_result["e_W4mag"]), + float(vizier_wise["W3mag"]), + float(vizier_wise["e_W4mag"]), ) + target_dict["WISE/WISE.W4"] = ( - float(vizier_result["W4mag"]), - float(vizier_result["e_W4mag"]), + float(vizier_wise["W4mag"]), + float(vizier_wise["e_W4mag"]), ) if write_json: - json_file = f"target_{self.gaia_release.lower()}_{self.gaia_source}.json" - print(f"\nStoring output: {json_file}") + json_file = f"target_{self.gaia_release.lower()}" \ + f"_{self.gaia_source}.json" + + print(f"\nStoring JSON output: {json_file}") with open(json_file, "w", encoding="utf-8") as open_file: json.dump(target_dict, open_file, indent=4) @@ -612,12 +702,12 @@ def find_calib( Simbad.add_votable_fields( "sptype", "ids", - "flux(J)", - "flux(H)", - "flux(K)", - "flux_error(J)", - "flux_error(H)", - "flux_error(K)", + # "flux(J)", + # "flux(H)", + # "flux(K)", + # "flux_error(J)", + # "flux_error(H)", + # "flux_error(K)", ) print(f"Radius of search cone = {search_radius} deg") @@ -721,61 +811,85 @@ def find_calib( else: cal_df.loc[gaia_item.index, "SpT"] = simbad_result["SP_TYPE"] - if np.ma.is_masked(simbad_result["FLUX_J"]): + # if np.ma.is_masked(simbad_result["FLUX_J"]): + # cal_df.loc[gaia_item.index, "2MASS/2MASS.J"] = np.nan + # else: + # cal_df.loc[gaia_item.index, "2MASS/2MASS.J"] = simbad_result[ + # "FLUX_J" + # ] + # + # if np.ma.is_masked(simbad_result["FLUX_H"]): + # cal_df.loc[gaia_item.index, "2MASS/2MASS.H"] = np.nan + # else: + # cal_df.loc[gaia_item.index, "2MASS/2MASS.H"] = simbad_result[ + # "FLUX_H" + # ] + # + # if np.ma.is_masked(simbad_result["FLUX_K"]): + # cal_df.loc[gaia_item.index, "2MASS/2MASS.Ks"] = np.nan + # else: + # cal_df.loc[gaia_item.index, "2MASS/2MASS.Ks"] = simbad_result[ + # "FLUX_K" + # ] + + radius = u.Quantity(1.0 * u.arcmin) + + vizier_result = vizier_obj.query_object( + f"GAIA {self.gaia_release} {gaia_item['source_id']}", + radius=radius, + catalog=["II/246/out", "II/328/allwise"], + ) + + if len(vizier_result) == 2: + # 2MASS + vizier_2mass = vizier_result["II/246/out"][0] + + # Check if the separation between the Gaia and + # the 2MASS coordinates is at most 200 mas + skip_source = vizier_2mass["_r"] > 0.2 + + if skip_source or np.ma.is_masked(vizier_2mass["Jmag"]): cal_df.loc[gaia_item.index, "2MASS/2MASS.J"] = np.nan else: - cal_df.loc[gaia_item.index, "2MASS/2MASS.J"] = simbad_result[ - "FLUX_J" - ] + cal_df.loc[gaia_item.index, "2MASS/2MASS.J"] = vizier_2mass["Jmag"] - if np.ma.is_masked(simbad_result["FLUX_H"]): + if skip_source or np.ma.is_masked(vizier_2mass["Hmag"]): cal_df.loc[gaia_item.index, "2MASS/2MASS.H"] = np.nan else: - cal_df.loc[gaia_item.index, "2MASS/2MASS.H"] = simbad_result[ - "FLUX_H" - ] + cal_df.loc[gaia_item.index, "2MASS/2MASS.H"] = vizier_2mass["Hmag"] - if np.ma.is_masked(simbad_result["FLUX_K"]): + if skip_source or np.ma.is_masked(vizier_2mass["Kmag"]): cal_df.loc[gaia_item.index, "2MASS/2MASS.Ks"] = np.nan else: - cal_df.loc[gaia_item.index, "2MASS/2MASS.Ks"] = simbad_result[ - "FLUX_K" - ] - - radius = u.Quantity(1.0 * u.arcmin) + cal_df.loc[gaia_item.index, "2MASS/2MASS.Ks"] = vizier_2mass["Kmag"] - vizier_result = vizier_obj.query_object( - f"GAIA {self.gaia_release} {gaia_item['source_id']}", - radius=radius, - catalog="II/328/allwise", - ) + # WISE - if len(vizier_result) == 1: - vizier_result = vizier_result["II/328/allwise"][0] + vizier_wise = vizier_result["II/328/allwise"][0] # Check if the separation between the Gaia and # the ALLWISE coordinates is at most 200 mas - skip_source = vizier_result["_r"] > 0.2 + skip_source = vizier_wise["_r"] > 0.2 - if skip_source or np.ma.is_masked(vizier_result["W1mag"]): + if skip_source or np.ma.is_masked(vizier_wise["W1mag"]): cal_df.loc[gaia_item.index, "WISE/WISE.W1"] = np.nan else: - cal_df.loc[gaia_item.index, "WISE/WISE.W1"] = vizier_result["W1mag"] + cal_df.loc[gaia_item.index, "WISE/WISE.W1"] = vizier_wise["W1mag"] - if skip_source or np.ma.is_masked(vizier_result["W2mag"]): + if skip_source or np.ma.is_masked(vizier_wise["W2mag"]): cal_df.loc[gaia_item.index, "WISE/WISE.W2"] = np.nan else: - cal_df.loc[gaia_item.index, "WISE/WISE.W2"] = vizier_result["W2mag"] + cal_df.loc[gaia_item.index, "WISE/WISE.W2"] = vizier_wise["W2mag"] - if skip_source or np.ma.is_masked(vizier_result["W3mag"]): + if skip_source or np.ma.is_masked(vizier_wise["W3mag"]): cal_df.loc[gaia_item.index, "WISE/WISE.W3"] = np.nan else: - cal_df.loc[gaia_item.index, "WISE/WISE.W3"] = vizier_result["W3mag"] + cal_df.loc[gaia_item.index, "WISE/WISE.W3"] = vizier_wise["W3mag"] - if skip_source or np.ma.is_masked(vizier_result["W4mag"]): + if skip_source or np.ma.is_masked(vizier_wise["W4mag"]): cal_df.loc[gaia_item.index, "WISE/WISE.W4"] = np.nan else: - cal_df.loc[gaia_item.index, "WISE/WISE.W4"] = vizier_result["W4mag"] + cal_df.loc[gaia_item.index, "WISE/WISE.W4"] = vizier_wise["W4mag"] # This query returns no sources? # gaia_query = f""" diff --git a/docs/tutorial.ipynb b/docs/tutorial.ipynb index 408f283..ac827b2 100644 --- a/docs/tutorial.ipynb +++ b/docs/tutorial.ipynb @@ -117,18 +117,38 @@ "XP continuous = True\n", "XP sampled = True\n", "RVS spectrum = False\n", + "Reading input file...me...\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + " 0/1 [00:00 Querying Simbad...\n", "\n", "Simbad ID = HD 206893\n", "Spectral type = F5V\n", "\n", + "-> Querying VizieR...\n", + "\n", + "2MASS source ID = 21452190-1246599\n", + "Separation between Gaia and 2MASS source = 1.4 mas\n", + "\n", "2MASS J mag = 5.869 +/- 0.023\n", "2MASS H mag = 5.687 +/- 0.034\n", "2MASS Ks mag = 5.593 +/- 0.021\n", "\n", - "-> Querying VizieR...\n", - "\n", "ALLWISE source ID = J214521.98-124659.9\n", "Separation between Gaia and ALLWISE source = 18.9 mas\n", "\n", @@ -137,7 +157,7 @@ "WISE W3 mag = 5.629 +/- 0.015\n", "WISE W4 mag = 5.481 +/- 0.043\n", "\n", - "Storing output: target_dr3_6843672087120107264.json\n" + "Storing JSON output: target_dr3_6843672087120107264.json\n" ] } ], @@ -183,7 +203,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "84d1d98fbcef458e98fa611e1243a19c", + "model_id": "af3fceb7cea7429b8b3e88f043b05add", "version_major": 2, "version_minor": 0 }, @@ -299,7 +319,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "b2f457739610474394b62a4b3edb0b05", + "model_id": "c1fe7e9fede142b18e1206f65e90980e", "version_major": 2, "version_minor": 0 }, diff --git a/requirements.txt b/requirements.txt index 0789336..97a18c0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,6 @@ astropy astroquery +gaiaxp matplotlib numpy pandas