0

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')
Iguananaut
  • 21,810
  • 5
  • 50
  • 63
Wolfgang
  • 329
  • 3
  • 13

1 Answers1

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.

keflavich
  • 18,278
  • 20
  • 86
  • 118
  • This is an excellent approach, thanks for your response. Thank you. However, I would like to point out that it shows an error saying: "ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()". To which all seems very confusing to me. – Wolfgang Jun 16 '15 at 08:07
  • That error occurs if you try to evaluate an array of booleans as a single boolean. For example, if you try `if np.ones(3):`, it will raise that exception. There was a mistake in my code: I wrote `&&` where I should have written `&`. Could you explain where you encountered that error? – keflavich Jun 16 '15 at 12:06
  • 1
    And to be clear, this is all just basic points about Numpy--beyond that fact that the array was stored in a FITS file none of it has to do with FITS or Astropy. I'd recommend finding a good guide on Numpy if you're going to work with it extensively. – Iguananaut Jun 16 '15 at 12:21
  • @keflavich: The error now states that newImData is not defined. I tried a few things, though I am still getting used the very particular notation in python, so I have not been very successful. – Wolfgang Jun 16 '15 at 23:50
  • @Iguananaut: This is very true. However, not having a terribly lengthy knowledge-base on coding using python has made it very difficult to come across some good sources that would pertain to this subject (especially since I wont be working with Numpy in the long run. – Wolfgang Jun 16 '15 at 23:55