Commit 424599d4 authored by karpov-sv's avatar karpov-sv
Browse files

Documentation for astrometric calibration added

parent 1c45bcc5
Astrometric calibration
=======================
Astrometric solution is necessary for any reasonable analysis of astronomical imaging data. It is usually being represented as an `Astropy World Coordinate System <https://docs.astropy.org/en/stable/wcs/index.html>`_ objects that may be directly loaded from FITS headers - if they contain it, of course!
.. code-block:: python
# Load FITS header from file
header = fits.getheader(filename)
# Construct WCS structure from this header
wcs = WCS(header)
We have some simple convenience routines that may be used to quickly extract some information from it - like image center and size (:func:`stdpipe.astrometry.get_frame_center`) or pixel scale (:func:`stdpipe.astrometry.get_pixscale`).
.. code-block:: python
# Get the center position, size and pixel scale for the image from WCS
center_ra,center_dec,center_sr = astrometry.get_frame_center(wcs=wcs, width=image.shape[1], height=image.shape[0])
pixscale = astrometry.get_pixscale(wcs=wcs)
print('Frame center is %.2f %.2f radius %.2f deg, %.2f arcsec/pixel' % (center_ra, center_dec, center_sr, pixscale*3600))
Initial astrometric solution (blind solving)
--------------------------------------------
If your FITS header does not have any initial astrometric solution, you may derive it using `Astrometry.Net <https://nova.astrometry.net>`_ blind matching algorithm. We have two routines for working with it.
First one, :func:`stdpipe.astrometry.blind_match_astrometrynet`, directly uses their on-line service through `Astroquery <https://astroquery.readthedocs.io/en/latest/astrometry_net/astrometry_net.html>`_ library and requires an API key to be present in your :file:`~/.astropy/config/astroquery.cfg` config file
.. note::
To get API key, you have to sign into your account at https://nova.astrometry.net and then generate the key on your profile page. Then, you should add it to your file :file:`~/.astropy/config/astroquery.cfg` by adding there
.. code-block:: ini
[astrometry_net]
# The Astrometry.net API key
api_key = 'your-api-key-goes-here'
If API key is not set, the function call will produce a warning message and fail.
You may also directly provide API key to the function call using `api_key` argument.
The function operates on the list of detected objects, so in principle network stress from its use is quite small (it will not try to upload your images to remote servers!). To speed up the solution, you may directly specify approximate center coordinates of the field center, as well as lower and upper bounds on the solution pixel scale.
Second routine, :func:`stdpipe.astrometry.blind_match_objects`, requires `local installation of Astrometry.Net binaries and index files <http://astrometry.net/use.html>`_. The routine itself has mostly the same signature as previous one, with some additional options to specify the configuration of the solver in finer details.
Example of using the code for solving the astrometry if not set in the header:
.. code-block:: python
if wcs is None or not wcs.is_celestial:
print('WCS is not set, blind solving for it...')
# Try to guess rough frame center
ra0 = header.get('RA')
dec0 = header.get('DEC')
sr0 = 1.0 # Should be enough?..
# Should be sufficient for most images
pixscale_low = 0.3
pixscale_upp = 3
# We will use no more than 500 brightest objects with S/N>10 to solve
wcs = astrometry.blind_match_astrometrynet(obj[:500], center_ra=ra0, center_dec=dec0, radius=sr0, scale_lower=pixscale_low, scale_upper=pixscale_upp, sn=10)
# .. or the same using local solver with custom config
wcs = astrometry.blind_match_objects(obj[:500], center_ra=ra0, center_dec=dec0, radius=sr0, scale_lower=pixscale_low, scale_upper=pixscale_upp, sn=10, verbose=True, _exe='/home/karpov/astrometry/bin/solve-field', config='/home/karpov/astrometry/etc/astrometry-2mass.cfg')
if wcs is not None and wcs.is_celestial:
print('Blind solving succeeded!')
.. autofunction:: stdpipe.astrometry.blind_match_astrometrynet
:noindex:
.. autofunction:: stdpipe.astrometry.blind_match_objects
:noindex:
Astrometric refinement
----------------------
Existing approximate astrometric solution may be further improved to better represent image distortions. It may be done using `SCAMP <https://github.com/astromatic/scamp>`_ code - we have a routine :func:`stdpipe.astrometry.refine_wcs_scamp` that conveniently wraps it, hiding most of trickier details of running it.
.. attention::
SCAMP, as well as SExtractor and SWarp, uses older and non-standard way of describing the distortions in the image - PV polynomials, while most of other softwares nowadays - including Astrometry.Net - prefer (standard) SIP polynomials. Python WCS loads both nicely, but there is no way to specify which one you wish to save back! Thus, converting one to another is not transparent, and you should also be aware which one you use! E.g. feeding SCAMP with initial WCS containing SIP polynomials will probably work, and it will return PV polynomial back that may be used for SWarping it. But directly running SWarp on Astrometry.Net - generated solutions will give you wrong results!
You may read the discussion of the problem e.g. there - https://github.com/evertrol/sippv
We also have another, less reliable and less tested routine :func:`stdpipe.astrometry.refine_wcs` that will operate with SIP distortions, and do so either in pure Python, or using `fit-wcs` executable from Astrometry.Net installation. We also have a higher-level routine that uniformly wraps all these functions - :func:`stdpipe.pipeline.refine_astrometry`.
Example:
.. code-block:: python
# Let's use SCAMP for astrometric refinement.
wcs = pipeline.refine_astrometry(obj, cat, 5*pixscale, wcs=wcs, method='scamp', cat_col_mag='rmag', verbose=True)
if wcs is not None:
# Update WCS info in the header
astrometry.clear_wcs(header, remove_comments=True, remove_underscored=True, remove_history=True)
header.update(wcs.to_header(relax=True))
.. autofunction:: stdpipe.astrometry.refine_wcs_scamp
:noindex:
.. autofunction:: stdpipe.pipeline.refine_astrometry
:noindex:
......@@ -255,8 +255,8 @@ latex_documents = [
(
master_doc,
"stdpipe.tex",
"stdpipe Documentation",
"The stdpipe Team",
"STDPipe Documentation",
"The STDPipe Team",
"manual",
),
]
......@@ -286,7 +286,7 @@ latex_documents = [
# One entry per manual page. List of tuples
# (source start file, name, description, authors, manual section).
man_pages = [(master_doc, "stdpipe", "stdpipe Documentation", [author], 1)]
man_pages = [(master_doc, "stdpipe", "STDPipe Documentation", [author], 1)]
# If true, show URL addresses after external links.
# man_show_urls = False
......@@ -301,10 +301,10 @@ texinfo_documents = [
(
master_doc,
"stdpipe",
"stdpipe Documentation",
"STDPipe Documentation",
author,
"stdpipe",
"One line description of project.",
"Simple Transient Detection Pipeline",
"Miscellaneous",
),
]
......
......@@ -11,8 +11,6 @@ from astropy.wcs.utils import fit_wcs_from_points
from astropy.coordinates import SkyCoord
from astropy.table import Table
from astroquery.astrometry_net import AstrometryNet
from scipy.stats import chi2
from . import utils
......@@ -90,6 +88,33 @@ def blind_match_objects(obj, order=2, update=False, sn=20, get_header=False,
center_ra=None, center_dec=None, radius=None,
scale_lower=None, scale_upper=None, scale_units='arcsecperpix',
config=None, extra={}, _workdir=None, _tmpdir=None, _exe=None, verbose=False):
"""Thin wrapper for blind plate solving using local Astrometry.Net and a list of detected objects.
It requires `solve-field` binary from Astrometry.Net and some index files to be locally available.
:param obj: List of objects on the frame that should contain at least `x`, `y` and `flux` columns.
:param order: Order for the SIP spatial distortion polynomial
:param update: If set, the object list will be updated in-place to contain correct `ra` and `dec` sky coordinates
:param sn: If provided, only objects with signal to noise ratio exceeding this value will be used for matching.
:param get_header: If True, function will return the FITS header object instead of WCS solution
:param width: Image width, to be used for guessing pixel coordinats of frame center. Optional.
:param height: Image height, to be used for guessing pixel coordinats of frame center. Optional.
:param center_ra: Approximate center RA of the field, degrees. Optional.
:param center_dec: Approximate center Dec of the field, degrees. Optional.
:param radius: If set, the server will look for solutions only within this radius from the center specified above.
:param scale_lower: Optional lower limit for the solution scale.
:param scale_upper: Optional upper limit for the solution scale.
:param scale_units: Units of the `scale_lower`/`scale_upper` parameters. May be one of `arcsecperpix`, `arcminwidth`, or `degwidth`.
:param config: Path to config file for `solve-field`, optional.
:param extra: Dictionary of additional parameters to be passed to `solve-field` binary, optional.
:param _workdir: If specified, all temporary files will be created in this directory, and will be kept intact after running `solve-field`. May be used for debugging exact inputs and outputs of the executable. Optional
:param _tmpdir: If specified, all temporary files will be created in a dedicated directory (that will be deleted after running the executable) inside this path.
:param _exe: Full path to `solve-field` executable. If not provided, the code tries to locate it automatically in your :envvar:`PATH`.
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function.
:returns: Either astrometric solution as astropy.wcs.WCS object, or FITS header if :code:`get_header=True`
"""
'''
Thin wrapper for blind plate solving using local Astrometry.Net and a list of detected objects.
'''
......@@ -234,11 +259,31 @@ def blind_match_astrometrynet(obj, order=2, update=False, sn=20, get_header=Fals
solve_timeout=600, api_key=None,
center_ra=None, center_dec=None, radius=None,
scale_lower=None, scale_upper=None, scale_units='arcsecperpix', **kwargs):
'''
Thin wrapper for remote plate solving using Astrometry.Net and a list of detected objects.
Most of the parameters are passed directly to astroquery.astrometrynet.AstrometryNet.solve_from_source_list routine.
"""Thin wrapper for remote plate solving using Astrometry.Net and a list of detected objects.
Most of the parameters are passed directly to `astroquery.astrometrynet.AstrometryNet.solve_from_source_list` routine.
API key may either be provided as an argument or specified in ~/.astropy/config/astroquery.cfg
'''
:param obj: List of objects on the frame that should contain at least `x`, `y` and `flux` columns.
:param order: Order for the SIP spatial distortion polynomial
:param update: If set, the object list will be updated in-place to contain correct `ra` and `dec` sky coordinates
:param sn: If provided, only objects with signal to noise ratio exceeding this value will be used for matching.
:param get_header: If True, function will return the FITS header object instead of WCS solution
:param width: Image width, to be used for guessing pixel coordinats of frame center. Optional.
:param height: Image height, to be used for guessing pixel coordinats of frame center. Optional.
:param solve_timeout: Timeout in seconds to wait for the solution from remote server.
:param api_key: API key, optional. If not set, it should be configured in your ~/.astropy/config/astroquery.cfg file
:param center_ra: Approximate center RA of the field, degrees. Optional.
:param center_dec: Approximate center Dec of the field, degrees. Optional.
:param radius: If set, the server will look for solutions only within this radius from the center specified above.
:param scale_lower: Optional lower limit for the solution scale.
:param scale_upper: Optional upper limit for the solution scale.
:param scale_units: Units of the `scale_lower`/`scale_upper` parameters. May be one of `arcsecperpix`, `arcminwidth`, or `degwidth`.
:returns: Either astrometric solution as astropy.wcs.WCS object, or FITS header if :code:`get_header=True`
"""
# Import required module - we postpone it until now to hide warnings for API key
from astroquery.astrometry_net import AstrometryNet
# Sort objects according to decreasing flux
aidx = np.argsort(-obj['flux'])
......@@ -360,6 +405,16 @@ def refine_wcs(obj, cat, order=2, match=True, sr=3/3600, update=False,
return wcs
def clear_wcs(header, remove_comments=False, remove_history=False, remove_underscored=False, copy=False):
"""Clears WCS related keywords from FITS header
:param header: Header to operate on
:param remove_comments: Whether to also remove COMMENT keywords
:param remove_history: Whether to also remove HISTORY keywords
:param remove_underscored: Whether to also remove all keywords starting with underscore (often made by e.g. Astrometry.Net)
:param copy: If True, do not change original FITS header
:returns: Modified FITS header
"""
if copy:
header = header.copy()
......@@ -440,11 +495,33 @@ def refine_wcs_scamp(obj, cat=None, wcs=None, header=None, sr=2/3600, order=3,
cat_mag_lim=99, sn=None, extra={},
get_header=False, update=False,
_workdir=None, _tmpdir=None, _exe=None, verbose=False):
"""
Wrapper for running SCAMP on user-provided object list and catalogue
"""Wrapper for running SCAMP on user-provided object list and catalogue to get refined astrometric solution.
:param obj: List of objects on the frame that should contain at least `x`, `y` and `flux` columns.
:param cat: Reference astrometric catalogue
:param wcs: Initial WCS
:param header: FITS header containing initial astrometric solution, optional.
:param sr: Matching radius in degrees
:param order: Polynomial order for PV distortion solution (1 or greater)
:param cat_col_ra: Catalogue column name for Right Ascension
:param cat_col_dec: Catalogue column name for Declination
:param cat_col_ra_err: Catalogue column name for Right Ascension error
:param cat_col_dec_err: Catalogue column name for Declination error
:param cat_col_mag: Catalogue column name for the magnitude in closest band
:param cat_col_mag_err: Catalogue column name for the magnitude error
:param cat_mag_lim: Magnitude limit for catalogue stars
:param sn: If provided, only objects with signal to noise ratio exceeding this value will be used for matching.
:param extra: Dictionary of additional parameters to be passed to SCAMP binary, optional.
:param get_header: If True, function will return the FITS header object instead of WCS solution
:param update: If set, the object list will be updated in-place to contain correct `ra` and `dec` sky coordinates
:param _workdir: If specified, all temporary files will be created in this directory, and will be kept intact after running SCAMP. May be used for debugging exact inputs and outputs of the executable. Optional
:param _tmpdir: If specified, all temporary files will be created in a dedicated directory (that will be deleted after running the executable) inside this path.
:param _exe: Full path to SCAMP executable. If not provided, the code tries to locate it automatically in your :envvar:`PATH`.
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function.
:returns: Refined astrometric solution, or FITS header if :code:`get_header=True`
"""
# Simple wrapper around print for logging in verbose mode only
# Simple wrapper around print for logging in verbose mode only
log = (verbose if callable(verbose) else print) if verbose else lambda *args,**kwargs: None
# Find the binary
......
......@@ -30,10 +30,30 @@ def refine_astrometry(obj, cat, sr=10/3600, wcs=None, order=0,
cat_col_ra_err='e_RAJ2000', cat_col_dec_err='e_DEJ2000',
n_iter=3, use_photometry=True, min_matches=5, method='astropy',
update=True, verbose=False, **kwargs):
"""
Higher-level astrometric refinement routine.
"""
"""Higher-level astrometric refinement routine that may use either SCAMP or pure Python based methods.
:param obj: List of objects on the frame that should contain at least `x`, `y` and `flux` columns.
:param cat: Reference astrometric catalogue
:param sr: Matching radius in degrees
:param wcs: Initial WCS
:param order: Polynomial order for SIP or PV distortion solution
:param cat_col_mag: Catalogue column name for the magnitude in closest band
:param cat_col_mag_err: Catalogue column name for the magnitude error
:param cat_col_ra: Catalogue column name for Right Ascension
:param cat_col_dec: Catalogue column name for Declination
:param cat_col_ra_err: Catalogue column name for Right Ascension error
:param cat_col_dec_err: Catalogue column name for Declination error
:param n_iter: Number of iterations for Python-based matching
:param use_photometry: Use photometry-assisted method in Python-based matching
:param min_matches: Minimal number of good matches in Python-based matching
:param method: May be either 'scamp' or 'astropy' or 'astrometrynet'
:param update: If set, the object list will be updated in-place to contain correct `ra` and `dec` sky coordinates
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function.
:param kwargs: All other parameters will be directly passed to :func:`~stdpipe.astrometry.refine_wcs_scamp`
:returns: Refined astrometric solution
"""
# Simple wrapper around print for logging in verbose mode only
log = (verbose if callable(verbose) else print) if verbose else lambda *args,**kwargs: None
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment