How to change the header after resizing a submap?

Question

I load a full map

from astropy.io import fits
from astropy.wcs import wcs

mapheader = fits.getheader(MapFile, 0)
mapdata = fits.getdata(MapFile, 0)
w = wcs.WCS(mapheader)

and I take a squared submap from it assuming the center is in RA,DEC in degrees this can be easily done using CutOut2D

from astropy.nddata import Cutout2D
from astropy import coordinates
from astropy import units as u

center = coordinates.SkyCoord(RA*u.deg, DEC*u.deg, frame='fk5')
submap = Cutout2D(mapdata, center, size=16*u.pix, wcs=w)
submapheader = submap.wcs.to_header()

The relevant difference in headers is that it moves the reference pixel “CRPIX”

If I resize the image, for example, I do an interpolation and pass the image from 16 pixels to 128 pixels

from scipy.misc import imresize
newsubmap = imresize(submap.data, (128,128), interp=‘cubic’)

how should I change the header to get a good projection on newsubmap?

I tried multiplying the reference pixel by the resizing factor, which is 128 in this example, but it's not simple as that

score 1 · Answer 1 · answered Nov 08 '16 at 19:12

scipy.misc.imresize in your case resizes the image to (128, 128). Given your statement:

I tried multiplying the reference pixel by the resizing factor, which is 128 in this example, but it's not simple as that.

I assume that is the first pitfall here. Are you really sure you want to resize to (128, 128) or do you want to "magnify" it by a factor of 128 or even by a factor of 1.28 (note that imresize uses fractional values!)

>>> from scipy.misc import imresize
>>> import numpy as np
>>> imresize(np.ones((1000, 1000)), (100, 100)).shape  # to (100, 100)
(100, 100)
>>> imresize(np.ones((1000, 1000)), 50).shape  # half the size. 50->50%
(500, 500)

Please make sure you're using imresize correctly.

So the next step is to reshape the WCS. Fortunatly WCS allows slicing, so this is quite easy if you know the original shape and the resizing factor.

Say you have a WCS like this:

>>> im.wcs
WCS Keywords
Number of WCS axes: 2
CTYPE : 'PIXEL'  'PIXEL'  
CRVAL : 2044.203  239.489  
CRPIX : 1022.1  119.7  
PC1_1 PC1_2  : 2.0  0.0  
PC2_1 PC2_2  : 0.0  2.0  
CDELT : 1.0  1.0

you can slice it given a step:

>>> im.wcs[::2, ::2]  # half the size
WCS Keywords
Number of WCS axes: 2
CTYPE : 'PIXEL'  'PIXEL'  
CRVAL : 2044.203  239.489  
CRPIX : 511.30000000000001  60.100000000000001  
PC1_1 PC1_2  : 2.0  0.0  
PC2_1 PC2_2  : 0.0  2.0  
CDELT : 2.0  2.0

Or slice it with steps smaller than 1 to increase it:

>>> im.wcs[::1/128, ::1/128]  # assuming you increase each axis by a factor of 128.
WCS Keywords
Number of WCS axes: 2
CTYPE : 'PIXEL'  'PIXEL'  
CRVAL : 2044.203  239.489  
CRPIX : 130765.3  15258.1  
PC1_1 PC1_2  : 2.0  0.0  
PC2_1 PC2_2  : 0.0  2.0  
CDELT : 0.0078125  0.0078125

Note that PC, CD and possible also SIP and other distortions are ignored then. You have to manually process them. However CDELT values will be processed, so simple FITS files will be processed correctly.

Note: I removed the NAXIS keywords, because they are subject to change for the next version so these wouldn't be reliable anyway. These are not processed currently and in the next version they will processed only if "start, stop and step" are integers or None.

How to change the header after resizing a submap?

1 Answers1