astrobase/astrokep.py
def detrend_centroid(lcd, detrend='legendre', sigclip=None, mingap=0.5):
    '''Detrends the x and y coordinate centroids for a Kepler light curve.

    Given an `lcdict` for a single quarter of Kepler data, returned by
    `read_kepler_fitslc`, this function returns this same dictionary,
    appending detrended centroid_x and centroid_y values.

    Here "detrended" means "finite, SAP quality flag set to 0, sigma clipped,
    timegroups selected based on `mingap` day gaps, then fit vs time by a
    legendre polynomial of lowish degree".

    Parameters
    ----------

    lcd : lcdict
        An `lcdict` generated by the `read_kepler_fitslc` function.

    detrend : {'legendre'}
        Method by which to detrend the LC. 'legendre' is the only thing
        implemented at the moment.

    sigclip : None or float or int or sequence of floats/ints
        Determines the type and amount of sigma-clipping done on the light curve
        to remove outliers. If None, no sigma-clipping is performed. If a two
        element sequence of floats/ints, the first element corresponds to the
        fainter sigma-clip limit, and the second element corresponds to the
        brighter sigma-clip limit.

    mingap : float
        Number of days by which to define "timegroups" (for individual fitting
        each of timegroup, and to eliminate "burn-in" of Kepler spacecraft. For
        long cadence data, 0.5 days is typical.

    Returns
    -------

    tuple
        This is of the form `(lcd, errflag)`, where:

        `lcd` : an `lcdict` with the new key `lcd['centroids']`, containing the
        detrended times, (centroid_x, centroid_y) values, and their errors.

        `errflag` : boolean error flag, could be raised at various points.

    '''

    qnum = npunique(lcd['quarter'])
    try:
        assert qnum.size == 1, 'lcd should be for a unique quarter'
        assert detrend == 'legendre'
        qnum = int(qnum)
    except Exception as e:
        errflag = True

    # Get finite, QUALITY_FLAG != 0 times, centroids, and their errors.
    # Fraquelli & Thompson (2012), or perhaps also newer papers, give the list
    # of exclusions for quality flags.
    nbefore = lcd['time'].size

    # "ctd" for centroid.
    times = lcd['time'][lcd['sap_quality'] == 0]
    # Kepler Archive Manual KDMC-10008-006, pg 18. MOM_CENTR1 is the *column*
    # value for the flux-weighted centroid, MOM_CENTR2 is the row value.
    ctd_x = lcd['mom_centr2'][lcd['sap_quality'] == 0]
    ctd_y = lcd['mom_centr1'][lcd['sap_quality'] == 0]
    ctd_x_err = lcd['mom_centr2_err'][lcd['sap_quality'] == 0]
    ctd_y_err = lcd['mom_centr1_err'][lcd['sap_quality'] == 0]

    find = npisfinite(times) & npisfinite(ctd_x) & npisfinite(ctd_y)
    find &= (npisfinite(ctd_x_err)) & (npisfinite(ctd_y_err))

    f_times, f_ctd_x, f_ctd_y = times[find], ctd_x[find], ctd_y[find]
    f_ctd_x_err, f_ctd_y_err = ctd_x_err[find], ctd_y_err[find]

    # Sigma clip whopping outliers. It'd be better to have a general purpose
    # function for this, but sigclip_magseries works.
    stimes_x, s_ctd_x, s_ctd_x_err = sigclip_magseries(
        f_times,
        f_ctd_x,
        f_ctd_x_err,
        magsarefluxes=True,
        sigclip=30.0
    )
    stimes_y, s_ctd_y, s_ctd_y_err = sigclip_magseries(
        f_times,
        f_ctd_y,
        f_ctd_y_err,
        magsarefluxes=True,
        sigclip=30.0
    )

    # Get times and centroids where everything is finite and sigma clipped.
    mask_x = npin1d(stimes_x, stimes_y)
    s_times, s_ctd_x, s_ctd_x_err = (
        stimes_x[mask_x],
        s_ctd_x[mask_x],
        s_ctd_x_err[mask_x]
    )
    mask_y = npin1d(stimes_y, stimes_x)
    tmp, s_ctd_y, s_ctd_y_err = (
        stimes_y[mask_y],
        s_ctd_y[mask_y],
        s_ctd_y_err[mask_y]
    )
    try:
        np.testing.assert_array_equal(s_times, tmp)
        assert len(s_ctd_y) == len(s_times)
        assert len(s_ctd_y_err) == len(s_times)
        assert len(s_ctd_x) == len(s_times)
        assert len(s_ctd_x_err) == len(s_times)
    except AssertionError:
        return lcd, True

    nqflag = s_times.size

    # Drop intra-quarter and interquarter gaps in the timeseries. These are the
    # same limits set by Armstrong et al (2014): split each quarter's
    # timegroups by whether points are within 0.5 day limits. Then drop points
    # within 0.5 days of any boundary.  Finally, since the interquarter burn-in
    # time is more like 1 day, drop a further 0.5 days from the edges of each
    # quarter.  A nicer way to implement this would be with numpy masks, but
    # this approach just constructs the full arrays for any given quarter.

    ngroups, groups = find_lc_timegroups(s_times, mingap=mingap)
    tmp_times, tmp_ctd_x, tmp_ctd_y = [], [], []
    tmp_ctd_x_err, tmp_ctd_y_err = [], []

    for group in groups:
        tg_times = s_times[group]
        tg_ctd_x = s_ctd_x[group]
        tg_ctd_y = s_ctd_y[group]
        tg_ctd_x_err = s_ctd_x_err[group]
        tg_ctd_y_err = s_ctd_y_err[group]
        try:
            sel = ((tg_times > npmin(tg_times)+mingap) &
                   (tg_times < npmax(tg_times)-mingap))
        except ValueError:
            # If tgtimes is empty, continue to next timegroup.
            continue

        tmp_times.append(tg_times[sel])
        tmp_ctd_x.append(tg_ctd_x[sel])
        tmp_ctd_y.append(tg_ctd_y[sel])
        tmp_ctd_x_err.append(tg_ctd_x_err[sel])
        tmp_ctd_y_err.append(tg_ctd_y_err[sel])

    s_times,s_ctd_x,s_ctd_y,s_ctd_x_err,s_ctd_y_err = (
        nparray([]),nparray([]),nparray([]),nparray([]),nparray([])
    )

    # N.b.: works fine with empty arrays.
    for ix, _ in enumerate(tmp_times):
        s_times = npappend(s_times, tmp_times[ix])
        s_ctd_x = npappend(s_ctd_x, tmp_ctd_x[ix])
        s_ctd_y = npappend(s_ctd_y, tmp_ctd_y[ix])
        s_ctd_x_err = npappend(s_ctd_x_err, tmp_ctd_x_err[ix])
        s_ctd_y_err = npappend(s_ctd_y_err, tmp_ctd_y_err[ix])

    # Extra inter-quarter burn-in of 0.5 days.
    try:

        s_ctd_x = s_ctd_x[(s_times > (npmin(s_times)+mingap)) &
                          (s_times < (npmax(s_times)-mingap))]

    except Exception as e:
        # Case: s_times is wonky, all across this quarter. (Implemented because
        # of a rare bug with a singleton s_times array).
        LOGERROR('DETREND FAILED, qnum {:d}'.format(qnum))
        return npnan, True

    s_ctd_y = s_ctd_y[(s_times > (npmin(s_times)+mingap)) &
                      (s_times < (npmax(s_times)-mingap))]
    s_ctd_x_err = s_ctd_x_err[(s_times > (npmin(s_times)+mingap)) &
                              (s_times < (npmax(s_times)-mingap))]
    s_ctd_y_err = s_ctd_y_err[(s_times > (npmin(s_times)+mingap)) &
                              (s_times < (npmax(s_times)-mingap))]
    # Careful to do this last...
    s_times = s_times[(s_times > (npmin(s_times)+mingap)) &
                      (s_times < (npmax(s_times)-mingap))]

    nafter = s_times.size

    LOGINFO(
        'CLIPPING (SAP), qnum: {:d}'.format(qnum) +
        '\nndet before qflag & sigclip: {:d} ({:.3g}),'.format(
            nbefore, 1.
        ) +
        '\nndet after qflag & finite & sigclip: {:d} ({:.3g})'.format(
            nqflag, nqflag/float(nbefore)
        ) +
        '\nndet after dropping pts near gaps: {:d} ({:.3g})'.format(
            nafter, nafter/float(nbefore)
        )
    )

    # DETREND: fit a "low" order legendre series (see
    # "legendredeg_vs_npts_per_timegroup_ctd.pdf"), and save it to the output
    # dictionary. Save the fit (residuals to be computed after).
    ctd_dtr = {}

    if detrend == 'legendre':
        mingap = 0.5  # days
        ngroups, groups = find_lc_timegroups(s_times, mingap=mingap)
        tmpctdxlegfit, tmpctdylegfit, legdegs = [], [], []
        for group in groups:
            tg_times = s_times[group]
            tg_ctd_x = s_ctd_x[group]
            tg_ctd_x_err = s_ctd_x_err[group]
            tg_ctd_y = s_ctd_y[group]
            tg_ctd_y_err = s_ctd_y_err[group]

            legdeg = _get_legendre_deg_ctd(len(tg_times))
            tg_ctd_x_fit, _, _ = _legendre_dtr(tg_times,tg_ctd_x,tg_ctd_x_err,
                                               legendredeg=legdeg)
            tg_ctd_y_fit, _, _ = _legendre_dtr(tg_times,tg_ctd_y,tg_ctd_y_err,
                                               legendredeg=legdeg)

            tmpctdxlegfit.append(tg_ctd_x_fit)
            tmpctdylegfit.append(tg_ctd_y_fit)
            legdegs.append(legdeg)

        fit_ctd_x, fit_ctd_y = nparray([]), nparray([])
        for ix, _ in enumerate(tmpctdxlegfit):
            fit_ctd_x = npappend(fit_ctd_x, tmpctdxlegfit[ix])
            fit_ctd_y = npappend(fit_ctd_y, tmpctdylegfit[ix])

    ctd_dtr = {'times':s_times,
               'ctd_x':s_ctd_x,
               'ctd_x_err':s_ctd_x_err,
               'fit_ctd_x':fit_ctd_x,
               'ctd_y':s_ctd_y,
               'ctd_y_err':s_ctd_y_err,
               'fit_ctd_y':fit_ctd_y}

    lcd['ctd_dtr'] = ctd_dtr

    return lcd, False
def detrend_centroid(lcd, detrend='legendre', sigclip=None, mingap=0.5):
    '''Detrends the x and y coordinate centroids for a Kepler light curve.

    Given an `lcdict` for a single quarter of Kepler data, returned by
    `read_kepler_fitslc`, this function returns this same dictionary,
    appending detrended centroid_x and centroid_y values.

    Here "detrended" means "finite, SAP quality flag set to 0, sigma clipped,
    timegroups selected based on `mingap` day gaps, then fit vs time by a
    legendre polynomial of lowish degree".

    Parameters
    ----------

    lcd : lcdict
        An `lcdict` generated by the `read_kepler_fitslc` function.

    detrend : {'legendre'}
        Method by which to detrend the LC. 'legendre' is the only thing
        implemented at the moment.

    sigclip : None or float or int or sequence of floats/ints
        Determines the type and amount of sigma-clipping done on the light curve
        to remove outliers. If None, no sigma-clipping is performed. If a two
        element sequence of floats/ints, the first element corresponds to the
        fainter sigma-clip limit, and the second element corresponds to the
        brighter sigma-clip limit.

    mingap : float
        Number of days by which to define "timegroups" (for individual fitting
        each of timegroup, and to eliminate "burn-in" of Kepler spacecraft. For
        long cadence data, 0.5 days is typical.

    Returns
    -------

    tuple
        This is of the form `(lcd, errflag)`, where:

        `lcd` : an `lcdict` with the new key `lcd['centroids']`, containing the
        detrended times, (centroid_x, centroid_y) values, and their errors.

        `errflag` : boolean error flag, could be raised at various points.

    '''

    qnum = npunique(lcd['quarter'])
    try:
        assert qnum.size == 1, 'lcd should be for a unique quarter'
        assert detrend == 'legendre'
        qnum = int(qnum)
    except Exception as e:
        errflag = True

    # Get finite, QUALITY_FLAG != 0 times, centroids, and their errors.
    # Fraquelli & Thompson (2012), or perhaps also newer papers, give the list
    # of exclusions for quality flags.
    nbefore = lcd['time'].size

    # "ctd" for centroid.
    times = lcd['time'][lcd['sap_quality'] == 0]
    # Kepler Archive Manual KDMC-10008-006, pg 18. MOM_CENTR1 is the *column*
    # value for the flux-weighted centroid, MOM_CENTR2 is the row value.
    ctd_x = lcd['mom_centr2'][lcd['sap_quality'] == 0]
    ctd_y = lcd['mom_centr1'][lcd['sap_quality'] == 0]
    ctd_x_err = lcd['mom_centr2_err'][lcd['sap_quality'] == 0]
    ctd_y_err = lcd['mom_centr1_err'][lcd['sap_quality'] == 0]

    find = npisfinite(times) & npisfinite(ctd_x) & npisfinite(ctd_y)
    find &= (npisfinite(ctd_x_err)) & (npisfinite(ctd_y_err))

    f_times, f_ctd_x, f_ctd_y = times[find], ctd_x[find], ctd_y[find]
    f_ctd_x_err, f_ctd_y_err = ctd_x_err[find], ctd_y_err[find]

    # Sigma clip whopping outliers. It'd be better to have a general purpose
    # function for this, but sigclip_magseries works.
    stimes_x, s_ctd_x, s_ctd_x_err = sigclip_magseries(
        f_times,
        f_ctd_x,
        f_ctd_x_err,
        magsarefluxes=True,
        sigclip=30.0
    )
    stimes_y, s_ctd_y, s_ctd_y_err = sigclip_magseries(
        f_times,
        f_ctd_y,
        f_ctd_y_err,
        magsarefluxes=True,
        sigclip=30.0
    )

    # Get times and centroids where everything is finite and sigma clipped.
    mask_x = npin1d(stimes_x, stimes_y)
    s_times, s_ctd_x, s_ctd_x_err = (
        stimes_x[mask_x],
        s_ctd_x[mask_x],
        s_ctd_x_err[mask_x]
    )
    mask_y = npin1d(stimes_y, stimes_x)
    tmp, s_ctd_y, s_ctd_y_err = (
        stimes_y[mask_y],
        s_ctd_y[mask_y],
        s_ctd_y_err[mask_y]
    )
    try:
        np.testing.assert_array_equal(s_times, tmp)
        assert len(s_ctd_y) == len(s_times)
        assert len(s_ctd_y_err) == len(s_times)
        assert len(s_ctd_x) == len(s_times)
        assert len(s_ctd_x_err) == len(s_times)
    except AssertionError:
        return lcd, True

    nqflag = s_times.size

    # Drop intra-quarter and interquarter gaps in the timeseries. These are the
    # same limits set by Armstrong et al (2014): split each quarter's
    # timegroups by whether points are within 0.5 day limits. Then drop points
    # within 0.5 days of any boundary.  Finally, since the interquarter burn-in
    # time is more like 1 day, drop a further 0.5 days from the edges of each
    # quarter.  A nicer way to implement this would be with numpy masks, but
    # this approach just constructs the full arrays for any given quarter.

    ngroups, groups = find_lc_timegroups(s_times, mingap=mingap)
    tmp_times, tmp_ctd_x, tmp_ctd_y = [], [], []
    tmp_ctd_x_err, tmp_ctd_y_err = [], []

    for group in groups:
        tg_times = s_times[group]
        tg_ctd_x = s_ctd_x[group]
        tg_ctd_y = s_ctd_y[group]
        tg_ctd_x_err = s_ctd_x_err[group]
        tg_ctd_y_err = s_ctd_y_err[group]
        try:
            sel = ((tg_times > npmin(tg_times)+mingap) &
                   (tg_times < npmax(tg_times)-mingap))
        except ValueError:
            # If tgtimes is empty, continue to next timegroup.
            continue

        tmp_times.append(tg_times[sel])
        tmp_ctd_x.append(tg_ctd_x[sel])
        tmp_ctd_y.append(tg_ctd_y[sel])
        tmp_ctd_x_err.append(tg_ctd_x_err[sel])
        tmp_ctd_y_err.append(tg_ctd_y_err[sel])

    s_times,s_ctd_x,s_ctd_y,s_ctd_x_err,s_ctd_y_err = (
        nparray([]),nparray([]),nparray([]),nparray([]),nparray([])
    )

    # N.b.: works fine with empty arrays.
    for ix, _ in enumerate(tmp_times):
        s_times = npappend(s_times, tmp_times[ix])
        s_ctd_x = npappend(s_ctd_x, tmp_ctd_x[ix])
        s_ctd_y = npappend(s_ctd_y, tmp_ctd_y[ix])
        s_ctd_x_err = npappend(s_ctd_x_err, tmp_ctd_x_err[ix])
        s_ctd_y_err = npappend(s_ctd_y_err, tmp_ctd_y_err[ix])

    # Extra inter-quarter burn-in of 0.5 days.
    try:

        s_ctd_x = s_ctd_x[(s_times > (npmin(s_times)+mingap)) &
                          (s_times < (npmax(s_times)-mingap))]

    except Exception as e:
        # Case: s_times is wonky, all across this quarter. (Implemented because
        # of a rare bug with a singleton s_times array).
        LOGERROR('DETREND FAILED, qnum {:d}'.format(qnum))
        return npnan, True

    s_ctd_y = s_ctd_y[(s_times > (npmin(s_times)+mingap)) &
                      (s_times < (npmax(s_times)-mingap))]
    s_ctd_x_err = s_ctd_x_err[(s_times > (npmin(s_times)+mingap)) &
                              (s_times < (npmax(s_times)-mingap))]
    s_ctd_y_err = s_ctd_y_err[(s_times > (npmin(s_times)+mingap)) &
                              (s_times < (npmax(s_times)-mingap))]
    # Careful to do this last...
    s_times = s_times[(s_times > (npmin(s_times)+mingap)) &
                      (s_times < (npmax(s_times)-mingap))]

    nafter = s_times.size

    LOGINFO(
        'CLIPPING (SAP), qnum: {:d}'.format(qnum) +
        '\nndet before qflag & sigclip: {:d} ({:.3g}),'.format(
            nbefore, 1.
        ) +
        '\nndet after qflag & finite & sigclip: {:d} ({:.3g})'.format(
            nqflag, nqflag/float(nbefore)
        ) +
        '\nndet after dropping pts near gaps: {:d} ({:.3g})'.format(
            nafter, nafter/float(nbefore)
        )
    )

    # DETREND: fit a "low" order legendre series (see
    # "legendredeg_vs_npts_per_timegroup_ctd.pdf"), and save it to the output
    # dictionary. Save the fit (residuals to be computed after).
    ctd_dtr = {}

    if detrend == 'legendre':
        mingap = 0.5  # days
        ngroups, groups = find_lc_timegroups(s_times, mingap=mingap)
        tmpctdxlegfit, tmpctdylegfit, legdegs = [], [], []
        for group in groups:
            tg_times = s_times[group]
            tg_ctd_x = s_ctd_x[group]
            tg_ctd_x_err = s_ctd_x_err[group]
            tg_ctd_y = s_ctd_y[group]
            tg_ctd_y_err = s_ctd_y_err[group]

            legdeg = _get_legendre_deg_ctd(len(tg_times))
            tg_ctd_x_fit, _, _ = _legendre_dtr(tg_times,tg_ctd_x,tg_ctd_x_err,
                                               legendredeg=legdeg)
            tg_ctd_y_fit, _, _ = _legendre_dtr(tg_times,tg_ctd_y,tg_ctd_y_err,
                                               legendredeg=legdeg)

            tmpctdxlegfit.append(tg_ctd_x_fit)
            tmpctdylegfit.append(tg_ctd_y_fit)
            legdegs.append(legdeg)

        fit_ctd_x, fit_ctd_y = nparray([]), nparray([])
        for ix, _ in enumerate(tmpctdxlegfit):
            fit_ctd_x = npappend(fit_ctd_x, tmpctdxlegfit[ix])
            fit_ctd_y = npappend(fit_ctd_y, tmpctdylegfit[ix])

    ctd_dtr = {'times':s_times,
               'ctd_x':s_ctd_x,
               'ctd_x_err':s_ctd_x_err,
               'fit_ctd_x':fit_ctd_x,
               'ctd_y':s_ctd_y,
               'ctd_y_err':s_ctd_y_err,
               'fit_ctd_y':fit_ctd_y}

    lcd['ctd_dtr'] = ctd_dtr

    return lcd, False
def detrend_centroid ( lcd , detrend = 'legendre' , sigclip = None , mingap = 0.5 ) : qnum = npunique ( lcd [ 'quarter' ] ) try : assert qnum . size == 1 , 'lcd should be for a unique quarter' assert detrend == 'legendre' qnum = int ( qnum ) except Exception as e : errflag = True # Get finite, QUALITY_FLAG != 0 times, centroids, and their errors. # Fraquelli & Thompson (2012), or perhaps also newer papers, give the list # of exclusions for quality flags. nbefore = lcd [ 'time' ] . size # "ctd" for centroid. times = lcd [ 'time' ] [ lcd [ 'sap_quality' ] == 0 ] # Kepler Archive Manual KDMC-10008-006, pg 18. MOM_CENTR1 is the *column* # value for the flux-weighted centroid, MOM_CENTR2 is the row value. ctd_x = lcd [ 'mom_centr2' ] [ lcd [ 'sap_quality' ] == 0 ] ctd_y = lcd [ 'mom_centr1' ] [ lcd [ 'sap_quality' ] == 0 ] ctd_x_err = lcd [ 'mom_centr2_err' ] [ lcd [ 'sap_quality' ] == 0 ] ctd_y_err = lcd [ 'mom_centr1_err' ] [ lcd [ 'sap_quality' ] == 0 ] find = npisfinite ( times ) & npisfinite ( ctd_x ) & npisfinite ( ctd_y ) find &= ( npisfinite ( ctd_x_err ) ) & ( npisfinite ( ctd_y_err ) ) f_times , f_ctd_x , f_ctd_y = times [ find ] , ctd_x [ find ] , ctd_y [ find ] f_ctd_x_err , f_ctd_y_err = ctd_x_err [ find ] , ctd_y_err [ find ] # Sigma clip whopping outliers. It'd be better to have a general purpose # function for this, but sigclip_magseries works. stimes_x , s_ctd_x , s_ctd_x_err = sigclip_magseries ( f_times , f_ctd_x , f_ctd_x_err , magsarefluxes = True , sigclip = 30.0 ) stimes_y , s_ctd_y , s_ctd_y_err = sigclip_magseries ( f_times , f_ctd_y , f_ctd_y_err , magsarefluxes = True , sigclip = 30.0 ) # Get times and centroids where everything is finite and sigma clipped. mask_x = npin1d ( stimes_x , stimes_y ) s_times , s_ctd_x , s_ctd_x_err = ( stimes_x [ mask_x ] , s_ctd_x [ mask_x ] , s_ctd_x_err [ mask_x ] ) mask_y = npin1d ( stimes_y , stimes_x ) tmp , s_ctd_y , s_ctd_y_err = ( stimes_y [ mask_y ] , s_ctd_y [ mask_y ] , s_ctd_y_err [ mask_y ] ) try : np . testing . assert_array_equal ( s_times , tmp ) assert len ( s_ctd_y ) == len ( s_times ) assert len ( s_ctd_y_err ) == len ( s_times ) assert len ( s_ctd_x ) == len ( s_times ) assert len ( s_ctd_x_err ) == len ( s_times ) except AssertionError : return lcd , True nqflag = s_times . size # Drop intra-quarter and interquarter gaps in the timeseries. These are the # same limits set by Armstrong et al (2014): split each quarter's # timegroups by whether points are within 0.5 day limits. Then drop points # within 0.5 days of any boundary.  Finally, since the interquarter burn-in # time is more like 1 day, drop a further 0.5 days from the edges of each # quarter.  A nicer way to implement this would be with numpy masks, but # this approach just constructs the full arrays for any given quarter. ngroups , groups = find_lc_timegroups ( s_times , mingap = mingap ) tmp_times , tmp_ctd_x , tmp_ctd_y = [ ] , [ ] , [ ] tmp_ctd_x_err , tmp_ctd_y_err = [ ] , [ ] for group in groups : tg_times = s_times [ group ] tg_ctd_x = s_ctd_x [ group ] tg_ctd_y = s_ctd_y [ group ] tg_ctd_x_err = s_ctd_x_err [ group ] tg_ctd_y_err = s_ctd_y_err [ group ] try : sel = ( ( tg_times > npmin ( tg_times ) + mingap ) & ( tg_times < npmax ( tg_times ) - mingap ) ) except ValueError : # If tgtimes is empty, continue to next timegroup. continue tmp_times . append ( tg_times [ sel ] ) tmp_ctd_x . append ( tg_ctd_x [ sel ] ) tmp_ctd_y . append ( tg_ctd_y [ sel ] ) tmp_ctd_x_err . append ( tg_ctd_x_err [ sel ] ) tmp_ctd_y_err . append ( tg_ctd_y_err [ sel ] ) s_times , s_ctd_x , s_ctd_y , s_ctd_x_err , s_ctd_y_err = ( nparray ( [ ] ) , nparray ( [ ] ) , nparray ( [ ] ) , nparray ( [ ] ) , nparray ( [ ] ) ) # N.b.: works fine with empty arrays. for ix , _ in enumerate ( tmp_times ) : s_times = npappend ( s_times , tmp_times [ ix ] ) s_ctd_x = npappend ( s_ctd_x , tmp_ctd_x [ ix ] ) s_ctd_y = npappend ( s_ctd_y , tmp_ctd_y [ ix ] ) s_ctd_x_err = npappend ( s_ctd_x_err , tmp_ctd_x_err [ ix ] ) s_ctd_y_err = npappend ( s_ctd_y_err , tmp_ctd_y_err [ ix ] ) # Extra inter-quarter burn-in of 0.5 days. try : s_ctd_x = s_ctd_x [ ( s_times > ( npmin ( s_times ) + mingap ) ) & ( s_times < ( npmax ( s_times ) - mingap ) ) ] except Exception as e : # Case: s_times is wonky, all across this quarter. (Implemented because # of a rare bug with a singleton s_times array). LOGERROR ( 'DETREND FAILED, qnum {:d}' . format ( qnum ) ) return npnan , True s_ctd_y = s_ctd_y [ ( s_times > ( npmin ( s_times ) + mingap ) ) & ( s_times < ( npmax ( s_times ) - mingap ) ) ] s_ctd_x_err = s_ctd_x_err [ ( s_times > ( npmin ( s_times ) + mingap ) ) & ( s_times < ( npmax ( s_times ) - mingap ) ) ] s_ctd_y_err = s_ctd_y_err [ ( s_times > ( npmin ( s_times ) + mingap ) ) & ( s_times < ( npmax ( s_times ) - mingap ) ) ] # Careful to do this last... s_times = s_times [ ( s_times > ( npmin ( s_times ) + mingap ) ) & ( s_times < ( npmax ( s_times ) - mingap ) ) ] nafter = s_times . size LOGINFO ( 'CLIPPING (SAP), qnum: {:d}' . format ( qnum ) + '\nndet before qflag & sigclip: {:d} ({:.3g}),' . format ( nbefore , 1. ) + '\nndet after qflag & finite & sigclip: {:d} ({:.3g})' . format ( nqflag , nqflag / float ( nbefore ) ) + '\nndet after dropping pts near gaps: {:d} ({:.3g})' . format ( nafter , nafter / float ( nbefore ) ) ) # DETREND: fit a "low" order legendre series (see # "legendredeg_vs_npts_per_timegroup_ctd.pdf"), and save it to the output # dictionary. Save the fit (residuals to be computed after). ctd_dtr = { } if detrend == 'legendre' : mingap = 0.5 # days ngroups , groups = find_lc_timegroups ( s_times , mingap = mingap ) tmpctdxlegfit , tmpctdylegfit , legdegs = [ ] , [ ] , [ ] for group in groups : tg_times = s_times [ group ] tg_ctd_x = s_ctd_x [ group ] tg_ctd_x_err = s_ctd_x_err [ group ] tg_ctd_y = s_ctd_y [ group ] tg_ctd_y_err = s_ctd_y_err [ group ] legdeg = _get_legendre_deg_ctd ( len ( tg_times ) ) tg_ctd_x_fit , _ , _ = _legendre_dtr ( tg_times , tg_ctd_x , tg_ctd_x_err , legendredeg = legdeg ) tg_ctd_y_fit , _ , _ = _legendre_dtr ( tg_times , tg_ctd_y , tg_ctd_y_err , legendredeg = legdeg ) tmpctdxlegfit . append ( tg_ctd_x_fit ) tmpctdylegfit . append ( tg_ctd_y_fit ) legdegs . append ( legdeg ) fit_ctd_x , fit_ctd_y = nparray ( [ ] ) , nparray ( [ ] ) for ix , _ in enumerate ( tmpctdxlegfit ) : fit_ctd_x = npappend ( fit_ctd_x , tmpctdxlegfit [ ix ] ) fit_ctd_y = npappend ( fit_ctd_y , tmpctdylegfit [ ix ] ) ctd_dtr = { 'times' : s_times , 'ctd_x' : s_ctd_x , 'ctd_x_err' : s_ctd_x_err , 'fit_ctd_x' : fit_ctd_x , 'ctd_y' : s_ctd_y , 'ctd_y_err' : s_ctd_y_err , 'fit_ctd_y' : fit_ctd_y } lcd [ 'ctd_dtr' ] = ctd_dtr return lcd , False
Detrends the x and y coordinate centroids for a Kepler light curve.

    Given an `lcdict` for a single quarter of Kepler data, returned by
    `read_kepler_fitslc`, this function returns this same dictionary,
    appending detrended centroid_x and centroid_y values.

    Here "detrended" means "finite, SAP quality flag set to 0, sigma clipped,
    timegroups selected based on `mingap` day gaps, then fit vs time by a
    legendre polynomial of lowish degree".

    Parameters
    ----------

    lcd : lcdict
        An `lcdict` generated by the `read_kepler_fitslc` function.

    detrend : {'legendre'}
        Method by which to detrend the LC. 'legendre' is the only thing
        implemented at the moment.

    sigclip : None or float or int or sequence of floats/ints
        Determines the type and amount of sigma-clipping done on the light curve
        to remove outliers. If None, no sigma-clipping is performed. If a two
        element sequence of floats/ints, the first element corresponds to the
        fainter sigma-clip limit, and the second element corresponds to the
        brighter sigma-clip limit.

    mingap : float
        Number of days by which to define "timegroups" (for individual fitting
        each of timegroup, and to eliminate "burn-in" of Kepler spacecraft. For
        long cadence data, 0.5 days is typical.

    Returns
    -------

    tuple
        This is of the form `(lcd, errflag)`, where:

        `lcd` : an `lcdict` with the new key `lcd['centroids']`, containing the
        detrended times, (centroid_x, centroid_y) values, and their errors.

        `errflag` : boolean error flag, could be raised at various points.
Detrends the x and y coordinate centroids for a Kepler light curve .