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 objects that may be directly loaded from FITS headers - if they contain it, of course!

# 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 (stdpipe.astrometry.get_frame_center()) or pixel scale (stdpipe.astrometry.get_pixscale()).

# 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))

Astrometric refinement¶

Existing approximate astrometric solution may be further improved to better represent image distortions. STDPipe provides several methods, all accessible through the higher-level wrapper stdpipe.pipeline.refine_astrometry():

method='quadhash' (default, recommended) - pure Python quad-hash pattern matching via stdpipe.astrometry_quad.refine_wcs_quadhash()
method='scamp' - external SCAMP binary via stdpipe.astrometry.refine_wcs_scamp()
method='astropy' - simple SIP fitting via AstroPy’s fit_wcs_from_points, through stdpipe.astrometry.refine_wcs_simple()
method='astrometrynet' - SIP fitting via Astrometry.Net’s fit-wcs binary, through stdpipe.astrometry.refine_wcs_simple()

stdpipe.pipeline.refine_astrometry(obj, cat, sr=0.002777777777777778, wcs=None, order=0, cat_col_mag='V', cat_col_mag_err=None, cat_col_ra='RAJ2000', cat_col_dec='DEJ2000', cat_col_ra_err='e_RAJ2000', cat_col_dec_err='e_DEJ2000', n_iter=3, use_photometry=True, min_matches=5, method='quadhash', update=True, refine_residual_field=False, residual_field_kwargs=None, verbose=False, **kwargs)[source]

Higher-level astrometric refinement routine.

Dispatches to quad-hash, SCAMP, astropy, or astrometrynet fitting depending on method.

Parameters:

objastropy.table.Table: List of objects on the frame, must contain at least x, y, and flux columns.
catastropy.table.Table or str: Reference astrometric catalogue.
srfloat, optional: Matching radius in degrees.
wcsastropy.wcs.WCS, optional: Initial WCS solution.
orderint, optional: Polynomial order for SIP or PV distortion solution. Default 0 for TAN; 2 is recommended for quadhash.
cat_col_magstr, optional: Catalogue column name for magnitude.
cat_col_mag_errstr, optional: Catalogue column name for magnitude error.
cat_col_rastr, optional: Catalogue column name for Right Ascension.
cat_col_decstr, optional: Catalogue column name for Declination.
cat_col_ra_errstr, optional: Catalogue column name for Right Ascension error.
cat_col_dec_errstr, optional: Catalogue column name for Declination error.
n_iterint, optional: Number of iterations for Python-based matching.
use_photometrybool, optional: Use photometry-assisted matching in Python-based methods.
min_matchesint, optional: Minimum number of good matches required.
methodstr, optional: Fitting method. One of 'quadhash' (default, 2–7× more accurate), 'scamp', 'astropy', or 'astrometrynet'.
updatebool, optional: If True, update ra and dec columns in obj in-place. When combined with refine_residual_field=True, x/y of obj are also corrected (their residual against the catalogue is subtracted) before ra/dec are recomputed from the refined positions.
refine_residual_fieldbool, optional: If True, after the WCS fit, fit a smooth (dx, dy) correction field on the matched obj–cat residuals using stdpipe.astrometry.refine_positions_from_catalog() and (if update=True) apply the inverse correction to obj['x'], obj['y'] in place.
residual_field_kwargsdict, optional: Forwarded to refine_positions_from_catalog(). Common keys: backend, image_shape, grid_shape, min_per_cell, smooth_sigma.
verbosebool or callable, optional: Whether to show verbose messages. May be boolean or a print-like callable.
**kwargs: All other keyword arguments are passed to the respective refinement function.

Returns:

astropy.wcs.WCS or None: Refined astrometric solution, or None on failure.

Quad-hash refinement (recommended)¶

The default and recommended method is quad-hash pattern matching, implemented in stdpipe.astrometry_quad.refine_wcs_quadhash(). It is a pure Python implementation with no external dependencies (only numpy, scipy, astropy).

Key features:

Pure Python - no external binaries needed
Robust pattern matching - uses geometric hashing of star quadrilaterals for reliable matching even with large initial WCS errors
High accuracy - typically 2-7x more accurate than SCAMP, with sub-arcsecond residuals
SIP distortion fitting - supports polynomial distortion orders 1-3
Projection-independent - works with any WCS projection (TAN, ZEA, SIN, STG, ARC, etc.), not just TAN
ZPN projection for wide fields - for images with FoV > 5°, the zenithal polynomial (ZPN) projection with PV radial terms plus optional SIP corrections gives lower residuals than TAN-SIP
Iterative refinement - affine re-matching and progressive sigma-clipping for robust outlier rejection

# Default method - quad-hash refinement
wcs = pipeline.refine_astrometry(obj, cat, 5*pixscale, wcs=wcs,
             cat_col_mag='rmag', verbose=True)

# Equivalent explicit call
from stdpipe.astrometry_quad import refine_wcs_quadhash
wcs = refine_wcs_quadhash(obj, cat, wcs=wcs, sr=10/3600,
             order=2, sn=5, verbose=True)

# Wide-field images (FoV > 5 degrees): use ZPN projection
wcs = refine_wcs_quadhash(obj, cat, wcs=wcs, sr=30/3600,
             order=2, projection='ZPN', pv_deg=5,
             sn=5, 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))

stdpipe.astrometry_quad.refine_wcs_quadhash(obj: Table, cat: Table, wcs: WCS = None, header=None, sr: float = 0.002777777777777778, order: int = 2, projection: str = None, pv_deg: int = 5, cat_col_ra: str = 'RAJ2000', cat_col_dec: str = 'DEJ2000', cat_col_mag: str = 'rmag', obj_col_mag: str = None, obj_col_flux: str = 'flux', sn: float = None, n_det: int = 150, n_ref: int = 600, max_quads_tested: int = 8000, allow_reflection: bool = True, use_wcs_prior: bool = True, scale_tolerance: float = 0.3, enforce_parity: bool = True, refine_use_all: bool = False, refine_n_det: int | None = None, refine_n_ref: int | None = None, refine_match_radius_arcsec: float | None = None, get_header: bool = False, update: bool = False, verbose: bool = False) → WCS[source]

Refine WCS using quad-hash pattern matching algorithm.

Pure Python implementation of astrometric refinement using quad-hash based pattern matching with no external dependencies (only numpy, scipy, astropy). Typically achieves sub-arcsecond accuracy; ~2–7× more accurate than SCAMP, especially for challenging conditions, though ~4× slower.

Parameters:

objastropy.table.Table

List of objects on the frame, must contain at least x, y, and either flux or mag columns.

catastropy.table.Table

Reference astrometric catalogue with RA, Dec, and magnitude columns.

wcsastropy.wcs.WCS, optional

Initial WCS solution (rough estimate).

headerastropy.io.fits.Header, optional

FITS header containing the initial astrometric solution (alternative to wcs).

srfloat, optional

Matching radius in degrees for stage 2 matching.

orderint, optional

SIP polynomial order for the distortion solution (0–5). For TAN projection, controls the SIP polynomial degree. For ZPN projection, controls additional SIP corrections on top of the PV radial model (0 means ZPN-only, 2 is recommended for wide-field images). Ignored for other projections.

projectionstr or None, optional

Target projection type for the output WCS. If None (default), uses the same projection as the input WCS. Supported values:

'TAN' — gnomonic (TAN) with SIP distortion polynomials
'ZPN' — zenithal polynomial with PV radial terms (+ optional SIP for non-radial corrections). Best for wide-field images (FoV > 5°) as ZPN naturally handles radial distortion that TAN-SIP struggles with.
Any other FITS WCS projection code supported by astropy ('STG', 'ARC', 'ZEA', 'SIN', etc.) — linear WCS fit only (no distortion polynomials).

The initial WCS (wcs) is always used for the pattern matching phase; the projection conversion only affects the final WCS fit.

pv_degint, optional

ZPN PV polynomial degree (PV2_0 … PV2_pv_deg). Only used when projection='ZPN'. Default is 5, which is sufficient for most optical systems.

cat_col_rastr, optional

Catalogue column name for Right Ascension. Default is 'RAJ2000'.

cat_col_decstr, optional

Catalogue column name for Declination. Default is 'DEJ2000'.

cat_col_magstr, optional

Catalogue column name for magnitude. Default is 'rmag'.

obj_col_magstr, optional

Object list column name for magnitude. Auto-detected if not provided.

obj_col_fluxstr, optional

Object list column name for flux. Used if magnitude column is absent. Default is 'flux'.

snfloat, optional

If provided, only objects with S/N exceeding this value will be used.

n_detint, optional

Number of brightest detections to use for matching. Default is 150.

n_refint, optional

Number of brightest catalogue stars to use for matching. Default is 600.

max_quads_testedint, optional

Maximum number of detection quads to test. Default is 8000.

allow_reflectionbool, optional

Allow reflected quad matches. Default is True.

use_wcs_priorbool, optional

Use the initial WCS to filter hypotheses by scale/parity. Default is True.

scale_tolerancefloat, optional

Fractional tolerance for scale filtering. Default is 0.30 (±30%).

enforce_paritybool, optional

Enforce parity consistency with the initial WCS. Default is True.

refine_use_allbool, optional

If True, expand the final refinement to a larger pool of objects/catalogue entries.

refine_n_detint or None, optional

Max number of detections in the expanded refinement pool. None means all.

refine_n_refint or None, optional

Max number of catalogue stars in the expanded refinement pool. None means all.

refine_match_radius_arcsecfloat or None, optional

Matching radius in arcsec for the expanded refinement pool. Defaults to the stage 2 radius.

get_headerbool, optional

If True, return the FITS header object instead of WCS solution.

updatebool, optional

If True, update object sky coordinates in-place using the refined WCS.

verbosebool or callable, optional

Whether to show verbose messages. May be boolean or a print-like callable.

Returns:

astropy.wcs.WCS or astropy.io.fits.Header or None: Refined WCS solution, or FITS header if get_header=True, or None on failure.

Examples

>>> from stdpipe.astrometry_quad import refine_wcs_quadhash
>>> wcs_refined = refine_wcs_quadhash(
...     obj, cat, wcs_init,
...     sr=10/3600,  # 10 arcsec matching radius
...     order=2,     # quadratic SIP distortion
...     sn=5,        # use only S/N > 5 objects
...     verbose=True
... )

For wide-field images (FoV > 5 degrees), ZPN projection with SIP corrections typically gives lower systematic residuals at field edges:

>>> wcs_refined = refine_wcs_quadhash(
...     obj, cat, wcs_init,
...     sr=30/3600,
...     order=2,           # SIP order on top of ZPN radial model
...     projection='ZPN',  # zenithal polynomial projection
...     pv_deg=5,          # PV radial polynomial degree
...     verbose=True
... )

SCAMP refinement¶

Alternatively, you may use SCAMP through stdpipe.astrometry.refine_wcs_scamp(). This requires external SCAMP binary to be installed.

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

# Use SCAMP for astrometric refinement
wcs = pipeline.refine_astrometry(obj, cat, 5*pixscale, wcs=wcs,
             method='scamp', cat_col_mag='rmag', verbose=True)

stdpipe.astrometry.refine_wcs_scamp(obj, cat=None, wcs=None, header=None, sr=0.0005555555555555556, order=3, cat_col_ra='RAJ2000', cat_col_dec='DEJ2000', cat_col_ra_err='e_RAJ2000', cat_col_dec_err='e_DEJ2000', cat_col_mag='rmag', cat_col_mag_err='e_rmag', cat_mag_lim=99, sn=None, position_maxerr=None, posangle_maxerr=None, pixscale_maxerr=None, match_flipped=False, extra={}, get_header=False, update=False, _workdir=None, _tmpdir=None, _exe=None, verbose=False)[source]

Wrapper for running SCAMP to get a refined astrometric solution.

Parameters:

objastropy.table.Table: List of objects on the frame, must contain at least x, y, and flux columns.
catastropy.table.Table or str, optional: Reference astrometric catalogue, or a SCAMP network catalogue name (e.g. 'GAIA-DR2').
wcsastropy.wcs.WCS, optional: Initial WCS solution.
headerastropy.io.fits.Header, optional: FITS header containing the initial astrometric solution.
srfloat, optional: Matching radius in degrees. Controls SCAMP’s CROSSID_RADIUS.
orderint, optional: Polynomial order for PV distortion solution (1 or greater).
cat_col_rastr, optional: Catalogue column name for Right Ascension.
cat_col_decstr, optional: Catalogue column name for Declination.
cat_col_ra_errstr, optional: Catalogue column name for Right Ascension error.
cat_col_dec_errstr, optional: Catalogue column name for Declination error.
cat_col_magstr, optional: Catalogue column name for magnitude.
cat_col_mag_errstr, optional: Catalogue column name for magnitude error.
cat_mag_limfloat or sequence of float, optional: Magnitude limit(s) for catalogue stars. A single value is used as an upper limit; two values [min, max] define a range.
snfloat or sequence of float, optional: S/N threshold(s) for SCAMP (passed as SN_THRESHOLDS).
position_maxerrfloat, optional: Maximum positional uncertainty of the initial WCS in arcminutes. Controls the search range for field center offset during pattern matching. Default is SCAMP’s built-in default of 1 arcmin.
posangle_maxerrfloat, optional: Maximum position angle uncertainty in degrees. Default is SCAMP’s built-in default of 5 degrees.
pixscale_maxerrfloat, optional: Maximum pixel scale uncertainty as a multiplicative factor (>1.0). E.g. 1.2 means ±20% around the initial scale. Default is SCAMP’s built-in default of 1.2.
match_flippedbool, optional: If True, also try matching with flipped (mirrored) axes. Default is False.
extradict, optional: Additional parameters to pass to the SCAMP binary.
get_headerbool, optional: If True, return the FITS header object instead of WCS solution.
updatebool, optional: If True, update object sky coordinates in-place using the new WCS.
_workdirstr, optional: If specified, all temporary files will be kept in this directory after the run.
_tmpdirstr, optional: If specified, temporary files will be created inside this path.
_exestr, optional: Full path to the SCAMP executable. Auto-detected from PATH if not provided.
verbosebool or callable, optional: Whether to show verbose messages. May be boolean or a print-like callable.

Returns:

astropy.wcs.WCS or astropy.io.fits.Header or None: Refined astrometric solution, or FITS header if get_header=True, or None on failure.

Simple SIP fitting (astropy / astrometrynet)¶

stdpipe.astrometry.refine_wcs_simple() provides basic SIP distortion fitting using either AstroPy’s fit_wcs_from_points (method='astropy') or Astrometry.Net’s fit-wcs binary (method='astrometrynet'). These methods perform a straightforward polynomial fit to pre-matched object/catalogue positions without the robust pattern matching of quad-hash or the full astrometric model of SCAMP.

# AstroPy-based SIP fitting
wcs = pipeline.refine_astrometry(obj, cat, 5*pixscale, wcs=wcs,
             method='astropy', order=2, cat_col_mag='rmag', verbose=True)

# Using Astrometry.Net's fit-wcs binary (requires local installation)
wcs = pipeline.refine_astrometry(obj, cat, 5*pixscale, wcs=wcs,
             method='astrometrynet', order=2, cat_col_mag='rmag', verbose=True)

stdpipe.astrometry.refine_wcs_simple(obj, cat, order=2, match=True, sr=0.0008333333333333334, update=False, cat_col_ra='RAJ2000', cat_col_dec='DEJ2000', method='astropy', _tmpdir=None, verbose=False)[source]

Refine the WCS using detected objects and a reference catalogue.

Parameters:

objastropy.table.Table: List of detected objects with at least x, y, ra, dec columns.
catastropy.table.Table: Reference astrometric catalogue.
orderint, optional: SIP polynomial order for distortion solution.
matchbool, optional: If True, perform nearest-neighbor matching between objects and catalogue. If False, assume they are already matched line by line.
srfloat, optional: Matching radius in degrees. Default is 3/3600 (3 arcseconds).
updatebool, optional: If True, update object sky coordinates in-place using the new WCS.
cat_col_rastr, optional: Catalogue column name for Right Ascension.
cat_col_decstr, optional: Catalogue column name for Declination.
methodstr, optional: Fitting method. Either 'astropy' or 'astrometrynet'.
_tmpdirstr, optional: If specified, temporary files will be created inside this path.
verbosebool or callable, optional: Whether to show verbose messages. May be boolean or a print-like callable.

Returns:

astropy.wcs.WCS or None: Refined WCS solution, or None on failure.

Residual-field correction¶

After the WCS fit, matched obj–cat pairs typically still show systematic per-position residuals: a few hundredths of a pixel from un-modelled distortions, a residual tilt the polynomial could not absorb, or PSF-position biases that the WCS solver cannot see. STDPipe can fit a smooth two-component (dx, dy) correction field on top of the WCS and either return it as a callable or apply it in place.

The primitive lives in stdpipe.astrometry:

fit_astrometric_residuals() – raw API: take observed and WCS-predicted pixel positions of matched sources, return a callable correct(x, y) -> (x_corrected, y_corrected).
refine_positions_from_catalog() – convenience that takes a match dict (the output of calibrate_photometry() or match()), pulls out the matched positions, and forwards to the primitive. Returns (correct, info) where info carries pre/post-correction median and 90th-percentile residual magnitudes.

Both default to a fast bilinear-grid backend (backend='grid') suitable for dense forced-photometry workloads (~600–1000x faster at prediction than LOESS); pass backend='loess' for higher-quality smoothing when prediction is needed at modest numbers of points. The underlying smoother is stdpipe.smoothing.fit_vector_field_2d().

# Build the correction from an already-matched obj/cat pair
from stdpipe import astrometry, pipeline

m = pipeline.calibrate_photometry(obj, cat, sr=2/3600, wcs=wcs,
                                  cat_col_mag='rmag', verbose=True)
correct, info = astrometry.refine_positions_from_catalog(
    m, obj, cat, wcs,
    image_shape=image.shape,
    grid_shape=(8, 6),       # binned-grid resolution
)
print('median |Δ| {:.3f} → {:.3f} px, q90 {:.3f} → {:.3f} px'.format(
    info['raw_median_dr_pix'], info['corrected_median_dr_pix'],
    info['raw_q90_dr_pix'], info['corrected_q90_dr_pix']))

# Use the correction to refine catalogue-projected positions e.g. for
# forced-position aperture photometry of every catalogue source.
x_pred, y_pred = wcs.all_world2pix(cat['ra'], cat['dec'], 0)
x_corr, y_corr = correct(x_pred, y_pred)

stdpipe.astrometry.fit_astrometric_residuals(x_obs, y_obs, x_pred, y_pred, *, image_shape=None, backend='grid', **smoother_kwargs)[source]

Fit a smooth (dx, dy) correction field for astrometric residuals.

Given matched pairs of observed and WCS-predicted pixel positions of the same sources, fit a smooth two-component field that maps any predicted position to the corresponding observed one. The returned callable applies the additive correction at arbitrary positions.

Parameters:

x_obs, y_obsarray-like, shape (N,): Observed pixel positions of matched sources (e.g. SEP centroids).
x_pred, y_predarray-like, shape (N,): WCS-predicted pixel positions of the same sources from a reference catalogue. x_obs - x_pred is the residual being modelled.
image_shape(H, W), optional: Forwarded to the grid backend; ignored by LOESS.
backend{“grid”, “loess”}: Smoothing backend; see stdpipe.smoothing.fit_vector_field_2d(). Default "grid" because the typical astrometric-residual workflow evaluates the correction at far more positions than were used to fit it (e.g. forced-position photometry of large catalogues), where grid lookup is ~600–1000x faster than LOESS prediction.
**smoother_kwargs: Forwarded to fit_vector_field_2d(). Grid backend keys: grid_shape, min_per_cell, smooth_sigma. LOESS backend keys: scales, k, robust_iters, …

Returns:

correctcallable: correct(x, y) -> (x_corrected, y_corrected): adds the smooth (dx, dy) field at (x, y). Has a .smoother attribute holding the underlying field model for diagnostics.

stdpipe.astrometry.refine_positions_from_catalog(match, obj, cat, wcs, *, obj_col_x='x', obj_col_y='y', cat_col_ra='ra', cat_col_dec='dec', image_shape=None, backend='grid', **smoother_kwargs)[source]

Build an astrometric residual correction from a match dict.

Wraps fit_astrometric_residuals() with the column-extraction boilerplate needed when the matches come from stdpipe.pipeline.calibrate_photometry() or stdpipe.photometry.match(). Returns the correction callable plus a small dictionary of pre/post-correction residual statistics.

Parameters:

matchdict: Cross-match dictionary with at least the keys idx (boolean mask of accepted matches), oidx (object indices into obj), and cidx (catalogue indices into cat).
obj, catastropy.table.Table: Object and catalogue tables matched by match.
wcsastropy.wcs.WCS: WCS used to project catalogue positions into pixel coordinates.
obj_col_x, obj_col_ystr: Column names for observed pixel positions in obj.
cat_col_ra, cat_col_decstr: Column names for sky positions in cat.
image_shape(H, W), optional: Forwarded to the smoothing backend.
backend{“grid”, “loess”}: Smoothing backend.
**smoother_kwargs: Forwarded to fit_astrometric_residuals().

Returns:

correctcallable: Same as fit_astrometric_residuals().
infodict: n_matched, n_used, raw_median_dr_pix, corrected_median_dr_pix, raw_q90_dr_pix, corrected_q90_dr_pix. The corrected statistics are evaluated at the same matched positions used for the fit (in-sample), so they over-state the gain; for an unbiased estimate, subset match['idx'] into a holdout mask before calling this.

Pipeline integration¶

refine_astrometry() accepts opt-in refine_residual_field=True and residual_field_kwargs= parameters. When enabled, after the WCS fit the function re-matches obj to cat positionally via spherical_match(), fits the residual field, and – when update=True – subtracts the smooth correction from obj['x']/obj['y'] in place before re-deriving obj['ra']/obj['dec'] from the refined positions. obj['x']/['y'] then lie on the WCS-consistent grid and downstream code can work with the WCS alone, without carrying the correction model.

wcs = pipeline.refine_astrometry(
    obj, cat, sr=2/3600, wcs=wcs, order=2, cat_col_mag='rmag',
    refine_residual_field=True,
    residual_field_kwargs=dict(
        image_shape=image.shape,
        grid_shape=(8, 6),
        min_per_cell=6, smooth_sigma=1.0,
    ),
    verbose=True,
)

A short diagnostic line is logged with the number of matches used and the median/90th-percentile residual magnitudes before and after the correction. The block silently skips when fewer than 50 positional matches survive, leaving obj untouched and the WCS unchanged.

Astrometric calibration¶

Initial astrometric solution (blind solving)¶

Astrometric refinement¶

Quad-hash refinement (recommended)¶

SCAMP refinement¶

Simple SIP fitting (astropy / astrometrynet)¶

Residual-field correction¶

Pipeline integration¶

STDPipe

Navigation

Related Topics