Python astropy.io.fits.open() Examples

def load_allstar_dr5():
    """
    Open LAMOST DR5 allstar

    :return: fits file opened by astropy
    :rtype: astropy.io.fits.hdu.hdulist.HDUList
    :History: 2018-Jun-17 - Written - Henry Leung (University of Toronto)
    """
    import os
    from astropy.io import fits

    file_name = "LAMO5_2MS_AP9_SD14_UC4_PS1_AW_Carlin_M.fits"
    _lamost_dr5_allsta_path = os.path.join(lamost_env(), "DR5", file_name)
    if not os.path.isfile(_lamost_dr5_allsta_path):
        raise FileNotFoundError(f'{file_name} file not found')
    return fits.open(_lamost_dr5_allsta_path) 

def test_fits_read_textfile_exists(self) -> None:

        with open(self.test_dir+'filenames.dat', 'w') as file_obj:
            file_obj.write(self.test_dir+'fits/images_0.fits\n')
            file_obj.write(self.test_dir+'fits/images_2.fits\n')

        module = FitsReadingModule(name_in='read10',
                                   input_dir=None,
                                   image_tag='files',
                                   overwrite=True,
                                   check=True,
                                   filenames=self.test_dir+'filenames.dat')

        self.pipeline.add_module(module)

        with pytest.raises(ValueError) as error:
            self.pipeline.run_module('read10')

        assert str(error.value) == f'The file {self.test_dir}fits/images_2.fits does not exist. ' \
                                   f'Please check that the path is correct.' 

def get_PAM_value(self, verbose=False):
        """
        Apply Pixel Area Map correction to WFC3 effective PSF model

        http://www.stsci.edu/hst/wfc3/pam/pixel_area_maps
        """
        confp = self.conf.conf
        if ('INSTRUMENT' in confp) & ('CAMERA' in confp):
            if '{0}-{1}'.format(confp['INSTRUMENT'], confp['CAMERA']) != 'WFC3-IR':
                return 1
        else:
            return 1

        try:
            pam = pyfits.open(os.getenv('iref')+'ir_wfc3_map.fits')
            pam_data = pam[1].data
            pam_value = pam_data[int(self.yc-self.pad), int(self.xc-self.pad)]
            pam.close()
        except:
            pam_value = 1

        if verbose:
            print('PAM correction at x={0}, y={1}: {2:.3f}'.format(self.xc-self.pad, self.yc-self.pad, pam_value))

        return pam_value 

def test_header_attribute(self) -> None:

        with fits.open(self.test_dir+'fits/images_0.fits') as hdu:
            header = hdu[0].header
            header['PARANG'] = 1.0
            hdu.writeto(self.test_dir+'fits/images_0.fits', overwrite=True)

        with fits.open(self.test_dir+'fits/images_1.fits') as hdu:
            header = hdu[0].header
            header['PARANG'] = 2.0
            header['HIERARCH ESO DET DIT'] = 0.1
            hdu.writeto(self.test_dir+'fits/images_1.fits', overwrite=True)

        module = FitsReadingModule(name_in='read5',
                                   input_dir=self.test_dir+'fits',
                                   image_tag='input',
                                   overwrite=True,
                                   check=True)

        self.pipeline.add_module(module)
        self.pipeline.run_module('read5') 

def test_odd_number_images(self) -> None:

        with fits.open(self.fitsfile) as hdulist:
            del hdulist[11]
            hdulist.writeto(self.fitsfile, overwrite=True)

        module = NearReadingModule(name_in='read8',
                                   input_dir=self.test_dir+'near',
                                   chopa_out_tag=self.positions[0],
                                   chopb_out_tag=self.positions[1])

        self.pipeline.add_module(module)

        with pytest.warns(UserWarning) as warning:
            self.pipeline.run_module('read8')

        assert len(warning) == 2

        assert warning[0].message.args[0] == f'FITS file contains odd number of images: ' \
                                             f'{self.fitsfile}'

        assert warning[1].message.args[0] == 'The number of chop cycles (5) is not equal to ' \
                                             'half the number of available HDU images (4).' 

def load_spectrum(filename, grid):
    """
    Load a single spectrum
    """
    file_in = pyfits.open(filename)
    wl = np.array(file_in[0].data[2])
    flux = np.array(file_in[0].data[0])
    ivar = np.array((file_in[0].data[1]))
    # correct for radial velocity of star
    redshift = file_in[0].header['Z']
    wl_shifted = wl - redshift * wl
    # resample
    flux_rs = (interpolate.interp1d(wl_shifted, flux))(grid)
    ivar_rs = (interpolate.interp1d(wl_shifted, ivar))(grid)
    ivar_rs[ivar_rs < 0] = 0. # in interpolating you can end up with neg
    return flux_rs, ivar_rs 

def test_frame_type_missing(self) -> None:

        with fits.open(self.fitsfile) as hdulist:
            hdulist[10].header.remove('ESO DET FRAM TYPE')
            hdulist.writeto(self.fitsfile, overwrite=True)

        module = NearReadingModule(name_in='read6',
                                   input_dir=self.test_dir+'near',
                                   chopa_out_tag=self.positions[0],
                                   chopb_out_tag=self.positions[1])

        self.pipeline.add_module(module)

        with pytest.raises(ValueError) as error:
            self.pipeline.run_module('read6')

        assert str(error.value) == 'Frame type not found in the FITS header. Image number: 9.' 

def test_frame_type_invalid(self) -> None:

        with fits.open(self.fitsfile) as hdulist:
            hdulist[10].header['ESO DET FRAM TYPE'] = 'Test'
            hdulist.writeto(self.fitsfile, overwrite=True)

        module = NearReadingModule(name_in='read5',
                                   input_dir=self.test_dir+'near',
                                   chopa_out_tag=self.positions[0],
                                   chopb_out_tag=self.positions[1])

        self.pipeline.add_module(module)

        with pytest.raises(ValueError) as error:
            self.pipeline.run_module('read5')

        assert str(error.value) == 'Frame type (Test) not a valid value. Expecting HCYCLE1 or ' \
                                   'HCYCLE2 as value for ESO DET FRAM TYPE.' 

def __init__(self, base_fname):
        """
        Args:
            base_fname (str): The map should be stored in two FITS files, named
                ``base_fname + '_' + X + '.fits'``, where ``X`` is ``'ngp'`` and
                ``'sgp'``.
        """
        self._data = {}

        for pole in self.poles:
            fname = '{}_{}.fits'.format(base_fname, pole)
            try:
                with fits.open(fname) as hdulist:
                    self._data[pole] = [hdulist[0].data, wcs.WCS(hdulist[0].header)]
            except IOError as error:
                print(dustexceptions.data_missing_message(self.map_name,
                                                          self.map_name_long))
                raise error 

def test_fine_channels(tmp_path):
    """
    Break read_mwa_corr_fits by submitting files with different fine channels.

    Test that error is raised if files with different numbers of fine channels
    are submitted.
    """
    mwa_uv = UVData()
    bad_fine = str(tmp_path / "bad_gpubox06_01.fits")
    with fits.open(filelist[2]) as mini6:
        mini6[1].data = np.concatenate((mini6[1].data, mini6[1].data))
        mini6.writeto(bad_fine)
    with pytest.raises(ValueError) as cm:
        mwa_uv.read([bad_fine, filelist[1]])
    assert str(cm.value).startswith("files submitted have different fine")
    del mwa_uv 

def test_diff_obs(tmp_path):
    """
    Break read_mwa_corr_fits by submitting files from different observations.

    Test that error is raised if files from different observations are
    submitted in the same file list.
    """
    mwa_uv = UVData()
    bad_obs = str(tmp_path / "bad2_gpubox06_01.fits")
    with fits.open(filelist[2]) as mini6:
        mini6[0].header["OBSID"] = "1131733555"
        mini6.writeto(bad_obs)
    with pytest.raises(ValueError) as cm:
        mwa_uv.read([bad_obs, filelist[0], filelist[1]])
    assert str(cm.value).startswith("files from different observations")
    del mwa_uv 

def main(infile, xypadsize):
    """Pad a fits image with zeros to padsize"""
    pad_xsize, pad_ysize = [int(s) for s in xypadsize.split(' ')]

    hdu = pyfits.open(infile)
    imdata = hdu[0].data[0, 0]
    (ysize, xsize) = imdata.shape
    if xsize > pad_xsize or ysize > pad_ysize:
        raise ValueError('pad_image: padded size is smaller than current size!')
    if xsize == pad_xsize and ysize == pad_ysize:
        return()

    xoffset = (pad_xsize - xsize) / 2
    yoffset = (pad_ysize - ysize) / 2
    newdata=np.zeros((1, 1, pad_ysize, pad_xsize))

    newdata[0, 0, yoffset:yoffset+ysize, xoffset:xoffset+xsize] = imdata
    hdu[0].data = newdata
    hdu[0].header['CRPIX1'] += xoffset
    hdu[0].header['CRPIX2'] += yoffset
    hdu.writeto(infile, clobber=True) 

def parse_input_healpix_data(input_data, field=0, hdu_in=None, nested=None):
    """
    Parse input HEALPIX data to return a Numpy array and coordinate frame object.
    """

    if isinstance(input_data, (TableHDU, BinTableHDU)):
        data = input_data.data
        header = input_data.header
        coordinate_system_in = parse_coord_system(header['COORDSYS'])
        array_in = data[data.columns[field].name].ravel()
        if 'ORDERING' in header:
            nested = header['ORDERING'].lower() == 'nested'
    elif isinstance(input_data, str):
        hdu = fits.open(input_data)[hdu_in or 1]
        return parse_input_healpix_data(hdu, field=field)
    elif isinstance(input_data, tuple) and isinstance(input_data[0], np.ndarray):
        array_in = input_data[0]
        coordinate_system_in = parse_coord_system(input_data[1])
    else:
        raise TypeError("input_data should either be an HDU object or a tuple of (array, frame)")

    return array_in, coordinate_system_in, nested 

def find_fits_img_size(filename):
    """
    Returns the size of a FITS image, given a valid FITS image file

    :param filename: The fully-qualified path of the FITS image file
    :type filename: str

    :returns: Tuple of horizontal/vertical dimensions
    :rtype: tuple
    """

    try:
        return settings.THUMBNAIL_MAX_SIZE
    except AttributeError:
        hdul = fits.open(filename)
        xsize = 0
        ysize = 0
        for hdu in hdul:
            try:
                xsize = max(xsize, hdu.header['NAXIS1'])
                ysize = max(ysize, hdu.header['NAXIS2'])
            except KeyError:
                pass
        return (xsize, ysize) 

def is_fits_image_file(file):
    """
    Checks if a file is a valid FITS image by checking if any header contains 'SCI' in the 'EXTNAME'.

    :param file: The file to be checked.
    :type file:

    :returns: True if the file is a FITS image, False otherwise
    :rtype: boolean
    """

    with file.open() as f:
        try:
            hdul = fits.open(f)
        except OSError:  # OSError is raised if file is not FITS format
            return False
        for hdu in hdul:
            if hdu.header.get('EXTNAME') == 'SCI':
                return True
    return False 

def add_all_spectra():

    from grizli.aws import db as grizli_db

    roots = grizli_db.from_sql("select root,count(root) as n from redshift_fit group BY root order by n DESC", engine)

    o = 1

    for root in roots['root'][::o]:
        existing = open('log').readlines()
        if root+'\n' in existing:
            print('Skip', root)
            continue

        fp = open('log', 'a')
        fp.write(root+'\n')
        fp.close()
        try:
            grizli_db.add_oned_spectra(root=root, engine=engine)
        except:
            pass 

def test_readwriteread_missing_info(tmp_path, casa_uvfits):
    uv_in = casa_uvfits
    uv_out = UVData()
    write_file = str(tmp_path / "outtest_casa.uvfits")
    write_file2 = str(tmp_path / "outtest_casa2.uvfits")

    # check missing telescope_name, timesys vs timsys spelling, xyz_telescope_frame=????
    uv_in.write_uvfits(write_file)
    with fits.open(write_file, memmap=True) as hdu_list:
        hdunames = uvutils._fits_indexhdus(hdu_list)
        vis_hdu = hdu_list[0]
        vis_hdr = vis_hdu.header.copy()

        vis_hdr.pop("TELESCOP")

        vis_hdu.header = vis_hdr

        ant_hdu = hdu_list[hdunames["AIPS AN"]]
        ant_hdr = ant_hdu.header.copy()

        time_sys = ant_hdr.pop("TIMSYS")
        ant_hdr["TIMESYS"] = time_sys
        ant_hdr["FRAME"] = "????"

        ant_hdu.header = ant_hdr

        hdulist = fits.HDUList(hdus=[vis_hdu, ant_hdu])
        hdulist.writeto(write_file2, overwrite=True)

    uv_out.read(write_file2)
    assert uv_out.telescope_name == "EVLA"
    assert uv_out.timesys == time_sys

    return 

def check_object_in_footprint(id, wcs_fits, cat, rd=None):
    """
    Check that a given object is within a WCS footprint

    Parameters
    ----------
    id : int
        Integer id of the catalog objects

    wcs_fits : str
        WCS filename, like 'ibfuw1psq_flt.01.wcs.fits'.

    cat : `~astropy.table.Table`
        Table object of the "ir_cat.fits" source detection table, with
        SExtractor columns NUMBER, X_WORLD, Y_WORLD.

    """
    import matplotlib.path
    import astropy.wcs as pywcs
    import astropy.io.fits as pyfits

    if rd is None:
        ix = cat['NUMBER'] == id
        ra, dec = cat['X_WORLD'][ix][0], cat['Y_WORLD'][ix][0]
    else:
        ra, dec = rd

    im = pyfits.open(wcs_fits)
    im[0].header['NAXIS'] = 2
    im[0].header['NAXIS1'] = im[0].header['CRPIX1']*2
    im[0].header['NAXIS2'] = im[0].header['CRPIX2']*2

    wcs = pywcs.WCS(im[0].header, fobj=im, relax=True)
    fp = matplotlib.path.Path(wcs.calc_footprint())
    has_point = fp.contains_point([ra, dec])
    im.close()

    return has_point 

def get_preview(self, size=THUMBNAIL_DEFAULT_SIZE, redraw=False):
        """
        Returns path to the thumbnail of this data product, and creates a thumbnail if none exists

       :Keyword Arguments:
            * size (`tuple`): Desired size of the thumbnail, as a 2-tuple of ints for width/height
            * redraw (`boolean`): True if the thumbnail will be recreated despite existing, False otherwise

        :returns: Path to the thumbnail image
        :rtype: str
        """
        if self.thumbnail:
            im = Image.open(self.thumbnail)
            if im.size != THUMBNAIL_DEFAULT_SIZE:
                redraw = True

        if not self.thumbnail or redraw:
            width, height = THUMBNAIL_DEFAULT_SIZE
            tmpfile = self.create_thumbnail(width=width, height=height)
            if tmpfile:
                outfile_name = os.path.basename(self.data.file.name)
                filename = outfile_name.split(".")[0] + "_tb.jpg"
                with open(tmpfile.name, 'rb') as f:
                    self.thumbnail.save(filename, File(f), save=True)
                    self.save()
        if not self.thumbnail:
            return ''
        return self.thumbnail.url 

def test_read_noversion_history(tmp_path):
    """
    Test that version info gets added to the history if it's missing
    """
    cal_in = UVCal()
    cal_out = UVCal()
    testfile = os.path.join(DATA_PATH, "zen.2457698.40355.xx.gain.calfits")
    write_file = str(tmp_path / "outtest_omnical.fits")
    write_file2 = str(tmp_path / "outtest_omnical2.fits")
    cal_in.read_calfits(testfile)

    cal_in.write_calfits(write_file, clobber=True)

    fname = fits.open(write_file)
    data = fname[0].data
    primary_hdr = fname[0].header
    hdunames = uvutils._fits_indexhdus(fname)
    ant_hdu = fname[hdunames["ANTENNAS"]]

    primary_hdr["HISTORY"] = ""

    prihdu = fits.PrimaryHDU(data=data, header=primary_hdr)
    hdulist = fits.HDUList([prihdu, ant_hdu])

    hdulist.writeto(write_file2, overwrite=True)
    hdulist.close()

    cal_out.read_calfits(write_file2)
    assert cal_in == cal_out 

def flag_file_init(tmp_path_factory):
    tmp_path = tmp_path_factory.mktemp("pyuvdata_corr_fits", numbered=True)
    spoof_file1 = str(tmp_path / "spoof_01_00.fits")
    spoof_file6 = str(tmp_path / "spoof_06_00.fits")
    # spoof box files of the appropriate size
    with fits.open(filelist[1]) as mini1:
        mini1[1].data = np.repeat(mini1[1].data, 8, axis=0)
        extra_dat = np.copy(mini1[1].data)
        for app_ind in range(2):
            mini1.append(fits.ImageHDU(extra_dat))
        mini1[2].header["MILLITIM"] = 500
        mini1[2].header["TIME"] = mini1[1].header["TIME"]
        mini1[3].header["MILLITIM"] = 0
        mini1[3].header["TIME"] = mini1[1].header["TIME"] + 1
        mini1.writeto(spoof_file1)

    with fits.open(filelist[2]) as mini6:
        mini6[1].data = np.repeat(mini6[1].data, 8, axis=0)
        extra_dat = np.copy(mini6[1].data)
        for app_ind in range(2):
            mini6.append(fits.ImageHDU(extra_dat))
        mini6[2].header["MILLITIM"] = 500
        mini6[2].header["TIME"] = mini6[1].header["TIME"]
        mini6[3].header["MILLITIM"] = 0
        mini6[3].header["TIME"] = mini6[1].header["TIME"] + 1
        mini6.writeto(spoof_file6)

    flag_testfiles = [spoof_file1, spoof_file6, filelist[0]]

    yield flag_testfiles 

def get_connection_info(config_file=None):
    """
    Read the database connection info
    """
    import yaml

    if config_file is None:
        config_file = os.path.join(os.path.dirname(__file__),
                                   '../data/db.yml')

        try:
            local_file = os.path.join(os.getenv('HOME'), 'db.local.yml')

            if os.path.exists(local_file):
                print('Use ~/db.local.yml')
                config_file = local_file
        except:
            pass

    fp = open(config_file)
    try:
        db_info = yaml.load(fp, Loader=yaml.FullLoader)
    except:
        db_info = yaml.load(fp)

    fp.close()

    return db_info 

def load_photutils_detection(self, seg_file=None, seg_cat=None,
                                 catalog_format='ascii.commented_header'):
        """
        Load segmentation image and catalog, either from photutils
        or SExtractor.

        If SExtractor, use `catalog_format='ascii.sextractor'`.

        """
        root = self.direct_file.split('.fits')[0]

        if seg_file is None:
            seg_file = root + '.detect_seg.fits'

        if not os.path.exists(seg_file):
            print('Segmentation image {0} not found'.format(segfile))
            return False

        self.seg = np.cast[np.float32](pyfits.open(seg_file)[0].data)

        if seg_cat is None:
            seg_cat = root + '.detect.cat'

        if not os.path.exists(seg_cat):
            print('Segmentation catalog {0} not found'.format(seg_cat))
            return False

        self.catalog = Table.read(seg_cat, format=catalog_format)
        self.catalog_file = seg_cat 

def test_file_extension(tmp_path):
    """
    Break read_mwa_corr_fits by submitting file with the wrong extension.

    Test that error is raised if a file with an extension that is not fits,
    metafits, or mwaf is submitted.
    """
    mwa_uv = UVData()
    bad_ext = str(tmp_path / "1131733552.meta")
    with fits.open(filelist[0]) as meta:
        meta.writeto(bad_ext)
    with pytest.raises(ValueError) as cm:
        mwa_uv.read(bad_ext, file_type="mwa_corr_fits")
    assert str(cm.value).startswith("only fits, metafits, and mwaf files supported")
    del mwa_uv 

def make_mosaic_footprints(field_root):
    """
    Make region files where wht images nonzero
    """
    files = glob.glob('{0}-f*dr?_wht.fits'.format(field_root))
    files.sort()

    fp = open('{0}_mosaic.reg'.format(field_root), 'w')
    fp.write('fk5\n')
    fp.close()

    for weight_image in files:
        filt = weight_image.split('_dr')[0].split('-')[-1]

        wave = filt[1:4]
        if wave[0] in '01':
            w = float(wave)*10
        else:
            w = float(wave)

        wint = np.clip(np.interp(np.log10(w/800), [-0.3, 0.3], [0, 1]), 0, 1)

        print(filt, w, wint)

        clr = utils.RGBtoHex(plt.cm.Spectral_r(wint))
        #plt.scatter([0],[0], color=clr, label=filt)

        reg = prep.drizzle_footprint(weight_image, shrink=10, ext=0, outfile=None, label=filt) + ' color={0}\n'.format(clr)

        fp = open('{0}_mosaic.reg'.format(field_root), 'a')
        fp.write(reg)
        fp.close() 

def test_writeread_healpix_no_corrdsys(cst_power_1freq_cut_healpix, tmp_path):
    beam_in = cst_power_1freq_cut_healpix
    beam_out = UVBeam()

    write_file = str(tmp_path / "outtest_beam.fits")
    write_file2 = str(tmp_path / "outtest_beam2.fits")
    beam_in.write_beamfits(write_file, clobber=True)

    # now remove coordsys but leave ctype 1
    fname = fits.open(write_file)
    data = fname[0].data
    primary_hdr = fname[0].header
    primary_hdr.pop("COORDSYS")
    hdunames = uvutils._fits_indexhdus(fname)
    hpx_hdu = fname[hdunames["HPX_INDS"]]
    bandpass_hdu = fname[hdunames["BANDPARM"]]

    prihdu = fits.PrimaryHDU(data=data, header=primary_hdr)
    hdulist = fits.HDUList([prihdu, hpx_hdu, bandpass_hdu])

    hdulist.writeto(write_file2, overwrite=True)
    hdulist.close()

    beam_out.read_beamfits(write_file2)
    assert beam_in == beam_out

    return 

def test_read_cst_write_read_fits_change_freq_units(cst_power_1freq, tmp_path):
    # set up power beam
    beam_in = cst_power_1freq
    beam_out = UVBeam()

    write_file = str(tmp_path / "outtest_beam.fits")
    write_file2 = str(tmp_path / "outtest_beam2.fits")
    beam_in.write_beamfits(write_file, clobber=True)

    # now change frequency units
    fname = fits.open(write_file)
    data = fname[0].data
    primary_hdr = fname[0].header
    primary_hdr["CUNIT3"] = "MHz"
    primary_hdr["CRVAL3"] = primary_hdr["CRVAL3"] / 1e6
    primary_hdr["CDELT3"] = primary_hdr["CRVAL3"] / 1e6
    hdunames = uvutils._fits_indexhdus(fname)
    bandpass_hdu = fname[hdunames["BANDPARM"]]

    prihdu = fits.PrimaryHDU(data=data, header=primary_hdr)
    hdulist = fits.HDUList([prihdu, bandpass_hdu])

    hdulist.writeto(write_file2, overwrite=True)
    hdulist.close()

    beam_out.read_beamfits(write_file2)
    assert beam_in == beam_out

    return 

def test_read_cst_write_read_fits_no_coordsys(cst_power_1freq, tmp_path):
    # set up power beam
    beam_in = cst_power_1freq
    beam_out = UVBeam()

    write_file = str(tmp_path / "outtest_beam.fits")
    write_file2 = str(tmp_path / "outtest_beam2.fits")
    beam_in.write_beamfits(write_file, clobber=True)

    # now remove coordsys but leave ctypes 1 & 2
    fname = fits.open(write_file)
    data = fname[0].data
    primary_hdr = fname[0].header
    primary_hdr.pop("COORDSYS")
    hdunames = uvutils._fits_indexhdus(fname)
    bandpass_hdu = fname[hdunames["BANDPARM"]]

    prihdu = fits.PrimaryHDU(data=data, header=primary_hdr)
    hdulist = fits.HDUList([prihdu, bandpass_hdu])

    hdulist.writeto(write_file2, overwrite=True)
    hdulist.close()

    beam_out.read_beamfits(write_file2)
    assert beam_in == beam_out

    return 

def test_read_cst_write_read_fits_intensity(cst_power_1freq, tmp_path):
    # set up power beam
    beam_in = cst_power_1freq
    beam_out = UVBeam()

    write_file = str(tmp_path / "outtest_beam.fits")
    write_file2 = str(tmp_path / "outtest_beam2.fits")
    beam_in.write_beamfits(write_file, clobber=True)

    # now replace 'power' with 'intensity' for btype
    fname = fits.open(write_file)
    data = fname[0].data
    primary_hdr = fname[0].header
    primary_hdr["BTYPE"] = "Intensity"
    hdunames = uvutils._fits_indexhdus(fname)
    bandpass_hdu = fname[hdunames["BANDPARM"]]

    prihdu = fits.PrimaryHDU(data=data, header=primary_hdr)
    hdulist = fits.HDUList([prihdu, bandpass_hdu])

    hdulist.writeto(write_file2, overwrite=True)
    hdulist.close()

    beam_out.read_beamfits(write_file2)
    assert beam_in == beam_out

    return 

def scale_frm_wcs(fn):
    from astropy.io import fits
    hdu = fits.open(fn)
    head = hdu[0].header
    return scale_frm_header(head)