From ae2108c84c1f9604aa2d03af9fcf8aa489d55919 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Thu, 13 Feb 2025 07:49:41 +0100 Subject: [PATCH 1/2] Support for fitting stationary background objects, simplified a bit the Gaia query, compatibility with latest astruquery, support for Python 3.12 --- .github/workflows/main.yml | 2 +- backtracks/backtracks.py | 104 +++++++++++++++++++++++++++++-------- backtracks/plotting.py | 16 +++++- backtracks/utils.py | 8 ++- pyproject.toml | 2 +- requirements.txt | 2 +- 6 files changed, 106 insertions(+), 28 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 671f59b..6de69e8 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -9,7 +9,7 @@ jobs: strategy: matrix: - python-version: ['3.9', '3.10', '3.11'] + python-version: ['3.9', '3.10', '3.11', '3.12'] steps: - uses: actions/checkout@v2 diff --git a/backtracks/backtracks.py b/backtracks/backtracks.py index fd5be3a..6a67b4a 100644 --- a/backtracks/backtracks.py +++ b/backtracks/backtracks.py @@ -164,7 +164,11 @@ def __init__(self, target_name: str, candidate_file: str, nearby_window: float = target_gaia = Table(hdu_list[1].data, masked=True) self.nearby = Table(hdu_list[2].data, masked=True) - self.gaia_id = target_gaia['SOURCE_ID'][0] + if 'source_id' in target_gaia.columns: + self.gaia_id = target_gaia['source_id'][0] + else: + self.gaia_id = target_gaia['SOURCE_ID'][0] + self.gaia_epoch = target_gaia['ref_epoch'][0] for col in target_gaia.columns.values(): @@ -332,14 +336,15 @@ def query_astrometry(self, nearby_window: float = 0.5): """ # resolve target in simbad - target_result_table = Simbad.query_object(self.target_name) + Simbad.add_votable_fields("ids") + simbad_query = Simbad.query_object(self.target_name) print(f'[BACKTRACKS INFO]: Resolved the target star \'{self.target_name}\' in Simbad!') - # target_result_table.pprint() - # get gaia ID from simbad + + # get Gaia source ID gaia_id = None - for target_id in Simbad.query_objectids(self.target_name)['ID']: - if f'Gaia {self.gaia_release}' in target_id: - gaia_id = int(target_id.replace(f'Gaia {self.gaia_release}', '')) + for id_name in simbad_query['ids'][0].split('|'): + if f'Gaia {self.gaia_release}' in id_name: + gaia_id = int(id_name.replace(f'Gaia {self.gaia_release}', '')) print('[BACKTRACKS INFO]: Resolved target\'s Gaia ID ' f'from Simbad, Gaia {self.gaia_release} {gaia_id}') @@ -348,18 +353,63 @@ def query_astrometry(self, nearby_window: float = 0.5): f"is not found in the selected catalog " f"({Gaia.MAIN_GAIA_TABLE}).") - coord = SkyCoord(ra=target_result_table['RA'][0], - dec=target_result_table['DEC'][0], - unit=(u.hourangle, u.degree), frame='icrs') - width = u.Quantity(50, u.arcsec) - height = u.Quantity(50, u.arcsec) - columns = ['source_id', 'SOURCE_ID', 'ra', 'dec', 'pmra', 'pmdec', 'parallax', 'radial_velocity', 'ref_epoch','ra_error','dec_error','parallax_error','pmra_error','pmdec_error','radial_velocity_error','ra_dec_corr','ra_parallax_corr','ra_pmra_corr','ra_pmdec_corr','dec_parallax_corr','dec_pmra_corr','dec_pmdec_corr','parallax_pmra_corr','parallax_pmdec_corr','pmra_pmdec_corr'] - target_gaia = Gaia.query_object_async(coordinate=coord, width=width, height=height, columns=columns) - if 'source_id' in target_gaia.columns: - target_gaia = target_gaia[target_gaia['source_id']==gaia_id] + # target coordinates + if 'RA' in simbad_query.columns: + coord = SkyCoord(ra=simbad_query['RA'][0], + dec=simbad_query['DEC'][0], + unit=(u.hourangle, u.degree), frame='icrs') else: - target_gaia = target_gaia[target_gaia['SOURCE_ID']==gaia_id] - self.gaia_epoch = target_gaia['ref_epoch'][0] + coord = SkyCoord(ra=simbad_query['ra'][0], + dec=simbad_query['dec'][0], + unit=(u.hourangle, u.degree), frame='icrs') + + columns = [ + 'source_id', + 'ra', + 'dec', + 'pmra', + 'pmdec', + 'parallax', + 'radial_velocity', + 'ref_epoch', + 'ra_error', + 'dec_error', + 'parallax_error', + 'pmra_error', + 'pmdec_error', + 'radial_velocity_error', + 'ra_dec_corr', + 'ra_parallax_corr', + 'ra_pmra_corr', + 'ra_pmdec_corr', + 'dec_parallax_corr', + 'dec_pmra_corr', + 'dec_pmdec_corr', + 'parallax_pmra_corr', + 'parallax_pmdec_corr', + 'pmra_pmdec_corr' + ] + + col_str = '' + for i, item in enumerate(columns): + if i != 0: + col_str += ',' + col_str += item + + # query target in Gaia + gaia_query = f""" + SELECT {col_str} + FROM gaia{self.gaia_release.lower()}.gaia_source + WHERE source_id = {gaia_id} + """ + gaia_job = Gaia.launch_job_async(gaia_query, dump_to_file=False, verbose=False) + target_gaia = gaia_job.get_results() + + if 'REF_EPOCH' in target_gaia.columns: + self.gaia_epoch = target_gaia['REF_EPOCH'][0] + else: + self.gaia_epoch = target_gaia['ref_epoch'][0] + print(f'[BACKTRACKS INFO]: gathered Gaia {self.gaia_release} data for {self.target_name}') print(f' * Gaia source ID = {gaia_id}') print(f' * Reference epoch = {self.gaia_epoch}') @@ -368,8 +418,8 @@ def query_astrometry(self, nearby_window: float = 0.5): print(f' * PM RA = {target_gaia["pmra"][0]:.2f} mas/yr') print(f' * PM Dec = {target_gaia["pmdec"][0]:.2f} mas/yr') print(f' * Parallax = {target_gaia["parallax"][0]:.2f} mas') - if isinstance(target_gaia["radial_velocity"][0], np.ma.core.MaskedConstant): - print(f' * RV = --') # this addresses the FutureWarning: Format strings passed to MaskedConstant are ignored, but in future may error or produce different behavior + if np.ma.is_masked(target_gaia["radial_velocity"][0]): + print(f' * RV = --') else: print(f' * RV = {target_gaia["radial_velocity"][0]:.2f} km/s') @@ -492,6 +542,13 @@ def prior_transform(self, u): # unpacking unit cube samples ra, dec, pmra, pmdec, par = param + elif len(param) == 2: + # unpacking unit cube samples + ra, dec = param + pmra = None + pmdec = None + par = None + else: # unpacking unit cube samples ra, dec, pmra, pmdec = param @@ -505,13 +562,14 @@ def prior_transform(self, u): # uniform priors for proper motion pmra = transform_uniform(pmra, self.mu_pmra-(10*self.sigma_pmra), self.mu_pmra+(10*self.sigma_pmra)) pmdec = transform_uniform(pmdec, self.mu_pmdec-(10*self.sigma_pmdec), self.mu_pmdec+(10*self.sigma_pmdec)) - else: + elif pmra is not None: # normal priors for proper motion pmra = transform_normal(pmra, self.mu_pmra, self.sigma_pmra) pmdec = transform_normal(pmdec, self.mu_pmdec, self.sigma_pmdec) if self.relax_par_priors: par = transform_uniform(par, 1e-2, self.paro) + # par = transform_uniform(par, 0., 2.) elif par is not None: # ndim = 5 or ndim = 11 # the PPF of Bailer-Jones 2015 eq. 17 @@ -524,8 +582,10 @@ def prior_transform(self, u): param = ra, dec, pmra, pmdec, par, ra_host, dec_host, pmra_host, pmdec_host, par_host, rv_host elif len(param) == 5: param = ra, dec, pmra, pmdec, par - else: + elif len(param) == 4: param = ra, dec, pmra, pmdec + else: + param = ra, dec return param diff --git a/backtracks/plotting.py b/backtracks/plotting.py index d0b08c2..9f751f5 100644 --- a/backtracks/plotting.py +++ b/backtracks/plotting.py @@ -112,9 +112,17 @@ def posterior(backtracks, fileprefix='./', filepost='.pdf'): labels = ["RA (deg, ep=2016)", "DEC (deg, ep=2016)", "pmra (mas/yr)", "pmdec (mas/yr)", "parallax (mas)"] levels = 1.0 - np.exp(-0.5 * np.arange(1, 3.1, 1) ** 2) # 1,2,3 sigma levels for a 2d gaussian # plot extended run (right) + if backtracks.results.samples.shape[1] == 2: + labels = labels[:2] + dims = range(2) + elif backtracks.results.samples.shape[1] == 4: + labels = labels[:4] + dims = range(4) + else: + dims = range(5) fig, ax = dyplot.cornerplot(backtracks.results, color='cornflowerblue', - dims=range(5), + dims=dims, # truths=np.zeros(ndim), # span=[(-4.5, 4.5) for i in range(ndim)], labels=labels, @@ -140,7 +148,7 @@ def posterior(backtracks, fileprefix='./', filepost='.pdf'): axis.set(xlabel=None, ylabel=None) for axis_col in ax[:,-1]: axis_col.xaxis.set_major_formatter(tick_formatter) - ax[-1,-1].set_xscale('log') + # ax[-1,-1].set_xscale('log') # ax[-1,-1].set_xlim(xmin=1e-2) target_name = backtracks.target_name.replace(' ', '_') @@ -397,6 +405,10 @@ def neighborhood(backtracks, fileprefix='./', filepost='.pdf'): nearby_table['parallax'] = backtracks.nearby['parallax'].data truths = backtracks.run_median[2:5] + if len(truths) == 0: + truths = np.array([0., 0., 0.]) + elif len(truths) == 2: + truths = np.append(truths, 0.) # 1,2,3 sigma levels for a 2d gaussian levels = 1.0 - np.exp(-0.5 * np.arange(1, 3.1, 1) ** 2) diff --git a/backtracks/utils.py b/backtracks/utils.py index 91d32f4..d71ff3c 100644 --- a/backtracks/utils.py +++ b/backtracks/utils.py @@ -37,7 +37,13 @@ def radecdists(backtracks, days, param): # for multiple epochs jd_start, jd_end, number = ephem_open() # can't we do this in the System class? - if len(param) == 4: + if len(param) == 2: + ra, dec = param + par = 0.0 + pmra = 0.0 + pmdec = 0.0 + host_icrs=backtracks.host_icrs + elif len(param) == 4: ra, dec, pmra, pmdec = param par=0 host_icrs=backtracks.host_icrs diff --git a/pyproject.toml b/pyproject.toml index 83b82de..7c886e3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,7 @@ authors = [ ] description = "Python package to fit relative astrometry with background star motion tracks." readme = "PYPI_SAFE_README.rst" -requires-python = ">=3.9,<3.12" +requires-python = ">=3.9,<3.13" dynamic = ["dependencies"] classifiers = [ "Programming Language :: Python :: 3", diff --git a/requirements.txt b/requirements.txt index f627060..8153f2a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ astropy -astroquery<=0.4.7 +astroquery corner dynesty matplotlib From 34b1c9f1ecd7d41c61db67f7c34765b46b9c2b12 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Thu, 13 Feb 2025 08:03:43 +0100 Subject: [PATCH 2/2] Removed leftover --- backtracks/backtracks.py | 1 - 1 file changed, 1 deletion(-) diff --git a/backtracks/backtracks.py b/backtracks/backtracks.py index 6a67b4a..fba94a7 100644 --- a/backtracks/backtracks.py +++ b/backtracks/backtracks.py @@ -569,7 +569,6 @@ def prior_transform(self, u): if self.relax_par_priors: par = transform_uniform(par, 1e-2, self.paro) - # par = transform_uniform(par, 0., 2.) elif par is not None: # ndim = 5 or ndim = 11 # the PPF of Bailer-Jones 2015 eq. 17