Source code for mufasa.utils.fits_utils
"""
Utilities for working with .fits files, including functions adapted from
FITS_tools to reduce dependencies on the aging package.
"""
import numpy as np
import astropy.wcs as pywcs
from astropy import coordinates
from astropy import units as u
[docs]
def get_pixel_mapping(header1, header2):
"""
Compute the pixel mapping between two FITS headers or WCS objects.
This function determines how pixel coordinates in one FITS header
or WCS (`header1`) map to pixel coordinates in another (`header2`).
It returns a grid that describes the transformation in terms of
world coordinates.
Parameters
----------
header1 : `~astropy.io.fits.Header` or `~astropy.wcs.WCS`
The header or WCS corresponding to the input image whose pixel
coordinates are being transformed.
header2 : `~astropy.io.fits.Header` or `~astropy.wcs.WCS`
The header or WCS corresponding to the output image onto which
the input image will be interpolated.
Returns
-------
grid : `~numpy.ndarray`
A 2D array of shape `(2, n, m)` where `n` and `m` are the dimensions
of `header2`. The array contains the mapped pixel coordinates from
`header1` in the following order:
- `grid[0, :, :]`: The y-coordinates in `header1`.
- `grid[1, :, :]`: The x-coordinates in `header1`.
Raises
------
TypeError
If either `header1` or `header2` is not an instance of
`astropy.io.fits.Header` or `astropy.wcs.WCS`.
NotImplementedError
If the coordinate system type (`CTYPE`) in the headers is not
recognized or if unit conversions between coordinate systems
are not implemented.
Notes
-----
- The function supports transformations between images that share
compatible celestial coordinate systems, such as `GLON`/`GLAT`
and `RA`/`DEC`.
- If the coordinate systems differ, they will be transformed to
align before computing the pixel mapping.
- The code here is borrowed from `FITS_tools.downsample`.
Examples
--------
>>> from astropy.io import fits
>>> from astropy.wcs import WCS
>>> header1 = fits.Header() # Define input header or WCS
>>> header2 = fits.Header() # Define output header or WCS
>>> grid = get_pixel_mapping(header1, header2)
>>> print(grid.shape)
(2, n, m) # Example output shape corresponding to header2 dimensions
"""
wcs1 = _load_wcs_from_header(header1)
wcs2 = _load_wcs_from_header(header2)
if not all([w1 == w2 for w1, w2 in zip(wcs1.wcs.ctype, wcs2.wcs.ctype)]):
allowed_coords = ('GLON', 'GLAT', 'RA', 'DEC')
if all([(any(word in w1 for word in allowed_coords) and
any(word in w2 for word in allowed_coords))
for w1, w2 in zip(wcs1.wcs.ctype, wcs2.wcs.ctype)]):
csys1 = _ctype_to_csys(wcs1.wcs)
csys2 = _ctype_to_csys(wcs2.wcs)
convert_coordinates = True
else:
raise NotImplementedError(
"Unit conversions between {0} and {1} have not yet been implemented."
.format(wcs1.wcs.ctype, wcs2.wcs.ctype)
)
else:
convert_coordinates = False
outshape = [wcs2.naxis2, wcs2.naxis1]
yy2, xx2 = np.indices(outshape)
lon2, lat2 = wcs2.wcs_pix2world(xx2, yy2, 0)
if convert_coordinates:
C2 = coordinates.SkyCoord(lon2, lat2, unit=(u.deg, u.deg), frame=csys2)
C1 = C2.transform_to(csys1)
lon2, lat2 = C1.spherical.lon.deg, C1.spherical.lat.deg
xx1, yy1 = wcs1.wcs_world2pix(lon2, lat2, 0)
grid = np.array([yy1.reshape(outshape), xx1.reshape(outshape)])
return grid
def _load_wcs_from_header(header):
"""
Load a WCS object from a FITS header or verify an existing WCS instance.
This function ensures that a `pywcs.WCS` object is created from the provided
FITS header or verifies that the input is already a valid WCS instance.
It also adds `naxis1` and `naxis2` attributes to the WCS object if they
are not already defined.
Parameters
----------
header : `~astropy.io.fits.Header` or `~pywcs.WCS`
The FITS header or WCS object to process.
Returns
-------
wcs : `~pywcs.WCS`
A WCS object with required attributes set.
Raises
------
TypeError
If the input is neither a valid FITS header nor a `pywcs.WCS` instance.
Notes
-----
- The code here is borrowed from `FITS_tools.downsample`.
"""
if issubclass(pywcs.WCS, header.__class__):
wcs = header
else:
try:
wcs = pywcs.WCS(header)
except:
raise TypeError("header must either be a astropy.io.fits.Header "
" or pywcs.WCS instance")
if not hasattr(wcs, 'naxis1'):
wcs.naxis1 = header['NAXIS1']
if not hasattr(wcs, 'naxis2'):
wcs.naxis2 = header['NAXIS2']
return wcs
def _ctype_to_csys(wcs):
"""
Convert WCS `CTYPE` information to a coordinate system name.
This function determines the coordinate system (e.g., `fk5`, `fk4`, or `galactic`)
based on the `CTYPE` attribute and the equinox in the provided WCS object.
Parameters
----------
wcs : `~pywcs.WCS`
A WCS object containing `CTYPE` and `equinox` information.
Returns
-------
csys : str
The name of the coordinate system. Possible values are:
- `'fk5'` for J2000 equinox.
- `'fk4'` for B1950 equinox.
- `'galactic'` for Galactic coordinates.
Raises
------
NotImplementedError
If the equinox is neither 2000 nor 1950 when `CTYPE` indicates
celestial coordinates.
Notes
-----
- This code here is borrowed from `FITS_tools.downsample`.
"""
ctype = wcs.ctype[0]
if 'RA' in ctype or 'DEC' in ctype:
if wcs.equinox == 2000:
return 'fk5'
elif wcs.equinox == 1950:
return 'fk4'
else:
raise NotImplementedError("Non-fk4/fk5 equinoxes are not allowed")
elif 'GLON' in ctype or 'GLAT' in ctype:
return 'galactic'
