How to crop .fits image and save world coordinates for building in Python astropias?

Question

This question attracted me for a while. I am trying to process some large amounts of data that are in the form of a .fits file (of the order of 11000x9000 pixels). What I need to do is create an “enlarged” RA / Dec coordinate graph (ideally using astropy.wcs) for many objects in the sky with outlines from one fits and greyscale (or heatmap from another) file.

My problem is that whenever I cut data from an image (into a region of interest), I lose touch with the coordinates of the sky. This means that the sliced ​​image is not in the right place.

I adapted the example from astrophysical documents to save you from the pain of my data. (Note: I want the outlines to cover more area than the image, no matter which solution for this should work for both data)

Here is the code I'm having problems with:

from matplotlib import pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
from astropy.utils.data import download_file
import numpy as np

fits_file = 'http://data.astropy.org/tutorials/FITS-images/HorseHead.fits'
image_file = download_file(fits_file, cache=True)
hdu = fits.open(image_file)[0]
wmap = WCS(hdu.header)
data = hdu.data

fig = plt.figure()
ax1 = fig.add_subplot(121, projection=wmap)
ax2 = fig.add_subplot(122, projection=wmap)
# Scale input image
bottom, top = 0., 12000.
data = (((top - bottom) * (data - data.min())) / (data.max() - data.min())) + bottom


'''First plot'''
ax1.imshow(data, origin='lower', cmap='gist_heat_r')

# Now plot contours
xcont = np.arange(np.size(data, axis=1))
ycont = np.arange(np.size(data, axis=0))
colors = ['forestgreen','green', 'limegreen']
levels = [2000., 7000., 11800.]

ax1.contour(xcont, ycont, data, colors=colors, levels=levels, linewidths=0.5, smooth=16)
ax1.set_xlabel('RA')
ax1.set_ylabel('Dec')
ax1.set_title('Full image')

''' Second plot ''' 
datacut = data[250:650, 250:650]
ax2.imshow(datacut, origin='lower', cmap=cmap)
ax2.contour(xcont, ycont, data, colors=colors, levels=levels, linewidths=0.5, smooth=16)

ax2.set_xlabel('RA')
ax2.set_ylabel('')
ax2.set_title('Sliced image')
plt.show()

I tried using the WCS cords of my chopped piece to fix this, but I'm not sure if I can transfer it anywhere!

pixcoords = wcs.wcs_pix2world(zip(*[range(250,650),range(250,650)]),1)

Answer 1

The good news: you can just chop your own astropy.WCS, which makes your task relationally trivial:

...

wmapcut = wmap[250:650, 250:650] # sliced here
datacut = data[250:650, 250:650]
ax2 = fig.add_subplot(122, projection=wmapcut) # use sliced wcs as projection
ax2.imshow(datacut, origin='lower', cmap='gist_heat_r')
# contour has to be sliced as well
ax2.contour(np.arange(datacut.shape[0]), np.arange(datacut.shape[1]), datacut, 
            colors=colors, levels=levels, linewidths=0.5, smooth=16)
...

WCS, (., , reproject)