Python HelperFunctions.astropy_smooth方法代码示例

# 需要导入模块: import HelperFunctions [as 别名]
# 或者: from HelperFunctions import astropy_smooth [as 别名]
def process_model(model, data, vsini_model=None, resolution=None, vsini_primary=None,
                  maxvel=1000.0, debug=False, oversample=1, logspace=True):
    # Read in the model if necessary
    if isinstance(model, str):
        if debug:
            print "Reading in the input model from %s" % model
        x, y = np.loadtxt(model, usecols=(0, 1), unpack=True)
        x = x * u.angstrom.to(u.nm)
        y = 10 ** y
        left = np.searchsorted(x, data[0].x[0] - 10)
        right = np.searchsorted(x, data[-1].x[-1] + 10)
        model = DataStructures.xypoint(x=x[left:right], y=y[left:right])
    elif not isinstance(model, DataStructures.xypoint):
        raise TypeError(
            "Input model is of an unknown type! Must be a DataStructures.xypoint or a string with the filename.")


    # Linearize the x-axis of the model (in log-spacing)
    if logspace:
        if debug:
            print "Linearizing model"
        xgrid = np.logspace(np.log10(model.x[0]), np.log10(model.x[-1]), model.size())
        model = FittingUtilities.RebinData(model, xgrid)
    else:
        xgrid = model.x

    # Broaden
    if vsini_model is not None and vsini_model > 1.0 * u.km.to(u.cm):
        if debug:
            print "Rotationally broadening model to vsini = %g km/s" % (vsini_model * u.cm.to(u.km))
        model = Broaden.RotBroad(model, vsini_model, linear=True)


    # Reduce resolution
    if resolution is not None and 5000 < resolution < 500000:
        if debug:
            print "Convolving to the detector resolution of %g" % resolution
        model = FittingUtilities.ReduceResolutionFFT(model, resolution)

    # Divide by the same smoothing kernel as we used for the data
    if vsini_primary is not None:
        """ old method
        smooth_factor = 0.5
        d_logx = np.log(xgrid[1] / xgrid[0])
        theta = GenericSmooth.roundodd(vsini_primary / 3e5 * smooth_factor / d_logx)
        print "Window size = {}\ndlogx = {}\nvsini = {}".format(theta, d_logx, vsini_primary)
        smooth = FittingUtilities.savitzky_golay(model.y, theta, 5)

        model.y = model.y - smooth
        """
        # model.y = HelperFunctions.HighPassFilter(model, vel=smooth_factor * vsini_primary * u.km.to(u.cm),
        #                                         linearize=True)
        #model.y -= HelperFunctions.astropy_smooth(model, vel=SMOOTH_FACTOR * vsini_primary, linearize=True)
        smoothed = HelperFunctions.astropy_smooth(model, vel=SMOOTH_FACTOR * vsini_primary, linearize=False)
        #minval = min(model.y)
        #model.y += abs(minval)
        #model.cont = abs(minval) * np.ones(model.size())
        model.y += model.cont.mean() - smoothed
        model.cont = np.ones(model.size()) * model.cont.mean()




    # Rebin subsets of the model to the same spacing as the data
    model_orders = []
    model_fcn = spline(model.x, model.y)
    if debug:
        model.output("Test_model.dat")
    for i, order in enumerate(data):
        if debug:
            sys.stdout.write("\rGenerating model subset for order %i in the input data" % (i + 1))
            sys.stdout.flush()
        # Find how much to extend the model so that we can get maxvel range.
        dlambda = order.x[order.size() / 2] * maxvel * 1.5 / 3e5
        left = np.searchsorted(model.x, order.x[0] - dlambda)
        right = np.searchsorted(model.x, order.x[-1] + dlambda)
        right = min(right, model.size() - 2)

        # Figure out the log-spacing of the data
        logspacing = np.log(order.x[1] / order.x[0])

        # Finally, space the model segment with the same log-spacing
        start = np.log(model.x[left])
        end = np.log(model.x[right])
        xgrid = np.exp(np.arange(start, end + logspacing, logspacing))

        segment = DataStructures.xypoint(x=xgrid, y=model_fcn(xgrid))
        # segment = FittingUtilities.RebinData(model[left:right + 1].copy(), xgrid)
        segment.cont = FittingUtilities.Continuum(segment.x, segment.y, lowreject=1.5, highreject=5, fitorder=2)
        model_orders.append(segment)

    print "\n"
    return model_orders

# 需要导入模块: import HelperFunctions [as 别名]
# 或者: from HelperFunctions import astropy_smooth [as 别名]
def Process_Data_serial(fname, badregions=[], interp_regions=[], extensions=True,
                 trimsize=1, vsini=None, logspacing=False, oversample=1.0, reject_outliers=True):
    """

    :param fname: The filename to read in (should be a fits file)
    :param badregions: a list of regions to exclude (contains strong telluric or stellar line residuals)
    :param interp_regions: a list of regions to interpolate over
    :param extensions: A boolean flag for whether the fits file is separated into extensions
    :param trimsize: The amount to exclude from both ends of every order (where it is very noisy)
    :param vsini: the primary star vsini. If given subtract an estimate of the primary star model obtained by
                  denoising and smoothing with a kernel size set by the vsini.
    :logspacing: If true, interpolate each order into a constant log-spacing.
    :return:
    """
    if isinstance(fname, list) and all([isinstance(f, DataStructures.xypoint) for f in fname]):
        orders = fname
    else:
        if extensions:
            orders = HelperFunctions.ReadExtensionFits(fname)

        else:
            orders = HelperFunctions.ReadFits(fname, errors=2)

    numorders = len(orders)
    for i, order in enumerate(orders[::-1]):
        # Trim data, and make sure the wavelength spacing is constant
        if trimsize > 0:
            order = order[trimsize:-trimsize]

        # Smooth the data
        if vsini is not None:
            # make sure the x-spacing is linear
            xgrid = np.linspace(order.x[0], order.x[-1], order.size())
            order = FittingUtilities.RebinData(order, xgrid)

            dx = order.x[1] - order.x[0]
            """ old method
            theta = GenericSmooth.roundodd(vsini / 3e5 * order.x.mean() / dx * smooth_factor)
            theta = max(theta, 21)
            #denoised = HelperFunctions.Denoise(order.copy())
            denoised = order.copy()
            smooth = FittingUtilities.savitzky_golay(denoised.y, theta, 5)
            order.y = order.y - smooth + order.cont.mean()
            """
            # order.y = order.cont.mean() + HelperFunctions.HighPassFilter(order,
            #                                                             vel=smooth_factor * vsini * u.km.to(u.cm))
            #smoothed = HelperFunctions.astropy_smooth(order, vel=SMOOTH_FACTOR * vsini, linearize=False)
            smoothed = HelperFunctions.astropy_smooth(order, vel=SMOOTH_FACTOR * vsini, linearize=True)
            order.y += order.cont.mean() - smoothed
            order.cont = np.ones(order.size()) * order.cont.mean()

        # Remove bad regions from the data
        for region in badregions:
            left = np.searchsorted(order.x, region[0])
            right = np.searchsorted(order.x, region[1])
            if left > 0 and right < order.size():
                print "Warning! Bad region covers the middle of order %i" % i
                print "Removing full order!"
                left = 0
                right = order.size()
            order.x = np.delete(order.x, np.arange(left, right))
            order.y = np.delete(order.y, np.arange(left, right))
            order.cont = np.delete(order.cont, np.arange(left, right))
            order.err = np.delete(order.err, np.arange(left, right))

        # Interpolate over interp_regions:
        for region in interp_regions:
            left = np.searchsorted(order.x, region[0])
            right = np.searchsorted(order.x, region[1])
            order.y[left:right] = order.cont[left:right]


        # Remove whole order if it is too small
        remove = False
        if order.x.size <= 1:
            remove = True
        else:
            velrange = 3e5 * (np.median(order.x) - order.x[0]) / np.median(order.x)
            if velrange <= 1050.0:
                remove = True
        if remove:
            print "Removing order %i" % (numorders - 1 - i)
            orders.pop(numorders - 1 - i)
        else:
            if reject_outliers:
                # Find outliers from e.g. bad telluric line or stellar spectrum removal.
                order.cont = FittingUtilities.Continuum(order.x, order.y, lowreject=3, highreject=3)
                outliers = HelperFunctions.FindOutliers(order, expand=10, numsiglow=5, numsighigh=5)
                # plt.plot(order.x, order.y / order.cont, 'k-')
                if len(outliers) > 0:
                    # plt.plot(order.x[outliers], (order.y / order.cont)[outliers], 'r-')
                    order.y[outliers] = order.cont[outliers]
                    order.cont = FittingUtilities.Continuum(order.x, order.y, lowreject=3, highreject=3)
                    order.y[outliers] = order.cont[outliers]

            # Save this order
            orders[numorders - 1 - i] = order.copy()

    # plt.show()

#.........这里部分代码省略.........