Search code examples
pythonmaskastropyfits

Saving masked image as FITS

I've constructed an image from some FITS files, and I want to save the resultant masked image as another FITS file. Here's my code:

import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt
#from astropy.nddata import CCDData
from ccdproc import CCDData

hdulist1 = fits.open('wise_neowise_w1-MJpersr.fits')
hdulist2 = fits.open('wise_neowise_w2-MJpersr.fits')

data1_raw = hdulist1[0].data
data2_raw = hdulist2[0].data

#   Hide negative values in order to take logs
#   Where {condition}==True, return data_raw, else return np.nan
data1 = np.where(data1_raw >= 0, data1_raw, np.nan)
data2 = np.where(data2_raw >= 0, data2_raw, np.nan)

#   Calculation and image subtraction
w1mag = -2.5 * (np.log10(data1) - 9.0)
w2mag = -2.5 * (np.log10(data2) - 9.0)
color = w1mag - w2mag

##   Find upper and lower 5th %ile of pixels
mask_percent = 5
masked_value_lower = np.nanpercentile(color, mask_percent)
masked_value_upper = np.nanpercentile(color, (100 - mask_percent))

##   Mask out the upper and lower 5% of pixels
##   Need to hide values outside the range [lower, upper]
color_masked = np.ma.masked_outside(color, masked_value_lower, masked_value_upper)
color_masked = np.ma.masked_invalid(color_masked)

plt.imshow(color)
plt.title('color')
plt.savefig('color.png', overwrite = True)
plt.imshow(color_masked)
plt.title('color_masked')
plt.savefig('color_masked.png', overwrite = True)

fits.writeto('color.fits',
             color,
             overwrite = True)
ccd = CCDData(color_masked, unit = 'adu')
ccd.write('color_masked.fits', overwrite = True))

hdulist1.close()
hdulist2.close()

When I use matplotlib.pyplot to imshow the images color and color_masked, they look as I expect: enter image description here enter image description here However, my two output files, color_masked.fits == color.fits. I think somehow I'm not quite understanding the masking process properly. Can anyone see where I've gone wrong?

Solution

  • astropy.io.fits only handles normal arrays and that means it just ignores/discards the mask of your MaskedArray.

    Depending on your use-case you have different options:

    Saving the file so other FITS programs recognize the mask

    I actually don't think that's possible. But some programs like DS9 can handle NaNs, so you could just set the masked values to NaN for the purpose of displaying them:

    data_naned = np.where(color_masked.mask, np.nan, color_masked)
fits.writeto(filename, data_naned, overwrite=True)

    They do still show up as "bright white spots" but they don't affect the color-scale.

    If you want to take this a step further you could replace the masked pixels using a convolution filter before writing them to a file. Not sure if there's one in astropy that only replaces masked pixels though.

    Saving the mask as extension so you can read them back

    You could use astropy.nddata.CCDData (available since astropy 2.0) to save it as FITS file with mask:

    from astropy.nddata import CCDData

ccd = CCDData(color_masked, unit='adu')
ccd.write('color_masked.fits', overwrite=True)

    Then the mask will be saved in an extension called 'MASK' and it can be read using CCDData as well:

    ccd2 = CCDData.read('color_masked.fits')

    The CCDData behaves like a masked array in normal NumPy operations but you could also convert it to a masked-array by hand:

    import numpy as np
marr = np.asanyarray(ccd2)