I'm having a problem wit fits file manipulation in the astropy package, and I'm in need of some help.
I essentially want to take an image I have in fits file format, and create a new file I need to start inputing correction factors to and a new image which can then be used with the correction factors and the original image to produce a correction image. Each of these will have the same dimensions.
Starting with this:
from astropy.io import fits
# Compute the size of the images (you can also do this manually rather than calling these keywords from the header):
#URL: /Users/UCL_Astronomy/Documents/UCL/PHASG199/M33_UVOT_sum/UVOTIMSUM/M33_sum_epoch1_um2_norm.img
nxpix_um2_ext1 = fits.open('...')[1]['NAXIS1']
nypix_um2_ext1 = fits.open('...')[1]['NAXIS2']
#nxpix_um2_ext1 = 4071 #hima_sk_um2[1].header['NAXIS1'] # IDL: nxpix_uw1_ext1 = sxpar(hima_sk_uw1_ext1,'NAXIS1')
#nypix_um2_ext1 = 4321 #hima_sk_um2[1].header['NAXIS2'] # IDL: nypix_uw1_ext1 = sxpar(hima_sk_uw1_ext1,'NAXIS2')
# Make a new image file with the same dimensions (and headers, etc) to save the correction factors:
coicorr_um2_ext1 = ??[nxpix_um2_ext1,nypix_um2_ext1]
# Make a new image file with the same dimensions (and headers, etc) to save the corrected image:
ima_sk_coicorr_um2_ext1 = ??[nxpix_um2_ext1,nypix_um2_ext1]
Can anyone give me the obvious knowledge I am missing to do this...the last two lines are just there to outline what is missing. I have included ?? to perhaps signal I need something else there perhaps fits.writeto() or something similar...
I think @VincePs answer is correct but I'll add some more information because I think you are not using the capabilities of astropy well here.
First of all Python is zero-based so the primary extension has the number 0. Maybe you got that wrong, maybe you don't - but it's uncommon to access the second HDU so I thought I better mention it.
hdu_num = 0 # Or use = 1 if you really want the second hdu.
First you do not need to open the same file twice, you can open it once and close it after extracting the relevant values:
with fits.open('...') as hdus:
nxpix_um2_ext1 = hdus[hdu_num]['NAXIS1']
nxpix_um2_ext1 = hdus[hdu_num]['NAXIS2']
# Continue without indentation and the file will be closed again.
or if you want to keep the whole header (for saving it later) and the data you can use:
with fits.open('...') as hdus:
hdr = hdus[hdu_num].header
data = hdus[hdu_num].data # I'll also take the data for comparison.
I'll continue with the second approach because I think it's a lot cleaner and you'll have all the data and header values ready.
new_data = np.zeros((hdr['NAXIS2'], hdr['NAXIS1']))
Please note that Python interprets the axis different than IRAF (and I think IDL, but I'm not sure) so you need axis2 as first and axis1 as second element.
So do a quick check that the shapes are the same:
print(new_data.shape)
print(data.shape)
If they are not equal I got confused about the axis in Python (again) but I don't think so. But instead of creating a new array based on the header values you can also create a new array by just using the old shape:
new_data_2 = np.zeros(data.shape)
That will ensure the dimensions and shape is identical. Now you have an empty image. If you rather like a copy then you can, but do not need to, explicitly copy the data (except if you opened the file explicitly in write/append/update mode then you should always copy it but that's not the default.)
new_data = data # or = data.copy() for explicitly copying.
Do your operations on it and if you want to save it again you can use what @VinceP suggested:
hdu = fits.PrimaryHDU(new_data, header=hdr) # This ensures the same header is written to the new file
hdulist = fits.HDUList([hdu])
hdulist.writeto('new.fits')
Please note that you don't have to alter the shape-related header keywords even if you changed the data's shape because during writeto astropy will update these (by default)