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

if band_dict['use_error_map']==True:

        # Load in error map and extract WCS
        in_fits_error = astropy.io.fits.getdata(in_fitspath_error)

        # Slice smaller map out of original map
        out_fits_error = in_fits_error.copy()
        out_fits_error = out_fits_error[y_min_border:,:]
        out_fits_error = out_fits_error[:,x_min_border:]
        out_fits_error = out_fits_error[:-y_max_border,:]
        out_fits_error = out_fits_error[:,:-x_max_border]

        # Save cutout to file
        out_fitspath_error = os.path.join( kwargs_dict['temp_dir_path'], 'Cutouts', source_dict['name'], source_dict['name']+'_'+band_dict['band_name']+'_Error.fits' )
        out_cutout_hdu_error = astropy.io.fits.PrimaryHDU(data=out_fits_error, header=out_header)
        out_cutout_hdulist_error = astropy.io.fits.HDUList([out_cutout_hdu_error])
        out_cutout_hdulist_error.writeto(out_fitspath_error, clobber=True)

    # Return the directory of the newly-created cutout
    out_fitsdir = os.path.split(out_fitspath)[0]
    band_dict['band_dir'] = out_fitsdir
    return band_dict

hdu0.header['JMAG'] = ('Unknown', '[mag] J band magnitude from 2MASS')
            hdu0.header['HMAG'] = ('Unknown', '[mag] H band magnitude from 2MASS')
            hdu0.header['KMAG'] = ('Unknown', '[mag] K band magnitude from 2MASS')
            hdu0.header['KEPMAG'] = ('Unknown', '[mag] Kepler magnitude (Kp) from KIC')
            hdu0.header['GRCOLOR'] = ('Unknown', '[mag] (g-r) color, SDSS bands')
            hdu0.header['JKCOLOR'] = ('Unknown', '[mag] (J-K) color, 2MASS bands')
            hdu0.header['GKCOLOR'] = ('Unknown', '[mag] (g-K) color, SDSS g - 2MASS K')
            hdu0.header['TEFF'] = ('Unknown', '[K] effective temperature from KIC')
            hdu0.header['LOGG'] = ('Unknown', '[cm/s2] log10 surface gravity from KIC')
            hdu0.header['FEH'] = ('Unknown', '[log10([Fe/H])] metallicity from KIC')
            hdu0.header['EBMINUSV'] = ('Unknown', '[mag] E(B-V) redenning from KIC')
            hdu0.header['AV'] = ('Unknown', '[mag] A_v extinction from KIC')
            hdu0.header['RADIUS'] = ('Unknown', '[solar radii] stellar radius from KIC')
            hdu0.header['TMINDEX'] = ('Unknown', 'unique 2MASS catalog ID from KIC')
            hdu0.header['SCPID'] = ('Unknown', 'unique SCP processing ID from KIC')
            hdulist = pyfits.HDUList(hdu0)
        except:
            errmsg = ('ERROR -- KEPCONVERT: cannot create primary extension in {}'
                      .format(outfile))
            kepmsg.err(logfile, errmsg, verbose)
        ## create the outfile HDU 1 extension
        try:
            fitscol = []
            for i in range(ncol):
                fitscol.append(pytfits.Column(name=colnames[i], format='D',
                                              array=work[i]))
            fitscols = pyfits.ColDefs(fitscol)
            hdu1 = pyfits.BinTableHDU.from_columns(fitscols)
            hdulist.append(hdu1)
            hdu1.header['INHERIT'] = (True, 'inherit primary keywords')
            hdu1.header['EXTNAME'] = ('LIGHTCURVE', 'name of extension')
            hdu1.header['EXTVER'] = (1, 'extension version number')

def _make_hdulist(self):
        hdu0 = fits.PrimaryHDU()

        # Set the header with defaults
        tmpl = self.get_header_template(0)
        for kw in tmpl:
            hdu0.header[kw] = (tmpl[kw], tmpl.comments[kw])
        # Override the primary extension defaults
        hdu0.header['ORIGIN'] = "Unofficial data product"
        hdu0.header['DATE'] = datetime.datetime.now().strftime("%Y-%m-%d")
        hdu0.header['CREATOR'] = "kadenza"
        hdu0.header['PROCVER'] = "{}".format(__version__)

        hdulist = fits.HDUList(hdu0)

        # open the Cadence Pixel File and Pixel Mapping File
        incpf = fits.open(self.cadence_pixel_file, memmap=True)
        inpmrf = fits.open(self.pixel_mapping_file, memmap=True)

        # grab the dates and times from the header
        card = incpf[0].header
        startkey = card['STARTIME']
        endkey = card['END_TIME']
        dateobs = card['DATE-OBS']
        timeobs = card['TIME-OBS']
        tstart = startkey + 2400000.5
        tstop = endkey + 2400000.5
        deadc = 0.92063492
        telapse = tstop - tstart
        exposure = telapse * deadc

Returns
    -------
    is_fits : bool
        Returns `True` if the given file is a FITS file.
    """
    if fileobj is not None:
        pos = fileobj.tell()
        sig = fileobj.read(30)
        fileobj.seek(pos)
        return sig == FITS_SIGNATURE
    elif filepath is not None:
        if filepath.lower().endswith(('.fits', '.fits.gz', '.fit', '.fit.gz',
                                      '.fts', '.fts.gz')):
            return True
    elif isinstance(args[0], (fits.HDUList, fits.ImageHDU, fits.PrimaryHDU)):
        return True
    else:
        return False

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 _hdulist(self, hdu0_keywords, ext_info):
        """Returns an astropy.io.fits.HDUList object."""
        return fits.HDUList([self._make_primary_hdu(hdu0_keywords=hdu0_keywords),
                             self._make_target_extension(ext_info=ext_info),
                             self._make_aperture_extension()])

ax.plot(xx, collapse, '.', label='collapsed')
            ax.plot(xx, collapse_smooth, '-', label='savgol{}'.format(degree))
            ax.plot(xx, collapse_smooth_s, '--', label='spline{}'.format(degree_s))
            ax.legend()
            plt.savefig('collapsed_smooth.png')
            plt.close()

        # Divide each fiber in rss_wl by spectrum
        gmean = col_good_mean.mean()
        data1 = rss_wl[0].data / collapse_smooth_s
        data1 /= gmean
        # Fill values with ones to avoid NaNs
        data2 = numpy.where(wlmap > 0, data1, 1.0)

        self.logger.warning("Copy all extensions for the moment")
        rss_wl2 = fits.HDUList([hdu.copy() for hdu in rss_wl])
        rss_wl2[0].data = data2
        return rss_wl2

"""Make images that can be used to create profiles.

    The following images can be created:
    * LON -- Longitude coordinate
    * LAT -- Latitude coordinate
    * DIST -- Distance to mask
    * SOLID_ANGLE -- Solid angle
    """
    from astropy.io import fits
    from gammapy.utils.fits import get_hdu
    from gammapy.image import SkyImage, SkyMask

    log.info('Reading {0}'.format(infile))
    hdu = get_hdu(infile)

    out_hdus = fits.HDUList()

    if make_coordinate_maps:
        image = SkyImage.empty_like(hdu)
        log.info('Computing LON and LAT images')
        lon, lat = image.coordinates()
        out_hdus.append(fits.ImageHDU(lon, hdu.header, 'LON'))
        out_hdus.append(fits.ImageHDU(lat, hdu.header, 'LAT'))

    if make_distance_map:
        excl = SkyMask.from_image_hdu(hdu)
        log.info('Computing DIST map')
        dist = excl.exclusion_distance
        out_hdus.append(fits.ImageHDU(dist, hdu.header, 'DIST'))

    log.info('Writing {0}'.format(outfile))
    out_hdus.writeto(outfile, clobber=overwrite)

    @qglue_parser(HDUList, priority=100)
    def _parse_data_hdulist(data, label):
        from glue.core.data_factories.fits import fits_reader
        return fits_reader(data, label=label)

for h in hdu_i:
                h.header['EXTVER'] = filt_i
                if hdul is None:
                    hdul = pyfits.HDUList([h])
                else:
                    hdul.append(h)

            print('Add to {0}.thumb.fits: {1}'.format(label, file))

            # Weight
            hdu_i = pyfits.open(file.replace('_sci', '_wht'))
            hdu_i[0].header['EXTNAME'] = 'WHT'
            for h in hdu_i:
                h.header['EXTVER'] = filt_i
                if hdul is None:
                    hdul = pyfits.HDUList([h])
                else:
                    hdul.append(h)

        hdul.writeto('{0}.thumb.fits'.format(label), overwrite=True,
                     output_verify='fix')

        for file in files:
            for f in [file, file.replace('_sci', '_wht')]:
                if os.path.exists(f):
                    print('Remove {0}'.format(f))
                    os.remove(f)

    # Segmentation figure
    thumb_file = '{0}.thumb.fits'.format(label)
    if (make_segmentation_figure) & (os.path.exists(thumb_file)) & (aws_prep_dir is not None):