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 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.set_axislabel('Right Ascension', fontsize=16) overlay.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) & (0 <= stars_y) & (stars_y < hdu.data.shape) ) ax.scatter( stars_x[valid_stars], stars_y[valid_stars], s=50, edgecolor='white', facecolor="none", transform=ax.get_transform("pixel"), )