How do I determine if a (RA,DEC) is in the field of view for a FITS file?

Astropy
, ,
1

I’ve got a FITS file from astrometry.net containing an image, and I’d like to filter a table containing (RA, DEC) values to just those that are in the field of view of the original image.

My current solution is to use the quiet flag in WCS.all_world2pix, but the overall approach seems pretty arcane. I was wondering if there is a more standard API to do this.

This seems related to Issue 11446 and PR 11693.

Here’s a complete example:

# Download the new-image.fits from:
# http://nova.astrometry.net/user_images/5631116#annotated

import matplotlib.pyplot as plt
import numpy as np

from astropy.coordinates import SkyCoord, EarthLocation
import astropy.units as u
from astropy.io import fits
from astropy.wcs import WCS



# Here are my coordinates. For my application, I need to filter
# coordinates to ones that lie in the field-of-view given by the 
# ``new-image.fits`` file.

coords = SkyCoord(
    [
        ( 28.04316907,  24.11731175), 
        ( 84.65821979,  -6.57395768),  # This is the only one in the field of view
        (314.92876576,  38.82292141),
        (229.14911128,   9.7130066 ), 
        (193.4991856 , -60.33541696),
    ], 
    frame="icrs", 
        
    unit=u.deg,
)



with fits.open("new-image.fits") as f:
    hdu = f[0]
    wcs = WCS(hdu.header, naxis=2)
    fig = plt.figure(figsize=(12, 12 * 2/3))
    ax = fig.add_subplot(111, projection=wcs)
    ax.imshow(np.moveaxis(hdu.data, 0, -1), origin="lower")

    overlay = ax.get_coords_overlay('icrs')
    overlay.grid(color='white', ls='dotted', lw=2)
    overlay[0].set_axislabel('Right Ascension', fontsize=16)
    overlay[1].set_axislabel('Declination', fontsize=16)

    # Hack: Use wcs.all_world2pix to filter...
    stars_x, stars_y = wcs.all_world2pix(
        coords.ra, 
        coords.dec, 
        0, 
        quiet=True,
    )

    # Is this really the only way?
    valid_stars = (
        (0 <= stars_x) & (stars_x < hdu.data.shape[2]) 
        & (0 <= stars_y) & (stars_y < hdu.data.shape[1])
    )

    ax.scatter(
        stars_x[valid_stars], 
        stars_y[valid_stars],
        s=50,
        edgecolor='white',
        facecolor="none",
        transform=ax.get_transform("pixel"),
    )
2

I’m not sure if there is a better way to do this, but for what it’s worth, the example you shared is the way I would do it!

By the way: do you have a 3D image cube, or should the valid_stars line be:

valid_stars = (
    (0 <= stars_x) & (stars_x < hdu.data.shape[1]) 
    & (0 <= stars_y) & (stars_y < hdu.data.shape[0])
)

?