Python astropy.io.fits.PrimaryHDU() Examples

def setup_class(self) -> None:

        self.limit = 1e-10
        self.test_dir = os.path.dirname(__file__) + '/'

        np.random.seed(1)
        images = np.random.normal(loc=0, scale=2e-4, size=(5, 11, 11))

        hdu = fits.PrimaryHDU()
        header = hdu.header
        header['INSTRUME'] = 'IMAGER'
        header['HIERARCH ESO DET EXP NO'] = 1
        header['HIERARCH ESO DET NDIT'] = 5
        header['HIERARCH ESO INS PIXSCALE'] = 0.01
        header['HIERARCH ESO ADA POSANG'] = 10.
        header['HIERARCH ESO ADA POSANG END'] = 20.
        header['HIERARCH ESO SEQ CUMOFFSETX'] = 5.
        header['HIERARCH ESO SEQ CUMOFFSETY'] = 5.
        hdu.data = images
        hdu.writeto(self.test_dir+'images.fits') 

def from_tree(cls, data, ctx):
        hdus = []
        first = True
        for hdu_entry in data:
            header = fits.Header([fits.Card(*x) for x in hdu_entry['header']])
            data = hdu_entry.get('data')
            if data is not None:
                try:
                    data = data.__array__()
                except ValueError:
                    data = None
            if first:
                hdu = fits.PrimaryHDU(data=data, header=header)
                first = False
            elif data.dtype.names is not None:
                hdu = fits.BinTableHDU(data=data, header=header)
            else:
                hdu = fits.ImageHDU(data=data, header=header)
            hdus.append(hdu)
        hdulist = fits.HDUList(hdus)
        return hdulist 

def test_insert_groupshdu_to_non_empty_list(self):
        """Tests inserting a Simple GroupsHDU to an empty HDUList."""

        hdul = fits.HDUList()
        hdu = fits.PrimaryHDU(np.arange(100, dtype=np.int32))
        hdul.insert(0, hdu)
        hdu = fits.GroupsHDU()

        with pytest.raises(ValueError):
            hdul.insert(1, hdu)

        info = [(0, 'PRIMARY', 1, 'GroupsHDU', 8, (), '',
                 '1 Groups  0 Parameters'),
                (1, '', 1, 'ImageHDU', 6, (100,), 'int32', '')]

        hdul.insert(0, hdu)

        assert hdul.info(output=False) == info

        hdul.writeto(self.temp('test-insert.fits'))

        assert fits.info(self.temp('test-insert.fits'), output=False) == info 

def test_is_fits_image_file_table_img(self, dp_mock):
        with tempfile.TemporaryDirectory() as tmpdir:
            with self.settings(MEDIA_ROOT=tmpdir):
                img_file = os.path.join(tmpdir, 'img.blah')  # file name should be irrelevant
                img = fits.PrimaryHDU(np.arange(100))
                img.header['EXTNAME'] = 'SCI'
                hdul = fits.HDUList([img])
                hdul.writeto(img_file)

                tabimg_file = os.path.join(tmpdir, 'both.fits')
                table = fits.BinTableHDU.from_columns([
                    fits.Column(name='col1', format='I', array=np.array([1, 2, 3])),
                    fits.Column(name='col2', format='I', array=np.array([4, 5, 6]))
                ])
                hdul = fits.HDUList([img, table])
                hdul.writeto(tabimg_file)
                self.data_product.data = tabimg_file
                self.assertTrue(is_fits_image_file(self.data_product.data.file)) 

def make_simple_hdulist(self):
        """
        Make a`~astropy.io.fits.HDUList` object with just a simple
        PrimaryHDU
        """
        p = pyfits.PrimaryHDU()
        p.header['ID'] = (self.id, 'Object ID')
        p.header['RA'] = (self.ra, 'R.A.')
        p.header['DEC'] = (self.dec, 'Decl.')
        p.header['NINPUT'] = (len(self.beams), 'Number of drizzled beams')
        p.header['HASLINES'] = ('', 'Lines in this file')

        for i, beam in enumerate(self.beams):
            p.header['FILE{0:04d}'.format(i+1)] = (beam.grism.parent_file,
                                                 'Parent filename')
            p.header['GRIS{0:04d}'.format(i+1)] = (beam.grism.filter,
                                                 'Beam grism element')
            p.header['PA{0:04d}'.format(i+1)] = (beam.get_dispersion_PA(),
                                                 'PA of dispersion axis')

        return pyfits.HDUList(p) 

def update_wcs_fits_log(file, wcs_ref, xyscale=[0, 0, 0, 1], initialize=True, replace=('.fits', '.wcslog.fits'), wcsname='SHIFT'):
    """
    Make FITS log when updating WCS
    """
    new_hdu = wcs_ref.to_fits(relax=True)[0]
    new_hdu.header['XSHIFT'] = xyscale[0]
    new_hdu.header['YSHIFT'] = xyscale[1]
    new_hdu.header['ROT'] = xyscale[2], 'WCS fit rotation, degrees'
    new_hdu.header['SCALE'] = xyscale[3], 'WCS fit scale'
    new_hdu.header['WCSNAME'] = wcsname

    wcs_logfile = file.replace(replace[0], replace[1])

    if os.path.exists(wcs_logfile):
        if initialize:
            os.remove(wcs_logfile)
            hdu = pyfits.HDUList([pyfits.PrimaryHDU()])
        else:
            hdu = pyfits.open(wcs_logfile)
    else:
        hdu = pyfits.HDUList([pyfits.PrimaryHDU()])

    hdu.append(new_hdu)
    hdu.writeto(wcs_logfile, overwrite=True, output_verify='fix') 

def _make_primary_hdu(self, target_id):
        """Returns the primary extension of a Target Pixel File."""
        hdu = fits.PrimaryHDU()
        # Copy the default keywords from a template file from the MAST archive
        tmpl = self.get_header_template(0)
        for kw in tmpl:
            hdu.header[kw] = (tmpl[kw], tmpl.comments[kw])
        # Override the defaults where necessary
        hdu.header['ORIGIN'] = "Unofficial data product"
        hdu.header['DATE'] = datetime.datetime.now().strftime("%Y-%m-%d")
        hdu.header['CREATOR'] = "kadenza"
        hdu.header['PROCVER'] = "{}".format(__version__)
        hdu.header['FILEVER'] = "0.0"
        hdu.header['OBJECT'] = target_name(target_id)
        hdu.header['KEPLERID'] = target_id
        hdu.header['CHANNEL'] = self.pixel_mapping.targets[target_id]['channel']
        hdu.header['MODULE'] = self.pixel_mapping.targets[target_id]['module']
        hdu.header['OUTPUT'] = self.pixel_mapping.targets[target_id]['output']
        # Empty a bunch of keywords rather than having incorrect info
        for kw in ["TIMVERSN", "CAMPAIGN", "DATA_REL", "TTABLEID",
                   "RA_OBJ", "DEC_OBJ"]:
            hdu.header[kw] = ""
        return hdu 

def write_res(filename, datas, extnames, header='', hdrref=None, clobber=False):
   if not header and hdrref: header = pyfits.getheader(hdrref)
   hdu = pyfits.PrimaryHDU(header=header)
   warnings.resetwarnings() # supress nasty overwrite warning http://pythonhosted.org/pyfits/users_guide/users_misc.html
   warnings.filterwarnings('ignore', category=UserWarning, append=True)
   hdu.writeto(filename, clobber=clobber, output_verify='fix')
   warnings.resetwarnings()
   warnings.filterwarnings('always', category=UserWarning, append=True)

   for i,extname in enumerate(extnames):
     data = datas[extname]
     if isinstance(data, np.ndarray):
        pyfits.append(filename, data)
     else:
        1/0

     pyfits.setval(filename, 'EXTNAME', value=extname, ext=i+1)
   #fitsio.write(filename, flux) 

def write_fits(data, header, file_name):
    """
    Combine data and a fits header to write a fits file.

    Parameters
    ----------
    data : numpy.ndarray
        The data to be written.

    header : astropy.io.fits.hduheader
        The header for the fits file.

    file_name : string
        The file to write

    Returns
    -------
    None
    """
    hdu = fits.PrimaryHDU(data)
    hdu.header = header
    hdulist = fits.HDUList([hdu])
    hdulist.writeto(file_name, overwrite=True)
    logging.info("Wrote {0}".format(file_name))
    return 

def saveScrns(self, DIR):
        """
        Saves the currently loaded phase screens to file,
        saving the r0 value in the fits header (in units of pixels). 
        Saved phase data is in radians @500nm

        Args:
            DIR (string): The directory to save the screens
        """

        for scrn in range(self.scrnNo):
            logger.info("Write Sreen {}....".format(scrn))
            hdu = fits.PrimaryHDU(self.wholeScrns[scrn])
            hdu.header["R0"] = self.wholeScrnR0 / self.pixel_scale
            hdulist = fits.HDUList([hdu])
            hdulist.writeto(DIR+"/scrn{}.fits".format(scrn))
            hdulist.close()

            logger.info("Done!") 

def __init__(self, primary_hdu=None, fits_extensions=None):

        hdu_list = []

        if primary_hdu is None:

            primary_hdu = fits.PrimaryHDU()

        else:

            assert isinstance(primary_hdu, fits.PrimaryHDU)

        hdu_list.append(primary_hdu)

        if fits_extensions is not None:

            fits_extensions = list(fits_extensions)

            hdu_list.extend([x.hdu for x in fits_extensions])

        # We embed instead of subclassing because the HDUList class has some weird interaction with the
        # __init__ and __new__ methods which makes difficult to do so (we couldn't figure it out)

        self._hdu_list = fits.HDUList(hdus=hdu_list) 

def pack(uncompressed_hdulist: fits.HDUList) -> fits.HDUList:
    if uncompressed_hdulist[0].data is None:
        primary_hdu = fits.PrimaryHDU(header=uncompressed_hdulist[0].header)
        hdulist = [primary_hdu]
    else:
        primary_hdu = fits.PrimaryHDU()
        compressed_hdu = fits.CompImageHDU(data=np.ascontiguousarray(uncompressed_hdulist[0].data),
                                           header=uncompressed_hdulist[0].header, quantize_level=64,
                                           dither_seed=2048, quantize_method=1)
        hdulist = [primary_hdu, compressed_hdu]

    for hdu in uncompressed_hdulist[1:]:
        if isinstance(hdu, fits.ImageHDU):
            compressed_hdu = fits.CompImageHDU(data=np.ascontiguousarray(hdu.data), header=hdu.header,
                                               quantize_level=64, quantize_method=1)
            hdulist.append(compressed_hdu)
        else:
            hdulist.append(hdu)
    return fits.HDUList(hdulist) 

def save_single_fits(image, name, key_dict=None, image2=None, image2type=None):
        # Save an array into the first extension of a fits file
        h0 = fits.PrimaryHDU()
        h1 = fits.ImageHDU(image, name='DATA')
        if image2 is not None:
            h2 = fits.ImageHDU(image2)
            if image2type is not None:
                h2.header['EXTNAME'] = image2type

        # if a keyword dictionary is provided, put the
        # keywords into the 0th and 1st extension headers
        if key_dict is not None:
            for key in key_dict:
                h0.header[key] = key_dict[key]
                h1.header[key] = key_dict[key]

        if image2 is None:
            hdulist = fits.HDUList([h0, h1])
        else:
            hdulist = fits.HDUList([h0, h1, h2])
        hdulist.writeto(name, overwrite=True) 

def test_read_wcs_not_creatable(tmpdir):
    # The following Header can't be converted to a WCS object. See also #6499.
    hdr_txt_example_WCS = textwrap.dedent('''
    SIMPLE  =                    T / Fits standard
    BITPIX  =                   16 / Bits per pixel
    NAXIS   =                    2 / Number of axes
    NAXIS1  =                 1104 / Axis length
    NAXIS2  =                 4241 / Axis length
    CRVAL1  =         164.98110962 / Physical value of the reference pixel X
    CRVAL2  =          44.34089279 / Physical value of the reference pixel Y
    CRPIX1  =                -34.0 / Reference pixel in X (pixel)
    CRPIX2  =               2041.0 / Reference pixel in Y (pixel)
    CDELT1  =           0.10380000 / X Scale projected on detector (#/pix)
    CDELT2  =           0.10380000 / Y Scale projected on detector (#/pix)
    CTYPE1  = 'RA---TAN'           / Pixel coordinate system
    CTYPE2  = 'WAVELENGTH'         / Pixel coordinate system
    CUNIT1  = 'degree  '           / Units used in both CRVAL1 and CDELT1
    CUNIT2  = 'nm      '           / Units used in both CRVAL2 and CDELT2
    CD1_1   =           0.20760000 / Pixel Coordinate translation matrix
    CD1_2   =           0.00000000 / Pixel Coordinate translation matrix
    CD2_1   =           0.00000000 / Pixel Coordinate translation matrix
    CD2_2   =           0.10380000 / Pixel Coordinate translation matrix
    C2YPE1  = 'RA---TAN'           / Pixel coordinate system
    C2YPE2  = 'DEC--TAN'           / Pixel coordinate system
    C2NIT1  = 'degree  '           / Units used in both C2VAL1 and C2ELT1
    C2NIT2  = 'degree  '           / Units used in both C2VAL2 and C2ELT2
    RADECSYS= 'FK5     '           / The equatorial coordinate system
    ''')
    with catch_warnings(FITSFixedWarning):
        hdr = fits.Header.fromstring(hdr_txt_example_WCS, sep='\n')
    hdul = fits.HDUList([fits.PrimaryHDU(np.ones((4241, 1104)), header=hdr)])
    filename = tmpdir.join('afile.fits').strpath
    hdul.writeto(filename)
    # The hdr cannot be converted to a WCS object because of an
    # InconsistentAxisTypesError but it should still open the file
    ccd = CCDData.read(filename, unit='adu')
    assert ccd.wcs is None 

def test_history_preserved_if_metadata_is_fits_header(tmpdir):
    fake_img = np.zeros([2, 2])
    hdu = fits.PrimaryHDU(fake_img)
    hdu.header['history'] = 'one'
    hdu.header['history'] = 'two'
    hdu.header['history'] = 'three'
    assert len(hdu.header['history']) == 3
    tmp_file = tmpdir.join('temp.fits').strpath
    hdu.writeto(tmp_file)

    ccd_read = CCDData.read(tmp_file, unit="adu")
    assert ccd_read.header['history'] == hdu.header['history'] 

def test_read_returns_image(tmpdir):
    # Test if CCData.read returns a image when reading a fits file containing
    # a table and image, in that order.
    tbl = Table(np.ones(10).reshape(5, 2))
    img = np.ones((5, 5))
    hdul = fits.HDUList(hdus=[fits.PrimaryHDU(), fits.TableHDU(tbl.as_array()),
                              fits.ImageHDU(img)])
    filename = tmpdir.join('table_image.fits').strpath
    hdul.writeto(filename)
    ccd = CCDData.read(filename, unit='adu')
    # Expecting to get (5, 5), the size of the image
    assert ccd.data.shape == (5, 5)


# https://github.com/astropy/astropy/issues/9664 

def test_passing_ImageHDU():
    """
    Passing ImageHDU or PrimaryHDU and comparing it with
    wcs initialized from header. For #4493.
    """
    path = get_pkg_data_filename('data/validate.fits')
    with fits.open(path) as hdulist:
        with pytest.warns(FITSFixedWarning):
            wcs_hdu = wcs.WCS(hdulist[0])
            wcs_header = wcs.WCS(hdulist[0].header)
            assert wcs_hdu.wcs.compare(wcs_header.wcs)
            wcs_hdu = wcs.WCS(hdulist[1])
            wcs_header = wcs.WCS(hdulist[1].header)
            assert wcs_hdu.wcs.compare(wcs_header.wcs) 

def test_append_primary_to_non_empty_list(self):
        """Tests appending a Simple PrimaryHDU to a non-empty HDUList."""

        with fits.open(self.data('arange.fits')) as hdul:
            hdu = fits.PrimaryHDU(np.arange(100, dtype=np.int32))
            hdul.append(hdu)

            info = [(0, 'PRIMARY', 1, 'PrimaryHDU', 7, (11, 10, 7), 'int32', ''),
                    (1, '', 1, 'ImageHDU', 6, (100,), 'int32', '')]

            assert hdul.info(output=False) == info

            hdul.writeto(self.temp('test-append.fits'))

        assert fits.info(self.temp('test-append.fits'), output=False) == info 

def test_validate_faulty_wcs():
    """
    From github issue #2053
    """
    h = fits.Header()
    # Illegal WCS:
    h['RADESYSA'] = 'ICRS'
    h['PV2_1'] = 1.0
    hdu = fits.PrimaryHDU([[0]], header=h)
    hdulist = fits.HDUList([hdu])
    # Check that this doesn't raise a NameError exception
    wcs.validate(hdulist) 

def test_append_extension_to_non_empty_list(self):
        """Tests appending a Simple ExtensionHDU to a non-empty HDUList."""

        with fits.open(self.data('tb.fits')) as hdul:
            hdul.append(hdul[1])

            info = [(0, 'PRIMARY', 1, 'PrimaryHDU', 11, (), '', ''),
                    (1, '', 1, 'BinTableHDU', 24, '2R x 4C', '[1J, 3A, 1E, 1L]', ''),
                    (2, '', 1, 'BinTableHDU', 24, '2R x 4C', '[1J, 3A, 1E, 1L]', '')]

            assert hdul.info(output=False) == info

            hdul.writeto(self.temp('test-append.fits'))

        assert fits.info(self.temp('test-append.fits'), output=False) == info 

def test_to_fits():
    w = wcs.WCS()
    header_string = w.to_header()
    wfits = w.to_fits()
    assert isinstance(wfits, fits.HDUList)
    assert isinstance(wfits[0], fits.PrimaryHDU)
    assert header_string == wfits[0].header[-8:] 

def test_append_groupshdu_to_non_empty_list(self):
        """Tests appending a Simple GroupsHDU to an empty HDUList."""

        hdul = fits.HDUList()
        hdu = fits.PrimaryHDU(np.arange(100, dtype=np.int32))
        hdul.append(hdu)
        hdu = fits.GroupsHDU()
        hdul.append(hdu) 

def testMaskLoad():

    sim = soapy.Sim(os.path.join(CONFIG_PATH, "sh_8x8.yaml"))
    sim.config.sim.simName = None
    sim.config.sim.logfile = None

    sim.aoinit()

    mask = numpy.ones((sim.config.sim.pupilSize, sim.config.sim.pupilSize))

    # save mask
    if os.path.isfile('testmask.fits'):
        os.remove('testmask.fits')
    
    hdu = fits.PrimaryHDU(mask)
    hdulist = fits.HDUList([hdu])
    hdulist.writeto('testmask.fits')
    hdulist.close()
    
    try:
        # attempt to load it
        sim.config.tel.mask = 'testmask.fits'
        sim.aoinit()

        # check its good
        p = sim.config.sim.simPad
        pad_mask = numpy.pad(mask, mode="constant", pad_width=((p,p),(p,p)))
        assert numpy.array_equal(sim.mask, pad_mask)
    except:
        raise
    finally:
        os.remove('testmask.fits') 

def test_insert_primary_to_empty_list(self):
        """Tests inserting a Simple PrimaryHDU to an empty HDUList."""
        hdul = fits.HDUList()
        hdu = fits.PrimaryHDU(np.arange(100, dtype=np.int32))
        hdul.insert(0, hdu)

        info = [(0, 'PRIMARY', 1, 'PrimaryHDU', 5, (100,), 'int32', '')]

        assert hdul.info(output=False) == info

        hdul.writeto(self.temp('test-insert.fits'))

        assert fits.info(self.temp('test-insert.fits'), output=False) == info 

def test_insert_extension_to_empty_list(self):
        """Tests inserting a Simple ImageHDU to an empty HDUList."""

        hdul = fits.HDUList()
        hdu = fits.ImageHDU(np.arange(100, dtype=np.int32))
        hdul.insert(0, hdu)

        info = [(0, 'PRIMARY', 1, 'PrimaryHDU', 4, (100,), 'int32', '')]

        assert hdul.info(output=False) == info

        hdul.writeto(self.temp('test-insert.fits'))

        assert fits.info(self.temp('test-insert.fits'), output=False) == info 

def test_insert_primary_to_non_empty_list(self):
        """Tests inserting a Simple PrimaryHDU to a non-empty HDUList."""

        with fits.open(self.data('arange.fits')) as hdul:
            hdu = fits.PrimaryHDU(np.arange(100, dtype=np.int32))
            hdul.insert(1, hdu)

            info = [(0, 'PRIMARY', 1, 'PrimaryHDU', 7, (11, 10, 7), 'int32', ''),
                    (1, '', 1, 'ImageHDU', 6, (100,), 'int32', '')]

            assert hdul.info(output=False) == info

            hdul.writeto(self.temp('test-insert.fits'))

        assert fits.info(self.temp('test-insert.fits'), output=False) == info 

def _load_itm_library(library_file):
    """Load ITM FITS file

    Parameters
    ----------
    library_path : str
        Path pointing to the location of the PSF library

    Returns
    -------
    library : photutils.griddedPSFModel
        Object containing PSF library
    """
    data = fits.getdata(library_file)
    hdr = fits.getheader(library_file)
    if data.shape == (2048, 2048):
        # Normalize the data
        data /= np.sum(data)

        # Add PSF location and oversampling keywords
        hdr['DET_YX0'] = ('(1023, 1023)', "The #0 PSF's (y,x) detector pixel position")
        hdr['OVERSAMP'] = (1, 'Oversampling factor for FFTs in computation')

        # Convert to HDUList and create library
        phdu = fits.PrimaryHDU(data, hdr)
        hdulist = fits.HDUList(phdu)
        library = to_griddedpsfmodel(hdulist)

        return library
    else:
        raise ValueError('Expecting ITM data of size (2048, 2048), not {}'.format(data.shape)) 

def test_insert_extension_to_non_empty_list(self):
        """Tests inserting a Simple ExtensionHDU to a non-empty HDUList."""

        with fits.open(self.data('tb.fits')) as hdul:
            hdul.insert(1, hdul[1])

            info = [(0, 'PRIMARY', 1, 'PrimaryHDU', 11, (), '', ''),
                    (1, '', 1, 'BinTableHDU', 24, '2R x 4C', '[1J, 3A, 1E, 1L]', ''),
                    (2, '', 1, 'BinTableHDU', 24, '2R x 4C', '[1J, 3A, 1E, 1L]', '')]

            assert hdul.info(output=False) == info

            hdul.writeto(self.temp('test-insert.fits'))

        assert fits.info(self.temp('test-insert.fits'), output=False) == info 

def savefits(self, array, ofile):
        """Save the given array in a fits file

        Parameters
        ----------
        array : numpy.ndarray
            2D array containing the data to be saved

        ofile : str
            Output filename to use
        """
        h0 = fits.PrimaryHDU()
        h1 = fits.ImageHDU()
        h1.data = array

        h1.header['INSTRUME'] = self.instrument
        h1.header['CDELT1'] = self.mosaic_scale_x / 3600.
        h1.header['CDELT2'] = self.mosaic_scale_y / 3600.
        h1.header['CRPIX1'] = self.crpix1
        h1.header['CRPIX2'] = self.crpix2
        h1.header['CRVAL1'] = self.center_ra
        h1.header['CRVAL2'] = self.center_dec
        h1.header['CTYPE1'] = self.mosaic_x_type
        h1.header['CTYPE2'] = self.mosaic_y_type
        h1.header['CUNIT1'] = 'deg'
        h1.header['CUNIT2'] = 'deg'
        h1.header['RA_REF'] = self.center_ra
        h1.header['DEC_REF'] = self.center_dec
        h1.header['ROLL_REF'] = self.mosaic_roll * 180./np.pi
        h1.header['WCSAXES'] = 2
        h1.header['EXTNAME'] = 'SCI    '
        h1.header['EXTVER'] = 1
        h1.header['PC1_1'] = self.cd11 / (self.mosaic_scale_x / 3600.)
        h1.header['PC1_2'] = self.cd12 / (self.mosaic_scale_x / 3600.)
        h1.header['PC2_1'] = self.cd21 / (self.mosaic_scale_y / 3600.)
        h1.header['PC2_2'] = self.cd22 / (self.mosaic_scale_y / 3600.)

        hlist = fits.HDUList([h0, h1])
        hlist.writeto(ofile, overwrite=True) 

def test_to_fits_1():
    """
    Test to_fits() with LookupTable distortion.
    """
    fits_name = get_pkg_data_filename('data/dist.fits')
    with pytest.warns(AstropyDeprecationWarning):
        w = wcs.WCS(fits_name)
    wfits = w.to_fits()
    assert isinstance(wfits, fits.HDUList)
    assert isinstance(wfits[0], fits.PrimaryHDU)
    assert isinstance(wfits[1], fits.ImageHDU)