Search code examples
pythonastropyastronomyfits

Combine two FITS files by taking a 'pizza wedge' slice?

I would like to combine two FITS files, by taking a slice out of one and inserting it into the other. The slice would be based on an angle measured from the centre pixel, see example image below:

enter image description here

Can this be done using Astropy? There are many questions on combining FITS on the site, but most of these are related to simply adding two files together, rather that combining segments like this.

Solution

  • Here is one recommended approach:

    1.Read in your two files

    Assuming the data is in an ImageHDU data array..

    from astropy.io import fits

# read the numpy arrays out of the files 
# assuming they contain ImageHDUs
data1 = fits.getdata('file1.fits')
data2 = fits.getdata('file2.fits')

    2. Cut out the sections and put them into a new numpy array

    Build up indices1 & indices2 for the desired sections in the new file... A simple numpy index to fill in the missing section into a new numpy array.

    After being inspired by https://stackoverflow.com/a/18354475/15531842

    The sector_mask function defined in the answer to get indices for each array using angular slices.

    mask = sector_mask(data1.shape, centre=(53,38), radius=100, angle_range=(280,340))
indices1 = ~mask
indices2 = mask

    Then these indices can be used to transfer the data into a new array.

    import numpy as np
newdata = np.zeros_like(data1)
newdata[indices1] = data1[indices1]
newdata[indices2] = data2[indices2]

    If the coordinate system is well known then it may be possible to use astropy's Cutout2D class, although I was not able to figure out how to fully use it. It wasn't clear if it could do an angular slice from the example given. See astropy example https://docs.astropy.org/en/stable/nddata/utils.html#d-cutout-using-an-angular-size

    3a. Then write out the new array as a new file.

    If special header information is not needed in the new file. Then the numpy array with the new image can be written out to a FITS file with one line of astropy code.

    # this is an easy way to write a numpy array to FITS file
# no header information is carried over
fits.writeto('file_combined.fits', data=newdata)

    3b. Carry the FITS header information over to the new file

    If there is a desire to carry over header information then an ImageHDU can be built from the numpy array and include the desired header as a dictionary.

    img_hdu = fits.ImageHDU(data=newdata, header=my_header_dict)

hdu_list = fits.HDUList()
hdu_list.append(fits.PrimaryHDU())
hdu_list.append(img_hdu)

hdu_list.writeto('file_combined.fits')