How does one update a FITS data array with PyFITS?

Question

I'm trying to create a new FITS file out of two older ones using PyFITS.

import pyfits
from sys import stdout
from sys import argv
import time

file1 = argv[1]
file2 = argv[2]

hdu1 = pyfits.open(file1)
hdu2 = pyfits.open(file2)
new0 = hdu1[0]
new1 = hdu1[0]

sci1 = hdu1[0].data
sci2 = hdu2[0].data

for r in range(0, len(sci1)):
    for c in range(0, len(sci1[r])):
       add = sci1[r][c] + sci2[r][c]
       new0.data[r][c] = add
for r in range(0, len(sci1)):
    for c in range(0, len(sci1[r])):
       print "(" + str(r) + ", " + str(c) + ") FirstVal = " + str(sci1[r][c]) + " || SecondVal = " + str(sci2[r][c])
       print "\t New File/Add = " + str(new0.data[r][c])

All it prints out is the first value, i.e. sci1[r][c]. This means that the variable isn't being modified at all. How can I make it modify? I'm very new to using FITS.

Answer 1

What you have done here is make sci1 a reference to new0.data which means the assignment to new0 also changes sci1, so it is modifying the intended variable but your print loop is printing the same object twice. If you want to have a copy instead of reference you have to use the objects copy method, in this case sci0 = new0.data.copy()

This is also not the way you are supposed to use numpy which pyfits uses to represent its images. Instead of loops you apply operations to full arrays which is in most cases easier to read and significantly faster. If you want to add two fits images represented as numpy arrays inplace:

new0.data += new1.data
print new0.data

or if you want to create a new image out of the sum of both inputs:

sum_image = new0.data + new1.data
# put it into an pyfits HDU (primary fits extension)
hdu = pyfits.PrimaryHDU(data=sum_image)