import os
import requests
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table
from stcal.alignment import compute_fiducial
from stcal.tweakreg._s3_catalog import S3_CATALOGS, get_s3_catalog
GSSS_CATALOGS = ["GAIAREFCAT", "GAIADR3", "GAIADR2", "GAIADR1"]
ASTROMETRIC_CAT_ENVVAR = "ASTROMETRIC_CATALOG_URL"
DEF_CAT_URL = "https://gsss.stsci.edu/webservices"
SERVICELOCATION = os.environ.get(ASTROMETRIC_CAT_ENVVAR, DEF_CAT_URL)
TIMEOUT = 30.0 # in seconds
__all__ = ["TIMEOUT", "compute_radius", "create_astrometric_catalog", "get_catalog"]
[docs]
def create_astrometric_catalog(
wcs,
epoch,
catalog="GAIADR3",
output="ref_cat.ecsv",
table_format="ascii.ecsv",
num_sources=None,
timeout=TIMEOUT,
):
"""Create an astrometric catalog that covers the inputs' field-of-view.
Parameters
----------
wcs : `~astropy.wcs.WCS`
WCS object specified by the user as generated by
`resample.resample_utils.make_output_wcs`. This will typically
have the same plate-scale and orientation as the first member in the
list of input images to make_output_wcs. Fortunately, for alignment,
this doesn't matter since no resampling of data will be performed.
epoch : float or None
Reference epoch used to update the coordinates for proper motion
(in decimal year). When ``None`` no proper motion correction will
be performed and all sources (even those without proper motion) will
be returned. When not ``None`` only sources with proper motion
will be returned.
catalog : str, optional
Name of catalog to extract astrometric positions for sources in the
input images' field-of-view. Default: GAIADR3. Options available are
documented on the catalog web page.
output : str, optional
Filename to give to the astrometric catalog read in from the master
catalog web service. If None, no file will be written out.
num_sources : int
Maximum number of brightest/faintest sources to return in catalog.
If `num_sources` is negative, return that number of the faintest
sources. By default, all sources are returned.
timeout : float
Maximum time to wait (in seconds) for the catalog service to respond.
Notes
-----
This function will point to astrometric catalog web service defined
through the use of the ASTROMETRIC_CATALOG_URL environment variable.
Returns
-------
ref_table : `~astropy.table.Table`
Astropy Table object of the catalog
"""
# start by creating a composite field-of-view for all inputs
radius, fiducial = compute_radius(wcs)
# perform query for this field-of-view
ref_dict = get_catalog(
fiducial[0], fiducial[1], epoch=epoch, search_radius=radius, catalog=catalog, timeout=timeout
)
if len(ref_dict) == 0:
return ref_dict
colnames = ("ra", "dec", "mag", "objID", "epoch")
ref_table = ref_dict[colnames]
# Add catalog name as meta data
ref_table.meta["catalog"] = catalog
# rename coordinate columns to be consistent with tweakwcs
ref_table.rename_column("ra", "RA")
ref_table.rename_column("dec", "DEC")
# sort table by magnitude, fainter to brightest
ref_table.sort("mag", reverse=True)
# If specified by the use through the 'num_sources' parameter,
# trim the returned catalog down to the brightest 'num_sources' sources
# Should 'num_sources' be a negative value, it will return the faintest
# 'num_sources' sources.
if num_sources is not None:
indx = -1 * num_sources
ref_table = ref_table[:indx] if num_sources < 0 else ref_table[indx:]
# Write out table to a file, if specified
if output is not None:
ref_table.write(output, format=table_format, overwrite=True)
return ref_table
[docs]
def compute_radius(wcs):
"""Compute the radius from the center to the furthest edge of the WCS."""
fiducial = compute_fiducial(
[
wcs,
]
)
img_center = SkyCoord(ra=fiducial[0] * u.degree, dec=fiducial[1] * u.degree)
wcs_foot = wcs.footprint()
img_corners = SkyCoord(ra=wcs_foot[:, 0] * u.degree, dec=wcs_foot[:, 1] * u.degree)
radius = img_center.separation(img_corners).max().value
return radius, fiducial
[docs]
def get_catalog(
right_ascension,
declination,
epoch=2016.0,
search_radius=0.1,
catalog="GAIADR3",
timeout=TIMEOUT,
):
"""Extract catalog from VO web service.
Parameters
----------
right_ascension : float
Right Ascension (RA) of center of field-of-view (in decimal degrees)
declination : float
Declination (Dec) of center of field-of-view (in decimal degrees)
epoch : float or None, optional
Reference epoch used to update the coordinates for proper motion
(in decimal year). When ``None`` no proper motion correction will
be performed and all sources (even those without proper motion) will
be returned. When not ``None`` only sources with proper motion
will be returned. Default: 2016.0
search_radius : float, optional
Search radius (in decimal degrees) from field-of-view center to use
for sources from catalog. Default: 0.1 degrees
catalog : str, optional
Name of catalog to query, as defined by web-service. Default: 'GAIADR3'
timeout : float, optional
Timeout in seconds to wait for the catalog web service to respond. Default: 30.0 s
Returns
-------
csv : `~astropy.table.Table`
CSV object of returned sources with all columns as provided by catalog
"""
if catalog in S3_CATALOGS:
return get_s3_catalog(
right_ascension,
declination,
epoch,
search_radius,
catalog,
timeout=timeout,
)
service_type = "vo/CatalogSearch.aspx"
headers = {"Content-Type": "text/csv"}
if epoch is None:
spec = f"RA={right_ascension}&DEC={declination}&SR={search_radius}&FORMAT=CSV&CAT={catalog}"
else:
spec = (
f"RA={right_ascension}&DEC={declination}&EPOCH={epoch}"
f"&SR={search_radius}&FORMAT=CSV&CAT={catalog}"
)
service_url = f"{SERVICELOCATION}/{service_type}?{spec}"
try:
rawcat = requests.get(service_url, headers=headers, timeout=timeout)
except requests.exceptions.ConnectionError as err:
raise requests.exceptions.ConnectionError(
"Could not connect to the VO API server. Try again later."
) from err
except requests.exceptions.Timeout as err:
raise requests.exceptions.Timeout("The request to the VO API server timed out.") from err
except requests.exceptions.RequestException as err:
raise requests.exceptions.RequestException("There was an unexpected error with the request.") from err
r_contents = rawcat.content.decode() # convert from bytes to a String
if r_contents.startswith("No data records"):
r_contents = "\n"
catalog = Table.read(r_contents, format="csv", comment="#")
if epoch and catalog:
# When provided with an epoch gsss returns corrected and non-corrected
# sources. Filter out the non-corrected ones.
catalog = catalog[~(catalog["pmra"].mask & catalog["pmdec"].mask)]
return catalog