import os
import requests
from astropy import table
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table
from stcal.alignment import compute_fiducial
ASTROMETRIC_CAT_ENVVAR = "ASTROMETRIC_CATALOG_URL"
DEF_CAT_URL = "http://gsss.stsci.edu/webservices"
SERVICELOCATION = os.environ.get(ASTROMETRIC_CAT_ENVVAR, DEF_CAT_URL)
TIMEOUT = 30.0 # in seconds
"""
Primary function for creating an astrometric reference catalog.
"""
__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
Reference epoch used to update the coordinates for proper motion
(in decimal year).
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, optional
Reference epoch used to update the coordinates for proper motion
(in decimal year). 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
"""
service_type = "vo/CatalogSearch.aspx"
spec_str = "RA={}&DEC={}&EPOCH={}&SR={}&FORMAT={}&CAT={}"
headers = {"Content-Type": "text/csv"}
fmt = "CSV"
spec = spec_str.format(
right_ascension,
declination,
epoch,
search_radius,
fmt,
catalog
)
service_url = f"{SERVICELOCATION}/{service_type}?{spec}"
try:
rawcat = requests.get(service_url, headers=headers, timeout=timeout)
except requests.exceptions.ConnectionError:
raise requests.exceptions.ConnectionError(
"Could not connect to the VO API server. Try again later."
)
except requests.exceptions.Timeout:
raise requests.exceptions.Timeout("The request to the VO API server timed out.")
except requests.exceptions.RequestException:
raise requests.exceptions.RequestException(
"There was an unexpected error with the request."
)
r_contents = rawcat.content.decode() # convert from bytes to a String
if r_contents.startswith("No data records"):
r_contents = "\n"
return Table.read(r_contents, format="csv", comment='#')