How to use the fitsio.FITS function in fitsio

'kernel' : kernel,
                'white_noise':1e-5,
                'optimize' : optimize
            }
        }
        logger = piff.config.setup_logger()
        interp3 = piff.Interp.process(config['interp'], logger)
        iterate(stars, interp3)
        validate(validate_stars, interp3)

    # Check that we can write interp to disk and read back in.
    if file_name is not None:
        testfile = os.path.join('output', file_name)
        with fitsio.FITS(testfile, 'rw', clobber=True) as f:
            interp.write(f, 'interp')
        with fitsio.FITS(testfile, 'r') as f:
            interp2 = piff.GPGeorgeInterp.read(f, 'interp')
        print("Revalidating after i/o.")
        X = np.vstack([training_data['u'], training_data['v']]).T
        for i in range(interp.nparams):
            print('A')
            np.testing.assert_allclose(interp.gps[i].get_matrix(X), interp2.gps[i].get_matrix(X))
            print('B')
            np.testing.assert_allclose(interp._init_theta[i], interp2._init_theta[i])
            print('C')
            np.testing.assert_allclose(interp.gps[i].get_parameter_vector(), interp2.gps[i].get_parameter_vector())
            #print('D')
            #print(interp.gps[i]._alpha)
            #print(interp2.gps[i]._alpha)
            #np.testing.assert_allclose(interp.gps[i]._alpha, interp2.gps[i]._alpha, rtol=1e-6, atol=1.e-7)
            print('E')
            np.testing.assert_allclose(interp.gps[i]._x, interp2.gps[i]._x)

'interp' : {
                'type' : 'GPGeorgeInterp',
                'kernel' : kernel,
                'white_noise':1e-5,
                'optimize' : optimize
            }
        }
        logger = piff.config.setup_logger()
        interp3 = piff.Interp.process(config['interp'], logger)
        iterate(stars, interp3)
        validate(validate_stars, interp3)

    # Check that we can write interp to disk and read back in.
    if file_name is not None:
        testfile = os.path.join('output', file_name)
        with fitsio.FITS(testfile, 'rw', clobber=True) as f:
            interp.write(f, 'interp')
        with fitsio.FITS(testfile, 'r') as f:
            interp2 = piff.GPGeorgeInterp.read(f, 'interp')
        print("Revalidating after i/o.")
        X = np.vstack([training_data['u'], training_data['v']]).T
        for i in range(interp.nparams):
            print('A')
            np.testing.assert_allclose(interp.gps[i].get_matrix(X), interp2.gps[i].get_matrix(X))
            print('B')
            np.testing.assert_allclose(interp._init_theta[i], interp2._init_theta[i])
            print('C')
            np.testing.assert_allclose(interp.gps[i].get_parameter_vector(), interp2.gps[i].get_parameter_vector())
            #print('D')
            #print(interp.gps[i]._alpha)
            #print(interp2.gps[i]._alpha)
            #np.testing.assert_allclose(interp.gps[i]._alpha, interp2.gps[i]._alpha, rtol=1e-6, atol=1.e-7)

def test_disk():
    # save and load
    r0 = 0.1
    sigma = 1.2
    g1 = -0.1
    g2 = 0.05
    gsparams = galsim.GSParams(minimum_fft_size=32)
    model = piff.Optical(r0=r0, sigma=sigma, g1=g1, g2=g2, lam=700.0, template='des',
                         pad_factor=0.5, oversampling=0.5, gsparams=gsparams)
    model_file = os.path.join('output','optics.fits')
    with fitsio.FITS(model_file, 'rw', clobber=True) as f:
        model.write(f, 'optics')
        model2 = piff.Optical.read(f, 'optics')

    for key in model.kwargs:
        assert key in model2.kwargs, 'key %r missing from model2 kwargs'%key
        assert model.kwargs[key] == model2.kwargs[key], 'key %r mismatch'%key
    for key in model.optical_psf_kwargs:
        assert key in model2.optical_psf_kwargs, 'key %r missing from model2 optical_psf_kwargs'%key
        assert model.optical_psf_kwargs[key] == model2.optical_psf_kwargs[key], 'key %r mismatch'%key
    for key in model.kolmogorov_kwargs:
        assert key in model2.kolmogorov_kwargs, 'key %r missing from model2 kolmogorov_kwargs'%key
        assert model.kolmogorov_kwargs[key] == model2.kolmogorov_kwargs[key], 'key %r mismatch'%key
    assert model.sigma == model2.sigma,'sigma mismatch'
    assert model.g1 == model2.g1,'g1 mismatch'
    assert model.g2 == model2.g2,'g2 mismatch'

T.height = []
    T.expnum = []
    T.extname = []
    T.calname = []
    T.seeing = []
    T.pixscale = []
    T.crval1 = []
    T.crval2 = []
    T.crpix1 = []
    T.crpix2 = []
    T.cd1_1 = []
    T.cd1_2 = []
    T.cd2_1 = []
    T.cd2_2 = []
    for i,(fn,expnum,seeing) in enumerate(zip(imfns, expnums, seeings)):
        F = fitsio.FITS(fn)
        primhdr = F[0].read_header()
        filter = primhdr['FILTER'].split('.')[0]
        exptime = primhdr['EXPTIME']
        pixscale = primhdr['PIXSCAL1'] / 3600.
        print('Pixscale:', pixscale * 3600, 'arcsec/pix')

        for hdu in range(1, len(F)):
        #for hdu in [13]:
            hdr = F[hdu].read_header()
            T.cpimage.append(fn)
            T.cpimage_hdu.append(hdu)
            T.filter.append(filter)
            T.exptime.append(exptime)
            args = []
            for k in ['CRVAL1', 'CRVAL2', 'CRPIX1', 'CRPIX2', 'CD1_1', 'CD1_2',
                      'CD2_1', 'CD2_2' ]:

raise

        self.outfilename = targetname + '_scan_' + str(self.cl.mapscans[0]) + '_' + str(self.cl.mapscans[-1]) + '_window' + str(window) + '_feed' + str(feed) + '_pol' + str(pol) + '.fits'

        self.log = Logging(self.cl, self.outfilename.rstrip('.fits'))

        if self.CLOBBER is False and os.path.exists(self.outfilename):
            self.log.doMessage('WARN', ' Will not overwrite existing pipeline output.\nConsider using \'--clobber\' option to overwrite.')
            sys.exit()

        # create a new table
        old_stdout = sys.stdout
        from cStringIO import StringIO
        sys.stdout = StringIO()
        # redirect stdout to not get clobber file warnings
        self.outfile = fitsio.FITS(self.outfilename, 'rw', clobber=True)
        sys.stdout = old_stdout

        dtype = self.infile[ext][0].dtype

        input_header = fitsio.read_header(self.infilename, ext)
        self.outfile.create_table_hdu(dtype=dtype, extname=input_header['EXTNAME'])
        self.outfile[-1].write_keys(input_header)
        self.outfile[-1]._update_info()
        self.outfile.update_hdu_list()

        # copy primary header from input to output file
        primary_header = fitsio.read_header(self.infilename, 0)
        self.outfile[0].write_keys(primary_header)

        self.outfile[0].write_key('PIPE_VER', PIPELINE_VERSION, comment="GBT Pipeline Version")

# Check for objects well off the edge.  We won't use them.
        big_bounds = image.bounds.expand(stamp_size)
        image_pos = [ pos for pos in image_pos if big_bounds.includes(pos) ]
        logger.debug("After remove those that are off the image, len = %s",len(image_pos))

        # Make the list of sky values:
        if sky_col is not None:
            if sky_col not in cat.dtype.names:
                raise ValueError("sky_col = %s is not a column in %s"%(sky_col,cat_file_name))
            sky = cat[sky_col]
        elif sky is not None:
            try:
                sky = float(sky)
            except ValueError:
                fits = fitsio.FITS(image_file_name)
                hdu = 1 if image_file_name.endswith('.fz') else 0
                header = fits[hdu].read_header()
                if sky not in header:
                    raise KeyError("Key %s not found in FITS header"%sky)
                sky = float(header[sky])
            sky = np.array([sky]*len(cat), dtype=float)
        else:
            sky = None

        # Make the list of gain values:
        # TODO: SV and Y1 DES images have two gain values, GAINA, GAINB.  It would be nice if we
        #       could get the right one properly.  OTOH, Y3+ will be in electrons, so gain=1 will
        #       the right value for all images.  So maybe not worth worrying about.
        if gain_col is not None:
            if gain_col not in cat.dtype.names:
                raise ValueError("gain_col = %s is not a column in %s"%(gain_col,cat_file_name))

def write_files(self, enhI, enhW, temp=False):
        if temp:
            mydir = self.get_dir()
            f,fn = tempfile.mkstemp(dir=mydir, suffix='.tmp')
            os.close(f)
        else:
            fn = self.get_imw_path()

        fits = fitsio.FITS(fn, 'rw', clobber=True)
        fits.write_image(enhI)
        fits.write_image(enhW)
        fits.close()
        return fn
        #imfn = self.get_image_path()

old = old + '.fz'
    print fn
    print old
    if fn == old:
        continue

    if fn in fits:
        F = fits[fn]
    else:
        pth = os.path.join(dirnm, fn.strip())
        if not os.path.exists(pth):
            p2 = pth.replace('.fits.fz', '.fits')
            if os.path.exists(p2):
                pth = p2
                #print 'Using', p2, 'instead of', pth
        F = fitsio.FITS(pth)
        fits[fn] = F

    if len(fits) > 100:
        fits = {}

    ff = F[C.cpimage_hdu[i]].read_header()
    ee = ff['EXTNAME'].strip()
    extname = C.extname[i].strip()
    if ee != extname:
        print 'HDU wrong for', fn
        if not (fn,extname) in hdus:
            for hdu in range(len(F)):
                hdr = F[hdu].read_header()
                if not 'EXTNAME' in hdr:
                    continue
                ext = hdr['EXTNAME'].strip()

def fits_table_driver(output, url, extension):
    """
    url -- fits file to open (local filename or URL)
    extension -- name or number of extension (HDU)
    """

    header = fitsio.read_header(url, ext=extension)
    n_entries = header["NAXIS2"]

    fitsfile = fitsio.FITS(url, iter_row_buffer=10000)
    table = fitsfile[extension]

    print("--INPUT --------------")
    print(fitsfile)
    print("")
    print(table)
    print("----------------------")

    for ii in range(n_entries):
        data = table[ii]
        output.send(dict(event=data))

hdu = 1 if image_file_name.endswith('.fz') else 0
                header = fits[hdu].read_header()
                if gain not in header:
                    raise KeyError("Key %s not found in FITS header"%gain)
                gain = float(header[gain])
            gain = np.array([gain]*len(cat), dtype=float)
        else:
            gain = [None] * len(cat)

        # Get the saturation level
        if satur is not None:
            try:
                satur = float(satur)
                logger.debug("Using given saturation value: %s",satur)
            except ValueError:
                fits = fitsio.FITS(image_file_name)
                hdu = 1 if image_file_name.endswith('.fz') else 0
                header = fits[hdu].read_header()
                if satur not in header:
                    raise KeyError("Key %s not found in FITS header"%satur)
                satur = float(header[satur])
                logger.debug("Using saturation from header: %s",satur)

        return image_pos, sky, gain, satur