World Coordinate System¶

class fitsy.Wcs(header, alt=' ')¶

Bases: object

World Coordinate System for an HDU.

Constructed via FitsFile.wcs(), ImageHdu.wcs(), or directly from a header. Supports celestial, spectral, time, phase and generic linear axes. Recognizes the SIP, TPV, TNX/ZPX, -TAB and DSS distortion conventions.

Examples

>>> with fitsy.open("image.fits") as f:
...     wcs = f[0].wcs()
...     kinds = wcs.axis_kinds()
...     world = wcs.pixel_to_world([512.0, 512.0])
...     ra = world[kinds.index("longitude")]

axis_kinds()¶

Kind of coordinate each axis carries, in axis order.

Use this to find an axis by meaning rather than by position. pixel_to_world() returns one value per axis in the same order, so entry i here names value i there.

Returns:: list of str – One entry per axis, each one of 'longitude', 'latitude', 'spectral', 'time', 'phase', 'stokes' or 'linear'.

Notes

The kind comes from the type half of CTYPEia, so an axis driven by a -TAB lookup still reports its coordinate type. is_tabular() reports the lookup itself.

Examples

Find the spectral axis of a cube and read its world value:

>>> with fitsy.open("cube.fits") as f:
...     wcs = f[0].wcs()
...     kinds = wcs.axis_kinds()
...     world = wcs.pixel_to_world([31.0, 23.0, 4.0])
...     freq = world[kinds.index("spectral")]

celestial_axes()¶

Indices of the celestial axes.

Returns:: tuple of int or None – (lon_axis, lat_axis) (zero-based), or None if no celestial pair is declared.

crval¶: Per-axis CRVAL reference values, in the unit given by cunit.

ctype¶: Per-axis CTYPE strings.

cunit¶: Per-axis CUNIT strings. Empty for an axis whose header carries no CUNIT card; that axis then uses the default unit for its axis type.

footprint()¶

World coordinates of the image’s corner pixels.

Corners are pixel centers, not the outer edge of the grid. For the outer edge, call pixel_to_world() with -0.5 and n - 0.5.

Corners come back in Gray-code order, so consecutive corners differ on one axis alone. A two-axis image therefore yields (0, 0), (nx-1, 0), (nx-1, ny-1), (0, ny-1), which walks the image counter-clockwise in pixel space and closes the ring.

Returns:: numpy.ndarray – Shape (2**k, naxis). k is the number of axes pixel_shape covers, which is naxis for a normal image. A two-axis image gives (4, 2) of (ra, dec) in degrees. A corner the WCS cannot transform comes back as nan. See the notes below.
Raises:: FitsError – In three cases: - pixel_shape is None. A fitted WCS has no image to take corners from. - An axis pixel_shape covers has length zero. - The WCS has more than 16 axes.

Notes

This reports corners, not an axis-aligned bounding box. A rotated image has corners outside the box its own minimum and maximum describe, and an RA axis that crosses zero makes such a box meaningless.

Corners go through pixel_to_world() in its batch form. A corner outside the projection’s domain therefore fills its row with nan instead of raising. A wide-field SIN or AZP image can put every corner outside that domain. The whole array is then nan. Test the result with numpy.isfinite. Pass one corner to pixel_to_world() to read the reason it failed.

WCSAXES may exceed NAXIS. A coordinate axis past the end of pixel_shape then has no length to take a corner from. That axis holds its reference pixel for every corner. The corner count follows the image, and every corner still carries a full world vector.

is_celestial¶: True when the WCS has a celestial axis pair.

is_tabular(axis)¶

Whether an axis takes its coordinate from a -TAB lookup.

Parameters:: axis (int) – Zero-based axis index.
Returns:: bool – True when axis axis is tabular. False for any other axis, and for an index this WCS does not have.

Notes

This is a property of the algorithm, not of the coordinate, so it is independent of axis_kinds(). A tabular axis needs its binary table loaded before it can transform. fitsy.FitsFile.wcs() and fitsy.ImageHdu.wcs() load it. fitsy.Wcs(header) cannot, because a header alone does not reach the table.

naxis¶: Number of WCS axes. May exceed the header’s NAXIS when WCSAXESa sets a higher axis count.

pixel_scale_at(px, py, origin=0)¶

Local pixel scale at (px, py).

Parameters:

px (float) – Pixel coordinates.
py (float) – Pixel coordinates.
origin (int, optional) – 0 (default) treats inputs as 0-based; 1 as 1-based FITS coordinates.

Returns:

tuple of float – Pixel scale in arcseconds per pixel along the two celestial axes.

Raises:

FitsError – If the WCS has no celestial axis pair, if it has unresolved -TAB axes, or if (px, py) or an adjacent pixel used for the finite difference falls outside the projection’s valid domain.

Notes

fitsy measures this by finite difference on the sphere, so the result includes projection distortion and any local skew. It is a great-circle distance, always positive, not the signed CDELT value. An image with flipped RA still reports a positive scale.

pixel_shape¶

Size of the image this WCS came from, in FITS axis order (NAXIS1 first), or None when unknown.

This is a snapshot of the NAXISn cards taken when the WCS was parsed. It is not part of the coordinate description: no transform reads it, and no transform checks a pixel coordinate against it. It is None for a WCS from fit_wcs(), which has no image, or for a header without NAXISn cards. A cropped or rebinned image leaves this value stale.

pixel_to_world(pix, origin=0)¶

Forward transform pixel coordinates to world coordinates.

Accepts one point or many. A length-naxis sequence is one point and returns a list. An (N, naxis) array is N points and returns an (N, naxis) array.

Parameters:

pix (array-like) – Shape (naxis,) for a single point, or (N, naxis) for a batch.
origin (int, optional) – 0 (default) treats pix as 0-based, 1 treats it as 1-based FITS coordinates.

Returns:

list of float or numpy.ndarray – World coordinates with units given by cunit. A list for a single point, an (N, naxis) array for a batch.

Raises:

ValueError – In three cases: - origin is neither 0 nor 1. - pix is neither 1-D nor 2-D. - A batch does not have exactly naxis columns. An (naxis, N) array is the transpose of a batch, not a batch of N points.
FitsError – In three cases: - A single point does not have naxis elements. - The WCS has unresolved -TAB axes. - pix falls outside the projection’s valid domain. This applies to a single point only.

Notes

The two forms differ on a point the WCS cannot transform. The single-point form raises, so the message names the reason. The batch form fills that point with nan and keeps going, because a wide field routinely mixes valid and invalid pixels.

to_header(alt=' ')¶

Serialize this WCS to a fresh Header.

Parameters:: alt (str, optional) – ' ' (default) for the primary description, or 'A' through 'Z' for an alternate.
Returns:: Header – A new header holding every WCS keyword fitsy’s reader understands: the linear pipeline, LONPOLE/LATPOLE and the projection’s PV parameters, SIP, TPV, TNX/ZPX, DSS plate solutions, spectral rest quantities, and the -TAB pointer cards.
Raises:: FitsError – If alt is not ' ' or 'A'-'Z'.

Notes

Parsing the returned header reproduces this WCS. Two things a bare header cannot carry:

NAXISn are written as zero placeholders, because a WCS carries no image dimensions. Merge the result into a header that already has the real values.
A -TAB axis names its lookup table by EXTNAME. That BINTABLE is a separate HDU and must be written alongside this header.

world_to_pixel(world, origin=0)¶

Inverse transform world coordinates to pixel coordinates.

Accepts one point or many, and mirrors pixel_to_world() in both shape handling and error handling.

Parameters:

world (array-like) – Shape (naxis,) for a single point, or (N, naxis) for a batch.
origin (int, optional) – 0 (default) returns 0-based pixel coordinates, 1 returns 1-based FITS coordinates.

Returns:

list of float or numpy.ndarray – Pixel coordinates in the chosen origin. A list for a single point, an (N, naxis) array for a batch.

Raises:

ValueError – In three cases: - origin is neither 0 nor 1. - world is neither 1-D nor 2-D. - A batch does not have exactly naxis columns.
FitsError – In three cases: - A single point does not have naxis elements. - The WCS has unresolved -TAB axes. - The inverse transform does not converge. This applies to a single point only.

Notes

As in pixel_to_world(), the batch form yields nan for a point that does not transform rather than raising.

class fitsy.WcsFit¶

Bases: object

Result of fit_wcs().

Carries the fitted Wcs and per-point residuals.

Variables:

wcs (Wcs) – The fitted world coordinate system.
rms_arcsec (float) – Root-mean-square residual across all reference points (arcsec).
max_arcsec (float) – Largest single-point residual (arcsec).

max_arcsec¶

residuals_arcsec¶: Per-point residuals as a numpy array of shape (N, 2), holding (delta_alpha * cos(delta), delta_dec) in arcseconds.

rms_arcsec¶

wcs¶

fitsy.fit_wcs(pixels, sky, projection='TAN', crpix=None, crval=None, sip_order=None, fit_sip_inverse=True, frame='equatorial', origin=0)¶

Fit a celestial WCS to (pixel, sky) reference correspondences.

Parameters:

pixels (array-like) – Shape (N, 2) array-like of pixel coordinates (numpy array, list of lists, tuple of tuples, etc.).
sky (array-like) – Shape (N, 2) array-like of (ra, dec) in degrees, or more generally (lon, lat) in the chosen frame (numpy array, list of lists, tuple of tuples, etc.).
projection (str, optional) – Three-letter projection code (Paper II Table 13), case-insensitive. One of AZP, SZP, TAN, STG, SIN, ARC, ZPN, ZEA, AIR, CYP, CEA, CAR, MER, SFL, PAR, MOL, AIT, COP, COE, COD, COO, BON, PCO, TSC, CSC, QSC, HPX or XPH. Default "TAN".
crpix (tuple of float, optional) – Pin the reference pixel. Interpreted in the same origin as pixels. Default None, which solves for the reference pixel as part of the fit.
crval (tuple of float, optional) – Pin the tangent point, in degrees. Default None, which uses the spherical centroid of the sky points.
sip_order (int, optional) – Order of a SIP polynomial distortion fit. Valid range is 2 to 9. Default None, which fits no SIP distortion. Orders 0 and 1 are rejected, because CRPIX and the CD matrix already absorb those terms.
fit_sip_inverse (bool, optional) – When sip_order is given, also fit the AP/BP inverse polynomial. Default True.
frame (str, optional) – Celestial frame for the sky coordinates, case-insensitive. 'equatorial' (default), 'icrs', 'fk5' and 'fk4' are synonyms; all four emit the RA--/DEC- CTYPE prefixes. The other accepted values are 'galactic', 'ecliptic', 'supergalactic' and 'helioecliptic'.
origin (int, optional) – 0 (default, numpy/C convention) treats pixels and crpix as 0-based; 1 treats them as 1-based FITS coordinates. The fitted WCS itself always carries 1-based CRPIX values per the FITS standard.

Returns:

WcsFit – Fitted WCS, residuals, and summary statistics.

Raises:

ValueError – If pixels or sky is not an (N, 2) array, if they have different numbers of rows, or if frame names none of the values above.
FitsError – If projection names none of the codes above, if fewer than 2 points are given (or fewer than 3 without a pinned crpix), or if the fit is ill-conditioned, for example from collinear points.

Examples

>>> import numpy as np, fitsy
>>> pix = np.array([[100.0, 100.0], [200.0, 100.0], [100.0, 200.0]])
>>> sky = np.array([[10.00, -5.00], [10.05, -5.00], [10.00, -4.95]])
>>> fit = fitsy.fit_wcs(pix, sky, projection="TAN")
>>> fit.rms_arcsec < 1e-6
True