Using astropy.fits and numpy to apply coincidence corrections to SWIFT fits image

Question

This question may be a little specialist, but hopefully someone might be able to help. I normally use IDL, but for developing a pipeline I'm looking to use python to improve running times.

My fits file handling setup is as follows:

import numpy as numpy
from astropy.io import fits

#Directory: /Users/UCL_Astronomy/Documents/UCL/PHASG199/M33_UVOT_sum/UVOTIMSUM/M33_sum_epoch1_um2_norm.img
with fits.open('...') as ima_norm_um2:
    #Open UVOTIMSUM file once and close it after extracting the relevant values:
    ima_norm_um2_hdr  = ima_norm_um2[0].header
    ima_norm_um2_data = ima_norm_um2[0].data
    #Individual dimensions for number of x pixels and number of y pixels:
    nxpix_um2_ext1 = ima_norm_um2_hdr['NAXIS1']
    nypix_um2_ext1 = ima_norm_um2_hdr['NAXIS2']
    #Compute the size of the images (you can also do this manually rather than calling these keywords from the header):
    #Call the header and data from the UVOTIMSUM file with the relevant keyword extensions:
    corrfact_um2_ext1 = numpy.zeros((ima_norm_um2_hdr['NAXIS2'], ima_norm_um2_hdr['NAXIS1']))
    coincorr_um2_ext1 = numpy.zeros((ima_norm_um2_hdr['NAXIS2'], ima_norm_um2_hdr['NAXIS1']))

#Check that the dimensions are all the same:
print(corrfact_um2_ext1.shape)
print(coincorr_um2_ext1.shape)
print(ima_norm_um2_data.shape)

# Make a new image file to save the correction factors:
hdu_corrfact = fits.PrimaryHDU(corrfact_um2_ext1, header=ima_norm_um2_hdr)
fits.HDUList([hdu_corrfact]).writeto('.../M33_sum_epoch1_um2_corrfact.img')

# Make a new image file to save the corrected image to:
hdu_coincorr = fits.PrimaryHDU(coincorr_um2_ext1, header=ima_norm_um2_hdr)
fits.HDUList([hdu_coincorr]).writeto('.../M33_sum_epoch1_um2_coincorr.img')

I'm looking to then apply the following corrections:

    # Define the variables from Poole et al. (2008) "Photometric calibration of the Swift ultraviolet/optical telescope":

alpha =  0.9842000
ft    =  0.0110329
a1    =  0.0658568
a2    = -0.0907142
a3    =  0.0285951
a4    =  0.0308063

for i in range(nxpix_um2_ext1 - 1): #do begin
    for j in range(nypix_um2_ext1 - 1): #do begin
        if (numpy.less_equal(i, 4) | numpy.greater_equal(i, nxpix_um2_ext1-4) | numpy.less_equal(j, 4) | numpy.greater_equal(j, nxpix_um2_ext1-4)): #then begin
            #UVM2
            corrfact_um2_ext1[i,j] == 0
            coincorr_um2_ext1[i,j] == 0
        else:
            xpixmin = i-4
            xpixmax = i+4
            ypixmin = j-4
            ypixmax = j+4
            #UVM2
            ima_UVM2sum = total(ima_norm_um2[xpixmin:xpixmax,ypixmin:ypixmax])
            xvec_UVM2 = ft*ima_UVM2sum
            fxvec_UVM2 = 1 + (a1*xvec_UVM2) + (a2*xvec_UVM2*xvec_UVM2) + (a3*xvec_UVM2*xvec_UVM2*xvec_UVM2) + (a4*xvec_UVM2*xvec_UVM2*xvec_UVM2*xvec_UVM2)
            Ctheory_UVM2 = - alog(1-(alpha*ima_UVM2sum*ft))/(alpha*ft)
            corrfact_um2_ext1[i,j] = Ctheory_UVM2*(fxvec_UVM2/ima_UVM2sum)
            coincorr_um2_ext1[i,j] = corrfact_um2_ext1[i,j]*ima_sk_um2[i,j]

The above snippet is where it is messing up, as I have a mixture of IDL syntax and python syntax. I'm just not sure how to convert certain aspects of IDL to python. For example, the ima_UVM2sum = total(ima_norm_um2[xpixmin:xpixmax,ypixmin:ypixmax]) I'm not quite sure how to handle.

I'm also missing the part where it will update the correction factor and coincidence correction image files, I would say. If anyone could have the patience to go over it with a fine tooth comb and suggest the neccessary changes I need that would be excellent.

The original normalised image can be downloaded here: Replace ... in above code with this file

MSeifert · Accepted Answer · 2016-02-15T15:52:08.667

2

One very important thing about numpy is that it does every mathematical or comparison function on an element-basis. So you probably don't need to loop through the arrays.

So maybe start where you convolve your image with a sum-filter. This can be done for 2D images by astropy.convolution.convolve or scipy.ndimage.filters.uniform_filter

I'm not sure what you want but I think you want a 9x9 sum-filter that would be realized by

from scipy.ndimage.filters import uniform_filter
ima_UVM2sum = uniform_filter(ima_norm_um2_data, size=9)

since you want to discard any pixel that are at the borders (4 pixel) you can simply slice them away:

ima_UVM2sum_valid = ima_UVM2sum[4:-4,4:-4]

This ignores the first and last 4 rows and the first and last 4 columns (last is realized by making the stop value negative)

now you want to calculate the corrections:

xvec_UVM2 = ft*ima_UVM2sum_valid
fxvec_UVM2 = 1 + (a1*xvec_UVM2) + (a2*xvec_UVM2**2) + (a3*xvec_UVM2**3) + (a4*xvec_UVM2**4)
Ctheory_UVM2 = - np.alog(1-(alpha*ima_UVM2sum_valid*ft))/(alpha*ft)

these are all arrays so you still do not need to loop.

But then you want to fill your two images. Be careful because the correction is smaller (we inored the first and last rows/columns) so you have to take the same region in the correction images:

corrfact_um2_ext1[4:-4,4:-4] = Ctheory_UVM2*(fxvec_UVM2/ima_UVM2sum_valid)
coincorr_um2_ext1[4:-4,4:-4]  = corrfact_um2_ext1[4:-4,4:-4] *ima_sk_um2

still no loop just using numpys mathematical functions. This means it is much faster (MUCH FASTER!) and does the same.

Maybe I have forgotten some slicing and that would yield a Not broadcastable error if so please report back.

Just a note about your loop: Python's first axis is the second axis in FITS and the second axis is the first FITS axis. So if you need to loop over the axis bear that in mind so you don't end up with IndexErrors or unexpected results.

edited Feb 15 '16 at 15:52

answered Feb 15 '16 at 15:34

MSeifert

145,886
38
333
352

Hi Michael. Many thanks for your lengthy reply. Very much appreciated. Just a quick query, should the line xvec_UVM2 = ft*ima_UVM2sum actually be xvec_UVM2 = ft*ima_UVM2sum_valid? (P.S. Will work this into my script and see how it goes.) – GCien Feb 15 '16 at 15:37
1

Yes, of course! I forgot to slice the first and last rows and columns when I wrote the answer and only included it later. So every subsequent ``ima_UVM2sum`` should be ``ima_UVM2sum_valid``. – MSeifert Feb 15 '16 at 15:39
Ideal. I'm getting 'module' object has no attribute 'ndimage', even though I have imported the scipy module... – GCien Feb 15 '16 at 15:42
1

try ``from scipy.ndimage.filters import uniform_filter`` and then replace ``scipy.ndimage.filters.uniform_filter`` by ``uniform_filter`` – MSeifert Feb 15 '16 at 15:43
Seems to have worked. And is there a way of clobbering the written fits file, such as fits.HDUList([hdu_coincorr]).writeto('...', clobber=Y)? – GCien Feb 15 '16 at 15:47
I've accepted the answer. So now, would you know how I go about updating the newly made correction factor and coincidence loss files with the newly calculated values? Is that simply using the astropy fits.update command? – GCien Feb 15 '16 at 15:51
1

You already wrote these correction images to disk (during ``writeto``), where do you want them to update? Or am I missing a step? – MSeifert Feb 15 '16 at 15:56
I have both image files as just arrays of 0s with the same dimesions as the loaded file, so I'm thinking should I add the 4 write.to lines in the top section of the script (in my question) to the very end, after all of these calculations have been made? (Sorry if these are really stupid questions!) – GCien Feb 15 '16 at 16:00
1

Ah sorry, yes, you have to move all the ``writeto`` related code after the calculation. :-) – MSeifert Feb 15 '16 at 16:02
I now have the same dimension image, but all with NaNs values...should I move this to a chat? – GCien Feb 15 '16 at 16:04
Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/103502/discussion-between-mseifert-and-michael-roberts). – MSeifert Feb 15 '16 at 16:05

Using astropy.fits and numpy to apply coincidence corrections to SWIFT fits image

1 Answers1