How to use the astropy.io.fits.open function in astropy

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)

info['sSFR'].format = '.1e'
    info['stellar_mass'].format = '.1e'
    
    
    if filter_bandpasses:
        arr = np.array(template_mags)
        for i, bp in enumerate(filter_bandpasses):
            info['mag_{0}'.format(bp.name)] = arr[:,i]
            info['mag_{0}'.format(bp.name)].format = '.3f'
            
    for c in info.colnames:
        info.rename_column(c, c.lower())
    
    # Emission line names
    #files=glob.glob('{0}*full.fits'.format(target))
    im = pyfits.open(files[0])
    h = im['COVAR'].header
    for i in range(24):
        key = 'FLUX_{0:03d}'.format(i)
        if key not in h:
            continue
        
        line = h.comments[key].split()[0]
        
        for root in ['flux','err','ew50','ewhw']:
            col = '{0}_{1}'.format(root, line)
            info.rename_column('{0}_{1:03d}'.format(root, i), col)
            if root.startswith('ew'):
                info[col].format = '.1f'
            else:
                info[col].format = '.1f'

def apStar_loader(file_obj, **kwargs):
    """
    Loader for APOGEE apStar files.

    Parameters
    ----------
    file_obj: str or file-like
        FITS file name or object (provided from name by Astropy I/O Registry).

    Returns
    -------
    data: Spectrum1D
        The spectrum that is represented by the data in this table.
    """
    hdulist = fits.open(file_obj, **kwargs)

    header = hdulist[0].header
    meta = {'header': header}
    wcs = WCS(hdulist[1].header)

    data = hdulist[1].data[0, :]  # spectrum in the first row of the first extension
    unit = Unit('1e-17 erg / (Angstrom cm2 s)')

    uncertainty = StdDevUncertainty(hdulist[2].data[0, :])

    # dispersion from the WCS but convert out of logspace
    # dispersion = 10**wcs.all_pix2world(np.arange(data.shape[0]), 0)[0]
    dispersion = 10**wcs.all_pix2world(np.vstack((np.arange(data.shape[0]),
                                                  np.zeros((data.shape[0],)))).T,
                                       0)[:, 0]
    dispersion_unit = Unit('Angstrom')

from astropy.io import fits as pyfits

import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.axes_grid1.axes_grid import AxesGrid
#from pywcsgrid2.axes_wcs import GridHelperWcs, AxesWcs
import pywcsgrid2

# read in the first image
xray_name="pspc_skyview.fits"
f_xray = pyfits.open(xray_name)
header_xray = f_xray[0].header

# the second image
radio_name="radio_21cm.fits"
f_radio = pyfits.open(radio_name)
header_radio = f_radio[0].header


# grid helper
grid_helper = pywcsgrid2.GridHelper(wcs=header_xray)

# AxesGrid to display tow images side-by-side
fig = plt.figure(1, (6,3.5))
grid = AxesGrid(fig, (0.15, 0.15, 0.8, 0.75), nrows_ncols=(1, 2),
                axes_pad=0.1, share_all=True,
                cbar_mode="each", cbar_location="top", cbar_pad=0,

def make_scaled_srcmap(self):
        """Make an exposure cube with the same binning as the counts map."""

        self.logger.info('Computing scaled source map.')

        bexp0 = pyfits.open(self._bexpmap_roi_file)
        bexp1 = pyfits.open(self.config['gtlike']['bexpmap'])
        srcmap = pyfits.open(self.config['gtlike']['srcmap'])

        if bexp0[0].data.shape != bexp1[0].data.shape:
            raise Exception('Wrong shape for input exposure map file.')

        bexp_ratio = bexp0[0].data / bexp1[0].data

        self.logger.info(
            'Min/Med/Max exposure correction: %f %f %f' % (np.min(bexp_ratio),
                                                           np.median(
                                                               bexp_ratio),
                                                           np.max(bexp_ratio)))

        for hdu in srcmap[1:]:

            if hdu.name == 'GTI': continue

def cos_identify(*args, **kwargs):
    """Check whether given file contains HST/COS spectral data."""
    with fits.open(args[0]) as hdu:
        if hdu[0].header['TELESCOP'] == 'HST' and hdu[0].header['INSTRUME'] == 'COS':
           return True

    return False

def get_fits_header(filename):
    """
    Reads a FITS header. A wrapper for common astropy.io.fits commands.

    :param filename: The name of the FITS file to retrieve the header from.

    :type filename: str

    :returns: Header instance -- The header from the primary HDU.
    """

    hdulist = pyfits.open(filename, memmap=1)

    htab = hdulist[0].header

    hdulist.close()

    return htab
# ------------------------------------------------------------------------------

def identify_wcs1d_fits(origin, *args, **kwargs):
    # check if file can be opened with this reader
    # args[0] = filename
    with fits.open(args[0]) as hdu:
        # check if number of axes is one
        return (hdu[0].header['NAXIS'] == 1 and
                hdu[0].header.get('WCSDIM', 1) == 1 and
                # check in CTYPE1 key for linear solution
                hdu[0].header.get('CTYPE1', '').upper() != 'MULTISPEC')

    return False

import matplotlib.pyplot as plt
from astropy.io import fits as pyfits
from mpl_toolkits.axes_grid1.axes_grid import ImageGrid
import pywcsgrid2
#from pywcsgrid2.axes_wcs import AxesWcs, GridHelperWcs
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset
import mpl_toolkits.axisartist as axisartist
import matplotlib.patheffects as patheffects

if 1:

    f = pyfits.open("pspc_skyview.fits")
    d = f[0].data
    h = f[0].header

    fig = plt.figure(1)

    grid = ImageGrid(fig, (1, 1, 1), nrows_ncols = (1, 1),
                     cbar_mode="single", cbar_pad="2%",
                     cbar_location="right",
                     axes_class=(pywcsgrid2.Axes, dict(header=h)))

    main_axes = grid[0]
    main_axes.locator_params(nbins=5)

    cb_axes = grid.cbar_axes[0] # colorbar axes

    im = main_axes.imshow(d, origin="lower", cmap=plt.cm.gray_r,

vol=self.sregion.shape
        for v in range(0,vol[2]-1):
            colored[:,:,v] = self.extent[4] + v*(-1)*abs(self.hdr['cdelt3'])
        field.image_data.point_data.add_array(colored.T.ravel())
        field.image_data.point_data.get_array(1).name = 'color'
        field.update()
        
        field2 = mlab.pipeline.set_active_attribute(field, point_scalars='scalar')
        contour = mlab.pipeline.contour(field2)
        contour2 = mlab.pipeline.set_active_attribute(contour, point_scalars='color')
        
        mlab.pipeline.surface(contour2, colormap='jet', opacity=self.opacity)
        
        ## Insert a continuum plot
        if self.contfile != '':
            im = fits.open(self.contfile)
            dat = im[0].data
            ##dat0 = dat[0]
            channel = dat[0]
            region = np.swapaxes(channel[self.ystart:self.yend,self.xstart:self.xend]*1000.,0,1)
            field = mlab.contour3d(region, colormap='gist_ncar')
            field.contour.minimum_contour = 5
        
        self.field = field2
        self.field.scene.render()
        self.labels()
        mlab.view(azimuth=0, elevation=0, distance='auto')
        mlab.show()