2

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?

Jim421616
  • 1,087
  • 1
  • 10
  • 29
  • Is your primary concern to **save the mask** or to use another program (unrelated to python/matplotlib) to display the FITS files with mask (e.g. IRAF, MIDAS or DS9)? – MSeifert Sep 19 '17 at 23:47
  • My next step is to fit a mathematical model to the un-masked portion of the image, so I guess I don't really need to save the mask. On the other hand, I'm doing it as a project for my MSc and I'm hoping to be able to show my supervisor my intermediate images in DS9. – Jim421616 Sep 19 '17 at 23:49
  • But for the purpose of this question: You want DS9 to display the image without masked pixels? – MSeifert Sep 19 '17 at 23:51
  • Oh, then the answer is yes :) – Jim421616 Sep 19 '17 at 23:52

1 Answers1

2

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)
MSeifert
  • 118,681
  • 27
  • 271
  • 293
  • Ok, thanks, MSeifert. Trying to install CCDData in Anaconda results in PackageNotFoundError. There is a package called ccdproc, but that doesn't seem to have the same functionality as CCDData. Has the package been replaced by something? – Jim421616 Sep 19 '17 at 23:17
  • @Jim421616 It's part of `astropy`, like the `astropy.io.fits` module you're using. If you have a recent version of astropy it should be importable using `from astropy.nddata import CCDData`. If you don't have a recent astropy version you can also update astropy with conda using `conda install astropy=2`. – MSeifert Sep 19 '17 at 23:20
  • @Jim421616 Historically it has been part of the `ccdproc` package so if you don't want (or can) upgrade you could use `conda install -c conda-forge ccdproc` and use `from ccdproc import CCDData` instead. – MSeifert Sep 19 '17 at 23:22
  • Those are very helpful comments; I've edited my post to include your ideas. color_masked is still showing the pixels I'm trying to mask though. I've tried importing CCDData from both astropy.nddata and from ccdproc (after following your installation instructions). Am I doing something wrong in using ccddata? – Jim421616 Sep 19 '17 at 23:40
  • 1
    @Jim421616 Thanks for the clarification. It seems like I originally misunderstood your question. I updated the answer. It's a bit "ugly" but based on my DS9 viewer it works. The masked values still show up as bright white spots but they don't affect the color-scale anymore which gives much more contrast. – MSeifert Sep 20 '17 at 00:05
  • Thank you! I guess I can use the mask again when I fit the model. – Jim421616 Sep 20 '17 at 00:10