def GetData(EPIC, season=None, cadence='lc', clobber=False, delete_raw=False,
            aperture_name='k2sff_15', saturated_aperture_name='k2sff_19',
            max_pixels=75, download_only=False, saturation_tolerance=-0.1,
            bad_bits=[1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 16, 17],
            get_hires=True,
            get_nearby=True, **kwargs):
    '''
    Returns a :py:obj:`DataContainer` instance with the
    raw data for the target.

    :param int EPIC: The EPIC ID number
    :param int season: The observing season (campaign). Default :py:obj:`None`
    :param str cadence: The light curve cadence. Default `lc`
    :param bool clobber: Overwrite existing files? Default :py:obj:`False`
    :param bool delete_raw: Delete the FITS TPF after processing it? \
           Default :py:obj:`False`
    :param str aperture_name: The name of the aperture to use. Select \
           `custom` to call :py:func:`GetCustomAperture`. Default `k2sff_15`
    :param str saturated_aperture_name: The name of the aperture to use if \
           the target is saturated. Default `k2sff_19`
    :param int max_pixels: Maximum number of pixels in the TPF. Default 75
    :param bool download_only: Download raw TPF and return? Default \
           :py:obj:`False`
    :param float saturation_tolerance: Target is considered saturated \
           if flux is within this fraction of the pixel well depth. \
           Default -0.1
    :param array_like bad_bits: Flagged :py:obj`QUALITY` bits to consider \
           outliers when computing the model. \
           Default `[1,2,3,4,5,6,7,8,9,11,12,13,14,16,17]`
    :param bool get_hires: Download a high resolution image of the target? \
           Default :py:obj:`True`
    :param bool get_nearby: Retrieve location of nearby sources? \
           Default :py:obj:`True`

    '''

    # Campaign no.
    if season is None:
        campaign = Season(EPIC)
        if hasattr(campaign, '__len__'):
            raise AttributeError(
                "Please choose a campaign/season for this target: %s."
                % campaign)
    else:
        campaign = season

    # Is there short cadence data available for this target?
    short_cadence = HasShortCadence(EPIC, season=campaign)
    if cadence == 'sc' and not short_cadence:
        raise ValueError("Short cadence data not available for this target.")

    # Local file name
    filename = os.path.join(EVEREST_DAT, 'k2', 'c%02d' % campaign,
                            ('%09d' % EPIC)[:4] + '00000', ('%09d' % EPIC)[4:],
                            'data.npz')

    # Download?
    if clobber or not os.path.exists(filename):

        # Get the TPF
        tpf = os.path.join(KPLR_ROOT, 'data', 'k2', 'target_pixel_files',
                           str(EPIC), 'ktwo%09d-c%02d_lpd-targ.fits.gz'
                           % (EPIC, campaign))
        sc_tpf = os.path.join(KPLR_ROOT, 'data', 'k2', 'target_pixel_files',
                              str(EPIC), 'ktwo%09d-c%02d_spd-targ.fits.gz'
                              % (EPIC, campaign))
        if clobber or not os.path.exists(tpf):
            kplr_client.k2_star(EPIC).get_target_pixel_files(fetch=True)

        with pyfits.open(tpf) as f:
            qdata = f[1].data

            # Get the TPF aperture
            tpf_aperture = (f[2].data & 2) // 2

            # Get the enlarged TPF aperture
            tpf_big_aperture = np.array(tpf_aperture)
            for i in range(tpf_big_aperture.shape[0]):
                for j in range(tpf_big_aperture.shape[1]):
                    if f[2].data[i][j] == 1:
                        for n in [(i - 1, j), (i + 1, j),
                                  (i, j - 1), (i, j + 1)]:
                            if n[0] >= 0 and n[0] < tpf_big_aperture.shape[0]:
                                if n[1] >= 0 and n[1] < \
                                        tpf_big_aperture.shape[1]:
                                    if tpf_aperture[n[0]][n[1]] == 1:
                                        tpf_big_aperture[i][j] = 1

        # Is there short cadence data?
        if short_cadence:
            with pyfits.open(sc_tpf) as f:
                sc_qdata = f[1].data

        # Get K2SFF apertures
        try:
            k2sff = kplr.K2SFF(EPIC, sci_campaign=campaign)
            k2sff_apertures = k2sff.apertures
            if delete_raw:
                os.remove(k2sff._file)
        except:
            k2sff_apertures = [None for i in range(20)]

        # Make a dict of all our apertures
        # We're not getting K2SFF apertures 0-9 any more
        apertures = {'tpf': tpf_aperture, 'tpf_big': tpf_big_aperture}
        for i in range(10, 20):
            apertures.update({'k2sff_%02d' % i: k2sff_apertures[i]})

        # Get the header info
        fitsheader = [pyfits.getheader(tpf, 0).cards,
                      pyfits.getheader(tpf, 1).cards,
                      pyfits.getheader(tpf, 2).cards]
        if short_cadence:
            sc_fitsheader = [pyfits.getheader(sc_tpf, 0).cards,
                             pyfits.getheader(sc_tpf, 1).cards,
                             pyfits.getheader(sc_tpf, 2).cards]
        else:
            sc_fitsheader = None

        # Get a hi res image of the target
        if get_hires:
            hires = GetHiResImage(EPIC)
        else:
            hires = None

        # Get nearby sources
        if get_nearby:
            nearby = GetSources(EPIC)
        else:
            nearby = []

        # Delete?
        if delete_raw:
            os.remove(tpf)
            if short_cadence:
                os.remove(sc_tpf)

        # Get the arrays
        cadn = np.array(qdata.field('CADENCENO'), dtype='int32')
        time = np.array(qdata.field('TIME'), dtype='float64')
        fpix = np.array(qdata.field('FLUX'), dtype='float64')
        fpix_err = np.array(qdata.field('FLUX_ERR'), dtype='float64')
        qual = np.array(qdata.field('QUALITY'), dtype=int)

        # Get rid of NaNs in the time array by interpolating
        naninds = np.where(np.isnan(time))
        time = Interpolate(np.arange(0, len(time)), naninds, time)

        # Get the motion vectors (if available!)
        pc1 = np.array(qdata.field('POS_CORR1'), dtype='float64')
        pc2 = np.array(qdata.field('POS_CORR2'), dtype='float64')
        if not np.all(np.isnan(pc1)) and not np.all(np.isnan(pc2)):
            pc1 = Interpolate(time, np.where(np.isnan(pc1)), pc1)
            pc2 = Interpolate(time, np.where(np.isnan(pc2)), pc2)
        else:
            pc1 = None
            pc2 = None

        # Do the same for short cadence
        if short_cadence:
            sc_cadn = np.array(sc_qdata.field('CADENCENO'), dtype='int32')
            sc_time = np.array(sc_qdata.field('TIME'), dtype='float64')
            sc_fpix = np.array(sc_qdata.field('FLUX'), dtype='float64')
            sc_fpix_err = np.array(sc_qdata.field('FLUX_ERR'), dtype='float64')
            sc_qual = np.array(sc_qdata.field('QUALITY'), dtype=int)
            sc_naninds = np.where(np.isnan(sc_time))
            sc_time = Interpolate(
                np.arange(0, len(sc_time)), sc_naninds, sc_time)
            sc_pc1 = np.array(sc_qdata.field('POS_CORR1'), dtype='float64')
            sc_pc2 = np.array(sc_qdata.field('POS_CORR2'), dtype='float64')
            if not np.all(np.isnan(sc_pc1)) and not np.all(np.isnan(sc_pc2)):
                sc_pc1 = Interpolate(
                    sc_time, np.where(np.isnan(sc_pc1)), sc_pc1)
                sc_pc2 = Interpolate(
                    sc_time, np.where(np.isnan(sc_pc2)), sc_pc2)
            else:
                sc_pc1 = None
                sc_pc2 = None
        else:
            sc_cadn = None
            sc_time = None
            sc_fpix = None
            sc_fpix_err = None
            sc_qual = None
            sc_pc1 = None
            sc_pc2 = None

        # Static pixel images for plotting
        pixel_images = [fpix[0], fpix[len(fpix) // 2], fpix[len(fpix) - 1]]

        # Atomically write to disk.
        # http://stackoverflow.com/questions/2333872/
        # atomic-writing-to-file-with-python
        if not os.path.exists(os.path.dirname(filename)):
            os.makedirs(os.path.dirname(filename))
        f = NamedTemporaryFile("wb", delete=False)
        np.savez_compressed(f, cadn=cadn, time=time, fpix=fpix,
                            fpix_err=fpix_err,
                            qual=qual, apertures=apertures,
                            pc1=pc1, pc2=pc2, fitsheader=fitsheader,
                            pixel_images=pixel_images, nearby=nearby,
                            hires=hires,
                            sc_cadn=sc_cadn, sc_time=sc_time, sc_fpix=sc_fpix,
                            sc_fpix_err=sc_fpix_err, sc_qual=sc_qual,
                            sc_pc1=sc_pc1, sc_pc2=sc_pc2,
                            sc_fitsheader=sc_fitsheader)
        f.flush()
        os.fsync(f.fileno())
        f.close()
        shutil.move(f.name, filename)

        if download_only:
            return

    # Load
    data = np.load(filename)
    apertures = data['apertures'][()]
    pixel_images = data['pixel_images']
    nearby = data['nearby']
    hires = data['hires'][()]

    if cadence == 'lc':
        fitsheader = data['fitsheader']
        cadn = data['cadn']
        time = data['time']
        fpix = data['fpix']
        fpix_err = data['fpix_err']
        qual = data['qual']
        pc1 = data['pc1']
        pc2 = data['pc2']
    elif cadence == 'sc':
        fitsheader = data['sc_fitsheader']
        cadn = data['sc_cadn']
        time = data['sc_time']
        fpix = data['sc_fpix']
        fpix_err = data['sc_fpix_err']
        qual = data['sc_qual']
        pc1 = data['sc_pc1']
        pc2 = data['sc_pc2']
    else:
        raise ValueError("Invalid value for the cadence.")

    # Select the "saturated aperture" to check if the star is saturated
    # If it is, we will use this aperture instead
    if saturated_aperture_name == 'custom':
        saturated_aperture = GetCustomAperture(data)
    else:
        if saturated_aperture_name is None:
            saturated_aperture_name = 'k2sff_19'
        saturated_aperture = apertures[saturated_aperture_name]
        if saturated_aperture is None:
            log.error("Invalid aperture selected. Defaulting to `tpf_big`.")
            saturated_aperture_name = 'tpf_big'
            saturated_aperture = apertures[saturated_aperture_name]

    # HACK: Some C05 K2SFF apertures don't match the target pixel file
    # pixel grid size. This is likely because they're defined on the M67
    # superstamp. For now, let's ignore these stars.
    if saturated_aperture.shape != fpix.shape[1:]:
        log.error("Aperture size mismatch!")
        return None

    # Compute the saturation flux and the 97.5th percentile
    # flux in each pixel of the saturated aperture. We're going
    # to compare these to decide if the star is saturated.
    satflx = SaturationFlux(EPIC, campaign=campaign) * \
        (1. + saturation_tolerance)
    f97 = np.zeros((fpix.shape[1], fpix.shape[2]))
    for i in range(fpix.shape[1]):
        for j in range(fpix.shape[2]):
            if saturated_aperture[i, j]:
                # Let's remove NaNs...
                tmp = np.delete(fpix[:, i, j], np.where(
                    np.isnan(fpix[:, i, j])))
                # ... and really bad outliers...
                if len(tmp):
                    f = SavGol(tmp)
                    med = np.nanmedian(f)
                    MAD = 1.4826 * np.nanmedian(np.abs(f - med))
                    bad = np.where((f > med + 10. * MAD) |
                                   (f < med - 10. * MAD))[0]
                    np.delete(tmp, bad)
                    # ... so we can compute the 97.5th percentile flux
                    i97 = int(0.975 * len(tmp))
                    tmp = tmp[np.argsort(tmp)[i97]]
                    f97[i, j] = tmp

    # Check if any of the pixels are actually saturated
    if np.nanmax(f97) <= satflx:
        log.info("No saturated columns detected.")
        saturated = False
    else:
        log.info("Saturated pixel(s) found. Switching to aperture `%s`." %
                 saturated_aperture_name)
        aperture_name = saturated_aperture_name
        saturated = True

    # Now grab the aperture we'll actually use
    if aperture_name == 'custom':
        aperture = GetCustomAperture(data)
    else:
        if aperture_name is None:
            aperture_name = 'k2sff_15'
        aperture = apertures[aperture_name]
        if aperture is None:
            log.error("Invalid aperture selected. Defaulting to `tpf_big`.")
            aperture_name = 'tpf_big'
            aperture = apertures[aperture_name]

    # HACK: Some C05 K2SFF apertures don't match the target pixel file
    # pixel grid size. This is likely because they're defined on the M67
    # superstamp. For now, let's ignore these stars.
    if aperture.shape != fpix.shape[1:]:
        log.error("Aperture size mismatch!")
        return None

    # Now we check if the aperture is too big. Can lead to memory errors...
    # Treat saturated and unsaturated stars differently.
    if saturated:

        # Need to check if we have too many pixels *after* collapsing columns.
        # Sort the apertures in decreasing order of pixels, but keep the apert.
        # chosen by the user first.
        aperture_names = np.array(list(apertures.keys()))
        npix_per_aperture = np.array(
            [np.sum(apertures[k]) for k in aperture_names])
        aperture_names = aperture_names[np.argsort(npix_per_aperture)[::-1]]
        aperture_names = np.append([aperture_name], np.delete(
            aperture_names, np.argmax(aperture_names == aperture_name)))

        # Loop through them. Pick the first one that satisfies
        # the `max_pixels` constraint
        for aperture_name in aperture_names:
            aperture = apertures[aperture_name]
            aperture[np.isnan(fpix[0])] = 0
            ncol = 0
            apcopy = np.array(aperture)
            for j in range(apcopy.shape[1]):
                if np.any(f97[:, j] > satflx):
                    apcopy[:, j] = 0
                    ncol += 1
            if np.sum(apcopy) + ncol <= max_pixels:
                break
        if np.sum(apcopy) + ncol > max_pixels:
            log.error(
                "No apertures available with fewer than %d pixels. Aborting."
                % max_pixels)
            return None

        # Now, finally, we collapse the saturated columns into single pixels
        # and make the pixel array 2D
        ncol = 0
        fpixnew = []
        ferrnew = []

        # HACK: K2SFF sometimes clips the heads/tails of saturated columns
        # That's really bad, since that's where all the information is. Let's
        # artificially extend the aperture by two pixels at the top and bottom
        # of each saturated column. This *could* increase contamination, but
        # it's unlikely since the saturated target is by definition really
        # bright
        ext = 0
        for j in range(aperture.shape[1]):
            if np.any(f97[:, j] > satflx):
                for i in range(aperture.shape[0]):
                    if (aperture[i, j] == 0) and \
                            (np.nanmedian(fpix[:, i, j]) > 0):
                        if (i + 2 < aperture.shape[0]) and \
                                aperture[i + 2, j] == 1:
                            aperture[i, j] = 2
                            ext += 1
                        elif (i + 1 < aperture.shape[0]) and \
                                aperture[i + 1, j] == 1:
                            aperture[i, j] = 2
                            ext += 1
                        elif (i - 1 >= 0) and aperture[i - 1, j] == 1:
                            aperture[i, j] = 2
                            ext += 1
                        elif (i - 2 >= 0) and aperture[i - 2, j] == 1:
                            aperture[i, j] = 2
                            ext += 1
        if ext:
            log.info("Extended saturated columns by %d pixel(s)." % ext)

        for j in range(aperture.shape[1]):
            if np.any(f97[:, j] > satflx):
                marked = False
                collapsed = np.zeros(len(fpix[:, 0, 0]))
                collapsed_err2 = np.zeros(len(fpix[:, 0, 0]))
                for i in range(aperture.shape[0]):
                    if aperture[i, j]:
                        if not marked:
                            aperture[i, j] = AP_COLLAPSED_PIXEL
                            marked = True
                        else:
                            aperture[i, j] = AP_SATURATED_PIXEL
                        collapsed += fpix[:, i, j]
                        collapsed_err2 += fpix_err[:, i, j] ** 2
                if np.any(collapsed):
                    fpixnew.append(collapsed)
                    ferrnew.append(np.sqrt(collapsed_err2))
                    ncol += 1
            else:
                for i in range(aperture.shape[0]):
                    if aperture[i, j]:
                        fpixnew.append(fpix[:, i, j])
                        ferrnew.append(fpix_err[:, i, j])
        fpix2D = np.array(fpixnew).T
        fpix_err2D = np.array(ferrnew).T
        log.info("Collapsed %d saturated column(s)." % ncol)

    else:

        # Check if there are too many pixels
        if np.sum(aperture) > max_pixels:

            # This case is simpler: we just pick the largest aperture
            # that's less than or equal to `max_pixels`
            keys = list(apertures.keys())
            npix = np.array([np.sum(apertures[k]) for k in keys])
            aperture_name = keys[np.argmax(npix * (npix <= max_pixels))]
            aperture = apertures[aperture_name]
            aperture[np.isnan(fpix[0])] = 0
            if np.sum(aperture) > max_pixels:
                log.error("No apertures available with fewer than " +
                          "%d pixels. Aborting." % max_pixels)
                return None
            log.warn(
                "Selected aperture is too big. Proceeding with aperture " +
                "`%s` instead." % aperture_name)

        # Make the pixel flux array 2D
        aperture[np.isnan(fpix[0])] = 0
        ap = np.where(aperture & 1)
        fpix2D = np.array([f[ap] for f in fpix], dtype='float64')
        fpix_err2D = np.array([p[ap] for p in fpix_err], dtype='float64')

    # Compute the background
    binds = np.where(aperture ^ 1)
    if RemoveBackground(EPIC, campaign=campaign) and (len(binds[0]) > 0):
        bkg = np.nanmedian(np.array([f[binds]
                                     for f in fpix], dtype='float64'), axis=1)
        # Uncertainty of the median:
        # http://davidmlane.com/hyperstat/A106993.html
        bkg_err = 1.253 * np.nanmedian(np.array([e[binds] for e in fpix_err],
                                       dtype='float64'), axis=1) \
            / np.sqrt(len(binds[0]))
        bkg = bkg.reshape(-1, 1)
        bkg_err = bkg_err.reshape(-1, 1)
    else:
        bkg = 0.
        bkg_err = 0.

    # Make everything 2D and remove the background
    fpix = fpix2D - bkg
    fpix_err = np.sqrt(fpix_err2D ** 2 + bkg_err ** 2)
    flux = np.sum(fpix, axis=1)
    ferr = np.sqrt(np.sum(fpix_err ** 2, axis=1))

    # Get NaN data points
    nanmask = np.where(np.isnan(flux) | (flux == 0))[0]

    # Get flagged data points -- we won't train our model on them
    badmask = []
    for b in bad_bits:
        badmask += list(np.where(qual & 2 ** (b - 1))[0])

    # Flag >10 sigma outliers -- same thing.
    tmpmask = np.array(list(set(np.concatenate([badmask, nanmask]))))
    t = np.delete(time, tmpmask)
    f = np.delete(flux, tmpmask)
    f = SavGol(f)
    med = np.nanmedian(f)
    MAD = 1.4826 * np.nanmedian(np.abs(f - med))
    bad = np.where((f > med + 10. * MAD) | (f < med - 10. * MAD))[0]
    badmask.extend([np.argmax(time == t[i]) for i in bad])

    # Campaign 2 hack: the first day or two are screwed up
    if campaign == 2:
        badmask.extend(np.where(time < 2061.5)[0])

    # TODO: Fix time offsets in first half of
    # Campaign 0. See note in everest 1.0 code

    # Finalize the mask
    badmask = np.array(sorted(list(set(badmask))))

    # Interpolate the nans
    fpix = Interpolate(time, nanmask, fpix)
    fpix_err = Interpolate(time, nanmask, fpix_err)

    # Return
    data = DataContainer()
    data.ID = EPIC
    data.campaign = campaign
    data.cadn = cadn
    data.time = time
    data.fpix = fpix
    data.fpix_err = fpix_err
    data.nanmask = nanmask
    data.badmask = badmask
    data.aperture = aperture
    data.aperture_name = aperture_name
    data.apertures = apertures
    data.quality = qual
    data.Xpos = pc1
    data.Ypos = pc2
    data.meta = fitsheader
    data.mag = fitsheader[0]['KEPMAG'][1]
    data.pixel_images = pixel_images
    data.nearby = nearby
    data.hires = hires
    data.saturated = saturated
    data.bkg = bkg

    return data