World coordinates

Python

Every fitsy.ImageHdu exposes a wcs() helper that returns a fitsy.Wcs instance (or None if the header carries no WCS). The same parser is reachable from fitsy.FitsFile.wcs(), which additionally resolves -TAB look-up axes.

"""WCS pixel <-> sky transforms on the bundled NGC 2403 image.

Run from the repo root:

    python examples/python/wcs.py
"""

import fitsy
import numpy as np

with fitsy.open("examples/data/ngc2403.fits.gz") as f:
    wcs = f.wcs(0)  # equivalent: f[0].wcs()

# Single pixel -> sky (0-based pixel coordinates).
ra, dec = wcs.pixel_to_celestial(724.0, 1086.0)
print(f"center:     RA={ra:.4f}  Dec={dec:.4f}")

# Sky -> pixel (round-trip).
px, py = wcs.celestial_to_pixel(ra, dec)
print(f"round-trip: ({px:.2f}, {py:.2f})")

# Batch transform: corners + center -> sky.
sky = wcs.pixel_to_celestial_many(
    np.array([[0.0, 0.0], [1447.0, 2171.0], [724.0, 1086.0]])
)
print("corners + center sky:")
print(sky)

Rust

FitsFile::wcs resolves -TAB axes automatically; ImageHdu::wcs is a lighter-weight alternative when you already have the HDU in hand and know there are no tabular axes.

//! WCS pixel-to-sky coordinate transforms using a real FITS image.
//!
//! Run from the repo root with:
//!
//!     cargo run --example wcs

use fitsy::{FitsFile, Hdu, Wcs};

fn main() -> Result<(), fitsy::FitsError> {
    let f = FitsFile::open("examples/data/ngc2403.fits.gz")?;

    // FitsFile::wcs(hdu_index, alt_char) resolves -TAB axes automatically.
    // Use ' ' (space) for the primary WCS; 'A'..'Z' for alternates.
    let wcs: Wcs = f.wcs(0, ' ')?.expect("no WCS in HDU 0");

    // Single pixel -> sky (0-based pixel coordinates).
    // The center of the first pixel is (0.0, 0.0).
    let (ra, dec) = wcs.pixel_to_celestial(724.0, 1086.0)?;
    println!("center:     RA={ra:.4}  Dec={dec:.4}");
    // center:     RA=114.2089  Dec=65.5917

    // Sky -> pixel (round-trip).
    let (px, py) = wcs.celestial_to_pixel(ra, dec)?;
    println!("round-trip: ({px:.2}, {py:.2})");
    // round-trip: (724.00, 1086.00)

    // Batch transform: corners + center -> sky.
    let pairs = vec![(0.0_f64, 0.0_f64), (1447.0, 2171.0), (724.0, 1086.0)];
    let sky = wcs.pixel_to_celestial_many(&pairs)?;
    println!("corners + center:");
    for ((px, py), (ra, dec)) in pairs.iter().zip(&sky) {
        println!("  ({px:.0}, {py:.0}) -> RA={ra:.4}  Dec={dec:.4}");
    }

    // Local pixel scale at the center (arcseconds per pixel, each axis).
    let (sx, sy) = wcs.pixel_scale_at(724.0, 1086.0)?;
    println!("pixel scale: {sx:.4}\" x {sy:.4}\"/px");

    // Full N-axis pixel_to_world / world_to_pixel (useful when the
    // image has non-celestial axes, e.g. spectral).
    let world = wcs.pixel_to_world(&[724.0, 1086.0])?;
    println!("world:  {world:?}");

    // Parsing directly from a Header skips -TAB resolution and is
    // lighter-weight when you know the image has no tabular axes.
    if let Hdu::Image(img) = f.hdu(0)? {
        let _wcs2 = Wcs::from_header(img.header(), ' ')?.expect("no WCS");
    }

    Ok(())
}

Pixel coordinate convention

Important

Both the Python and Rust APIs default to 0-based pixel coordinates (numpy / C convention): the center of the first pixel is 0.0, and the center of the last pixel along NAXISn is float(NAXISn) - 1. A numpy index [row, col] maps directly to pixel_to_celestial(col, row).

Pass origin=1 (Python) to use the FITS 1-based convention (matching CRPIX in the header). The Rust API is always 0-based; subtract 1 from FITS-native coordinates before calling.

# The pixel at numpy index [row, col] = [128, 256]
ra, dec = wcs.pixel_to_celestial(256.0, 128.0)

Batch transforms

Use pixel_to_celestial_many() and celestial_to_pixel_many() for (N, 2) numpy inputs. The Rust equivalents take &[(f64, f64)] slices.

Fitting a WCS

fitsy.fit_wcs() (Python) and fitsy::wcs::fit_celestial_wcs (Rust) solve for a celestial WCS given pixel <-> sky correspondences.

"""Fit a celestial WCS from pixel/sky correspondences.

Run from the repo root:

    python examples/python/fit_wcs.py
"""

import fitsy
import numpy as np

pix = np.array(
    [
        [100.0, 100.0],
        [200.0, 100.0],
        [100.0, 200.0],
        [200.0, 200.0],
    ]
)
sky = np.array(
    [
        [10.00, -5.00],
        [10.05, -5.00],
        [10.00, -4.95],
        [10.05, -4.95],
    ]
)

fit = fitsy.fit_wcs(pix, sky, projection="TAN")
print(f'rms = {fit.rms_arcsec:.3f}"  max = {fit.max_arcsec:.3f}"')

# `fit.wcs` is a fully usable Wcs.
ra, dec = fit.wcs.pixel_to_celestial(150.0, 150.0)
print(f"center: RA={ra:.4f}  Dec={dec:.4f}")

# Serialize back to a header dict for writing.
header = fit.wcs.to_header()
print("CRVAL1:", header["CRVAL1"])

//! Fit a celestial WCS from pixel / sky reference correspondences.
//!
//! Run from the repo root with:
//!
//!     cargo run --example fit_wcs

use fitsy::wcs::{WcsFitOptions, fit_celestial_wcs};

fn main() -> Result<(), fitsy::FitsError> {
    // Four corner correspondences: pixel (0-based) and sky in degrees.
    let pixels = vec![
        (100.0_f64, 100.0),
        (200.0, 100.0),
        (100.0, 200.0),
        (200.0, 200.0),
    ];
    let sky = vec![
        (10.00_f64, -5.00),
        (10.05, -5.00),
        (10.00, -4.95),
        (10.05, -4.95),
    ];

    // Default: TAN projection, CRPIX solved as a free parameter, no SIP.
    let opts = WcsFitOptions::default();
    let fit = fit_celestial_wcs(&pixels, &sky, &opts)?;

    println!(
        "rms = {:.3}\"  max = {:.3}\"",
        fit.rms_arcsec, fit.max_arcsec
    );
    for (i, (dx, dy)) in fit.residuals_arcsec.iter().enumerate() {
        println!("  point {i}: ({dx:+.3}\", {dy:+.3}\")");
    }

    // The fitted Wcs is fully usable for forward / inverse transforms.
    let (ra, dec) = fit.wcs.pixel_to_celestial(150.0, 150.0)?;
    println!("center: RA={ra:.4}  Dec={dec:.4}");

    Ok(())
}