FITS Image Array Pixel out of bounds

Question

I am attempting to set a threshold for pixels in a FITS file using the code below. However, I get an error that says that:

IndexError: Index 2620 is out of bounds for axis 0 with size 2620

Any ideas on how to fix this?

This is the code:

from astropy.io import fits
import numpy as np
hdulist = fits.open("12ratio.fits")
origImData = hdulist[0].data
newImData = origImData*0

for x, y in np.nditer(origImData.shape):
        curPixel = origImData[x, y]
        if curPixel > 0.28 or curPixel < 3.11:
                newImData[x, y] = curPixel
        else:
                newImData[x, y] = 0
newhdu = fits.PrimaryHDU(newImData)
newhdulist = fits.HDUList([newhdu])
newhdulist.writeto('modifiedratio12.fits')

Answer 1

The problem is your iteration scheme. You're iterating over the shape of the array, not the pixel indices.

In general there is a better approach for masking:

mask = (origImData > 0.28) & (origImData < 3.11)
newImData[mask] = origImData[mask]

will accomplish what you seem to be aiming for.

Or just as effective:

from astropy.io import fits
import numpy as np
hdulist = fits.open("12ratio.fits")
origImData = hdulist[0].data
mask = (origImData > 0.28) & (origImData < 3.11)
origImData[~mask] = 0

newhdu = fits.PrimaryHDU(data=origImData, header=hdulist[0].header)
newhdu.writeto('modifiedratio12.fits')

Here I've used the numpy array boolean not operation (~) to invert the mask. I've also skipped making an HDUList since that is not needed for a single-extension FITS file.

There are two errors in the code you've shown:

for x, y in np.nditer(origImData.shape):

will assign x and y the values of origImData.shape, i.e. x=2640 and y=(something else), it will not assign them the values 0,1,...,2640 as you intend.

The comparison operation:

    if curPixel > 0.28 or curPixel < 3.11:

will always be true for real numbers: all numbers are either greater than 0.28 or less than 3.11.