Extracting WCS Coordinates from FITS

Question

I am trying to change the coordinate system of a FITS file from its original Equatorial coordinate system to Galactic coordinates (in degrees)in order to manipulate the resulting FITS image using these coordinates.

For this, I would need to extract an array which contains the Equatorial positions of each pixel in order to transform them into the desired Galactic Coordinates. This is where my knowledge is limited, and cannot seem to figure out how to extract that array.

Ultimately, I want to slice the image based on latitudes in the manner:

import pyfits

f = pyfits.open("im1.fits")
h, data = f[0].header, f[0].data

#Here I would include the transformation from Equatorial to Galactic of
the position array doing something like:
coord = SkyCoord(ra, dec, frame=Galactic, unit='deg')

#This would do the slicing
from astropy.nddata import Cutout2D
from astropy import units as u
from astropy import coordinates

#Given a longitude and latitude
size = [mylon, mylat]
cut = Cutout2D(f, center, size, wcs=wcs)

Answer 1

You most likely want to regrid the data rather than transform all of the positions; if the data are on an equatorial grid, it is not possible to slice them along Galactic coordinates. reproject is the best astropy-based tool for this job. You'll need to set up a Galactic header and reproject to that:

import reproject
galheader = fits.Header.fromtextfile('gal.hdr')
myfitsfile = fits.open('im1.fits')
newim, weights = reproject.reproject_interp(myfitsfile, galheader)

You can also use reproject.reproject_exact, which uses a different reprojection algorithm.