Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Gaia targets #154

Open
wants to merge 12 commits into
base: v3.1.0
Choose a base branch
from
Open

Gaia targets #154

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
f7d0701
first version, with gaia targets loaded inside getCalspec.
corentinravoux Jun 19, 2024
e112d2f
adding new gaiaspec package
corentinravoux Jul 17, 2024
3166760
importation of gaia module with try
corentinravoux Jul 17, 2024
4f31e6d
getting ra dec and z for gaia
corentinravoux Jul 17, 2024
c959105
small fix
corentinravoux Jul 17, 2024
4eb59e4
adding gaiaspec in fit spectrum and spectrogram
corentinravoux Jul 23, 2024
1dd0210
making gaiaspec loading more explicit
corentinravoux Aug 21, 2024
0a15587
adapt to gaiaspec changes
corentinravoux Aug 21, 2024
ac7e604
Fix to not download everytime the gaia cache spectrum + clearer code
corentinravoux Aug 23, 2024
5aa42b2
Merge branch 'v3.1.0' of github.com:LSSTDESC/Spectractor into gaia_ta…
jeremyneveu Aug 23, 2024
9958b70
missing else
jeremyneveu Aug 23, 2024
410f5a3
add sed error
jeremyneveu Sep 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
271 changes: 182 additions & 89 deletions spectractor/extractor/targets.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,12 @@

from getCalspec import getCalspec

try:
from gaiaspec import getGaia
except ModuleNotFoundError:
getGaia = None



def load_target(label, verbose=False):
"""Load the target properties according to the type set by parameters.OBS_OBJECT_TYPE.
Expand Down Expand Up @@ -75,10 +81,12 @@ def __init__(self, label, verbose=False):
self.type = None
self.wavelengths = []
self.spectra = []
self.spectra_err = []
self.verbose = verbose
self.emission_spectrum = False
self.hydrogen_only = False
self.sed = None
self.sed_err = None
self.lines = None
self.radec_position = None
self.radec_position_after_pm = None
Expand Down Expand Up @@ -231,7 +239,6 @@ def __init__(self, label, verbose=False):
"""
Target.__init__(self, label, verbose=verbose)
self.my_logger = set_logger(self.__class__.__name__)
self.simbad_table = None
self.load()

def load(self):
Expand All @@ -257,30 +264,53 @@ def load(self):
# ``Simbad...`` methods secretly makes an instance, which stays around,
# has a connection go stale, and then raises an exception seemingly
# at some random time later
simbadQuerier = SimbadClass()
patchSimbadURL(simbadQuerier)

simbadQuerier.add_votable_fields('flux(U)', 'flux(B)', 'flux(V)', 'flux(R)', 'flux(I)', 'flux(J)', 'sptype',
'parallax', 'pm', 'z_value')
if not getCalspec.is_calspec(self.label) and getCalspec.is_calspec(self.label.replace(".", " ")):
self.label = self.label.replace(".", " ")
astroquery_label = self.label
if getCalspec.is_calspec(self.label):
calspec = getCalspec.Calspec(self.label)
astroquery_label = calspec.Astroquery_Name
self.simbad_table = simbadQuerier.query_object(astroquery_label)

if self.simbad_table is not None:
if self.verbose or True:
self.my_logger.info(f'\n\tSimbad:\n{self.simbad_table}')
self.radec_position = SkyCoord(self.simbad_table['RA'][0] + ' ' + self.simbad_table['DEC'][0], unit=(u.hourangle, u.deg))
else:
raise RuntimeError(f"Target {self.label} not found in Simbad")
self.get_radec_position_after_pm(date_obs="J2000")
if not np.ma.is_masked(self.simbad_table['Z_VALUE']):
self.redshift = float(self.simbad_table['Z_VALUE'])

if getGaia is None:
is_gaiaspec = False
else:
is_gaiaspec = getGaia.is_gaiaspec(self.label)
if is_gaiaspec:
# In the case where the gaia spectrum is not available in the
# package data, get source info from simbadQuerier
gaia_sources = getGaia.get_gaia_sources()
source = gaia_sources[gaia_sources == self.label]
table_coordinates = [{"PMRA": source["pmra"].iloc[0],
"PMDEC": source["pmdec"].iloc[0],
"PLX_VALUE": source["parallax"].iloc[0],}]

self.radec_position = SkyCoord(ra = source["ra"].iloc[0],
dec = source["dec"].iloc[0] ,
unit=u.deg)
self.redshift = 0
date_reference="J2016"
else:
simbadQuerier = SimbadClass()
patchSimbadURL(simbadQuerier)
simbadQuerier.add_votable_fields('flux(U)', 'flux(B)', 'flux(V)', 'flux(R)', 'flux(I)', 'flux(J)', 'sptype',
'parallax', 'pm', 'z_value')
if not getCalspec.is_calspec(self.label) and getCalspec.is_calspec(self.label.replace(".", " ")):
self.label = self.label.replace(".", " ")
astroquery_label = self.label
if getCalspec.is_calspec(self.label):
calspec = getCalspec.Calspec(self.label)
astroquery_label = calspec.Astroquery_Name
table_coordinates = simbadQuerier.query_object(astroquery_label)

if table_coordinates is not None:
if self.verbose or True:
self.my_logger.info(f'\n\tSimbad:\n{table_coordinates}')
self.radec_position = SkyCoord(table_coordinates['RA'][0] + ' ' + table_coordinates['DEC'][0], unit=(u.hourangle, u.deg))
else:
raise RuntimeError(f"Target {self.label} not found in Simbad")
if not np.ma.is_masked(table_coordinates['Z_VALUE']):
self.redshift = float(table_coordinates['Z_VALUE'])
else:
self.redshift = 0
date_reference="J2000"

self.get_radec_position_after_pm(table_coordinates,
date_obs="J2000",
date_reference = date_reference)
self.load_spectra()

def load_spectra(self):
Expand All @@ -301,74 +331,29 @@ def load_spectra(self):
"""
self.wavelengths = [] # in nm
self.spectra = []
self.spectra_err = []
# first try if it is a Calspec star
is_calspec = getCalspec.is_calspec(self.label)
if getGaia is None:
is_gaiaspec = False
is_gaia_full = False
else:
is_gaiaspec = getGaia.is_gaiaspec(self.label)
is_gaia_full = False
if is_gaiaspec == False:
is_gaia_full = getGaia.is_gaia_full(self.label)
if is_calspec:
calspec = getCalspec.Calspec(self.label)
self.emission_spectrum = False
self.hydrogen_only = False
self.lines = Lines(HYDROGEN_LINES + ATMOSPHERIC_LINES + STELLAR_LINES,
redshift=self.redshift, emission_spectrum=self.emission_spectrum,
hydrogen_only=self.hydrogen_only)
spec_dict = calspec.get_spectrum_numpy()
# official units in spectractor are nanometers for wavelengths and erg/s/cm2/nm for fluxes
spec_dict["WAVELENGTH"] = spec_dict["WAVELENGTH"].to(u.nm)
spec_dict["FLUX"] = spec_dict["FLUX"].to(u.erg / u.second / u.cm**2 / u.nm)
self.wavelengths.append(spec_dict["WAVELENGTH"].value)
self.spectra.append(spec_dict["FLUX"].value)
self.load_calspec()
elif is_gaiaspec|is_gaia_full:
self.load_gaia()
# TODO DM-33731: the use of self.label in parameters.STAR_NAMES:
# below works for running but breaks a test so needs fixing for DM
elif 'HD' in self.label: # or self.label in parameters.STAR_NAMES: # it is a star
self.emission_spectrum = False
self.hydrogen_only = False
self.lines = Lines(ATMOSPHERIC_LINES + HYDROGEN_LINES + STELLAR_LINES,
redshift=self.redshift, emission_spectrum=self.emission_spectrum,
hydrogen_only=self.hydrogen_only)
self.load_emission_spectrum(hydrogen_only_flag=False)
elif 'PNG' in self.label:
self.emission_spectrum = True
self.lines = Lines(ATMOSPHERIC_LINES + ISM_LINES + HYDROGEN_LINES,
redshift=self.redshift, emission_spectrum=self.emission_spectrum,
hydrogen_only=self.hydrogen_only)
self.load_emission_spectrum(hydrogen_only_flag=True)
else: # maybe a quasar, try with NED query
from astroquery.ned import Ned
try:
hdulists = Ned.get_spectra(self.label) #, show_progress=False)
except Exception as err:
raise err
if len(hdulists) > 0:
self.emission_spectrum = True
self.hydrogen_only = False
if self.redshift > 0.2:
self.hydrogen_only = True
parameters.LAMBDA_MIN *= 1 + self.redshift
parameters.LAMBDA_MAX *= 1 + self.redshift
self.lines = Lines(ATMOSPHERIC_LINES+ISM_LINES+HYDROGEN_LINES,
redshift=self.redshift, emission_spectrum=self.emission_spectrum,
hydrogen_only=self.hydrogen_only)
for k, h in enumerate(hdulists):
if h[0].header['NAXIS'] == 1:
self.spectra.append(h[0].data)
else:
for d in h[0].data:
self.spectra.append(d)
wave_n = len(h[0].data)
if h[0].header['NAXIS'] == 2:
wave_n = len(h[0].data.T)
wave_step = h[0].header['CDELT1']
wave_start = h[0].header['CRVAL1'] - (h[0].header['CRPIX1'] - 1) * wave_step
wave_end = wave_start + wave_n * wave_step
waves = np.linspace(wave_start, wave_end, wave_n)
is_angstrom = False
for key in list(h[0].header.keys()):
if 'angstrom' in str(h[0].header[key]).lower():
is_angstrom = True
if is_angstrom:
waves *= 0.1
if h[0].header['NAXIS'] > 1:
for i in range(h[0].header['NAXIS'] + 1):
self.wavelengths.append(waves)
else:
self.wavelengths.append(waves)
self.load_ned()
self.build_sed()
self.my_logger.debug(f"\n\tTarget label: {self.label}"
f"\n\tCalspec? {is_calspec}"
Expand All @@ -378,21 +363,125 @@ def load_spectra(self):
if self.lines is not None and len(self.lines.lines) > 0:
self.my_logger.debug(f"\n\tLines: {[l.label for l in self.lines.lines]}")

def get_radec_position_after_pm(self, date_obs):
if self.simbad_table is not None:
target_pmra = self.simbad_table[0]['PMRA'] * u.mas / u.yr
def load_calspec(self):
calspec = getCalspec.Calspec(self.label)
self.emission_spectrum = False
self.hydrogen_only = False
self.lines = Lines(
HYDROGEN_LINES + ATMOSPHERIC_LINES + STELLAR_LINES,
redshift=self.redshift,
emission_spectrum=self.emission_spectrum,
hydrogen_only=self.hydrogen_only,
)
spec_dict = calspec.get_spectrum_numpy()
# official units in spectractor are nanometers for wavelengths and erg/s/cm2/nm for fluxes
spec_dict["WAVELENGTH"] = spec_dict["WAVELENGTH"].to(u.nm)
for key in ["FLUX", "STATERROR", "SYSERROR"]:
spec_dict[key] = spec_dict[key].to(u.erg / u.second / u.cm**2 / u.nm)
self.wavelengths.append(spec_dict["WAVELENGTH"].value)
self.spectra.append(spec_dict["FLUX"].value)
self.spectra_err.append(np.sqrt(spec_dict["STATERROR"].value**2+spec_dict["SYSERROR"].value**2))

def load_gaia(self):
gaia = getGaia.Gaia(self.label)
self.emission_spectrum = False
self.hydrogen_only = False
self.lines = Lines(
HYDROGEN_LINES + ATMOSPHERIC_LINES + STELLAR_LINES,
redshift=self.redshift,
emission_spectrum=self.emission_spectrum,
hydrogen_only=self.hydrogen_only,
)
spec_dict = gaia.get_spectrum_numpy()
# official units in spectractor are nanometers for wavelengths and erg/s/cm2/nm for fluxes
spec_dict["WAVELENGTH"] = spec_dict["WAVELENGTH"].to(u.nm)
for key in ["FLUX", "STATERROR", "SYSERROR"]:
spec_dict[key] = spec_dict[key].to(u.erg / u.second / u.cm**2 / u.nm)
self.wavelengths.append(spec_dict["WAVELENGTH"].value)
self.spectra.append(spec_dict["FLUX"].value)
self.spectra_err.append(np.sqrt(spec_dict["STATERROR"].value**2+spec_dict["SYSERROR"].value**2))

def load_emission_spectrum(self, hydrogen_only_flag):
if hydrogen_only_flag:
self.emission_spectrum = True
self.lines = Lines(
ATMOSPHERIC_LINES + ISM_LINES + HYDROGEN_LINES,
redshift=self.redshift,
emission_spectrum=self.emission_spectrum,
hydrogen_only=self.hydrogen_only,
)
else:
self.emission_spectrum = False
self.hydrogen_only = False
self.lines = Lines(
ATMOSPHERIC_LINES + HYDROGEN_LINES + STELLAR_LINES,
redshift=self.redshift,
emission_spectrum=self.emission_spectrum,
hydrogen_only=self.hydrogen_only,
)

def load_ned(self):
from astroquery.ned import Ned
try:
hdulists = Ned.get_spectra(self.label) # , show_progress=False)
except Exception as err:
raise err
if len(hdulists) > 0:
self.emission_spectrum = True
self.hydrogen_only = False
if self.redshift > 0.2:
self.hydrogen_only = True
parameters.LAMBDA_MIN *= 1 + self.redshift
parameters.LAMBDA_MAX *= 1 + self.redshift
self.lines = Lines(
ATMOSPHERIC_LINES + ISM_LINES + HYDROGEN_LINES,
redshift=self.redshift,
emission_spectrum=self.emission_spectrum,
hydrogen_only=self.hydrogen_only,
)
for k, h in enumerate(hdulists):
if h[0].header["NAXIS"] == 1:
self.spectra.append(h[0].data)
else:
for d in h[0].data:
self.spectra.append(d)
wave_n = len(h[0].data)
if h[0].header["NAXIS"] == 2:
wave_n = len(h[0].data.T)
wave_step = h[0].header["CDELT1"]
wave_start = (
h[0].header["CRVAL1"] - (h[0].header["CRPIX1"] - 1) * wave_step
)
wave_end = wave_start + wave_n * wave_step
waves = np.linspace(wave_start, wave_end, wave_n)
is_angstrom = False
for key in list(h[0].header.keys()):
if "angstrom" in str(h[0].header[key]).lower():
is_angstrom = True
if is_angstrom:
waves *= 0.1
if h[0].header["NAXIS"] > 1:
for i in range(h[0].header["NAXIS"] + 1):
self.wavelengths.append(waves)
else:
self.wavelengths.append(waves)


def get_radec_position_after_pm(self, table_coordinates, date_obs="J2000", date_reference="J2000"):
if table_coordinates is not None:
target_pmra = table_coordinates[0]['PMRA'] * u.mas / u.yr
if np.isnan(target_pmra):
target_pmra = 0 * u.mas / u.yr
target_pmdec = self.simbad_table[0]['PMDEC'] * u.mas / u.yr
target_pmdec = table_coordinates[0]['PMDEC'] * u.mas / u.yr
if np.isnan(target_pmdec):
target_pmdec = 0 * u.mas / u.yr
target_parallax = self.simbad_table[0]['PLX_VALUE'] * u.mas
target_parallax = table_coordinates[0]['PLX_VALUE'] * u.mas
if target_parallax == 0 * u.mas:
target_parallax = 1e-4 * u.mas
target_coord = SkyCoord(ra=self.radec_position.ra, dec=self.radec_position.dec,
distance=Distance(parallax=target_parallax),
pm_ra_cosdec=target_pmra, pm_dec=target_pmdec, frame='icrs', equinox="J2000",
obstime="J2000")
obstime=date_reference)
self.radec_position_after_pm = target_coord.apply_space_motion(new_obstime=Time(date_obs))
return self.radec_position_after_pm
else:
Expand All @@ -418,9 +507,13 @@ def build_sed(self, index=0):
if len(self.spectra) == 0:
self.sed = interp1d(parameters.LAMBDAS, np.zeros_like(parameters.LAMBDAS), kind='linear', bounds_error=False,
fill_value=0.)
self.sed_err = interp1d(parameters.LAMBDAS, np.zeros_like(parameters.LAMBDAS), kind='linear', bounds_error=False,
fill_value=0.)
else:
self.sed = interp1d(self.wavelengths[index], self.spectra[index], kind='linear', bounds_error=False,
fill_value=0.)
self.sed_err = interp1d(self.wavelengths[index], self.spectra_err[index], kind='linear', bounds_error=False,
fill_value=10*np.max(self.spectra_err[index]))

def plot_spectra(self):
""" Plot the spectra stored in the self.spectra list.
Expand Down
17 changes: 15 additions & 2 deletions spectractor/fit/fit_spectrogram.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,10 @@
from spectractor.simulation.atmosphere import Atmosphere, AtmosphereGrid
from spectractor.fit.fitter import (FitWorkspace, FitParameters, run_minimisation, run_minimisation_sigma_clipping,
write_fitparameter_json)
try:
from gaiaspec import getGaia
except ModuleNotFoundError:
getGaia = None

plot_counter = 0

Expand Down Expand Up @@ -57,7 +61,16 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,

"""
if not getCalspec.is_calspec(spectrum.target.label):
raise ValueError(f"{spectrum.target.label=} must be a CALSPEC star according to getCalspec package.")
if getGaia is not None:
is_gaiaspec = getGaia.is_gaiaspec(spectrum.target.label)
is_gaia_full = False
if is_gaiaspec == False:
is_gaia_full = getGaia.is_gaia_full(spectrum.target.label)
if not is_gaiaspec:
if not is_gaia_full:
raise ValueError(f"{spectrum.target.label=} must be a CALSPEC or GAIA star.")
else:
raise ValueError(f"{spectrum.target.label=} must be a CALSPEC star according to getCalspec package.")
self.spectrum = spectrum
self.filename = spectrum.filename.replace("spectrum", "spectrogram")
self.diffraction_orders = np.arange(spectrum.order, spectrum.order + 3 * np.sign(spectrum.order), np.sign(spectrum.order))
Expand Down Expand Up @@ -119,7 +132,7 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
self.atm_params_indices = np.array([params.get_index(label) for label in ["VAOD", "angstrom_exp", "ozone [db]", "PWV [mm]"]])
# A2 is free only if spectrogram is a simulation or if the order 2/1 ratio is not known and flat
if "A2" in params.labels:
params.fixed[params.get_index(f"A{self.diffraction_orders[1]}")] = False #"A2_T" not in self.spectrum.header
params.fixed[params.get_index(f"A{self.diffraction_orders[1]}")] = False #not getCalspec.is_calspec(spectrum.target.label) #"A2_T" not in self.spectrum.header
if "A3" in params.labels:
params.fixed[params.get_index(f"A{self.diffraction_orders[2]}")] = "A3_T" not in self.spectrum.header
params.fixed[params.get_index(r"shift_x [pix]")] = False # Delta x
Expand Down
Loading
Loading