Python sep.extract函数代码示例

def find_sources_with_sep(img):
    """Return sources (x, y) sorted by brightness. Use SEP package.
    """
    import sep
    if isinstance(img, np.ma.MaskedArray):
        image = img.filled(fill_value=np.median(img)).astype('float32')
    else:
        image = img.astype('float32')

    bkg = sep.Background(image)
    thresh = 3. * bkg.globalrms
    try:
        sources = sep.extract(image - bkg.back(), thresh)
    except Exception as e:
        buff_message = 'internal pixel buffer full'
        if e.message[0:26] == buff_message:
            sep.set_extract_pixstack(600000)
        try:
            sources = sep.extract(image - bkg.back(), thresh)
        except Exception as e:
            if e.message[0:26] == buff_message:
                sep.set_extract_pixstack(900000)
                sources = sep.extract(image - bkg.back(), thresh)

    sources.sort(order='flux')
    return np.array([[asrc['x'], asrc['y']] for asrc in sources[::-1]])

        def compare_image(the_image):
            """Return the fraction of sources found in the original image"""
            # pixel comparison is not good, doesn't work. Compare catalogs.
            if isinstance(the_image, np.ma.MaskedArray):
                full_algn = the_image.filled(fill_value=np.median(the_image))\
                    .astype('float32')
            else:
                full_algn = the_image.astype('float32')
            full_algn[the_image == 0] = np.median(the_image)
            import sep
            bkg = sep.Background(full_algn)
            thresh = 3.0 * bkg.globalrms
            allobjs = sep.extract(full_algn - bkg.back(), thresh)
            allxy = np.array([[obj['x'], obj['y']] for obj in allobjs])

            from scipy.spatial import KDTree
            ref_coordtree = KDTree(self.star_new_pos)

            # Compare here srcs list with self.star_ref_pos
            num_sources = 0
            for asrc in allxy:
                found_source = ref_coordtree.query_ball_point(asrc, 3)
                if found_source:
                    num_sources += 1
            fraction_found = float(num_sources) / float(len(allxy))
            return fraction_found

	def __sepFindFWHM(self,tries):
		from astropy.io import fits
		import math
		import traceback
		focpos=[]
		fwhm=[]
		fwhm_min=None
		fwhm_MinimumX=None
		keys = list(tries.keys())
		keys.sort()
		ln2=math.log(2)
		for k in keys:
			try:
				fwhms=[]
				ff=fits.open(tries[k])
				# loop on images..
				for i in range(1,len(ff)):
					data=ff[i].data
					bkg=sep.Background(numpy.array(data,numpy.float))
					sources=sep.extract(data-bkg, 5.0 * bkg.globalrms)
					for s in sources:
						fwhms.append(2 * math.sqrt(ln2 * (s[12]**2 + s[13]**2)))
				im_fwhm=numpy.median(fwhms)
				# find median from fwhms measurements..
				self.log('I','offset {0} fwhm {1} with {2} stars'.format(k,im_fwhm,len(fwhms)))
				focpos.append(k)
				fwhm.append(im_fwhm)
				if (fwhm_min is None or im_fwhm < fwhm_min):
					fwhm_MinimumX = k
					fwhm_min = im_fwhm
			except Exception as ex:
				self.log('W','offset {0}: {1} {2}'.format(k,ex,traceback.format_exc()))
		return focpos,fwhm,fwhm_min,fwhm_MinimumX

def sep_phot(data, ap, th):
    """
    Preforms photometry by SEP, similar to source extractor
    """

    # Measure a spatially variable background of some image data (np array)
    try:
        bkg = sep.Background(data)  # , mask=mask, bw=64, bh=64, fw=3, fh=3) # optional parameters
    except ValueError:
        data = data.byteswap(True).newbyteorder()
        bkg = sep.Background(data)  # , mask=mask, bw=64, bh=64, fw=3, fh=3) # optional parameters

    # Directly subtract the background from the data in place
    bkg.subfrom(data)

    # for the background subtracted data, detect objects in data given some threshold
    thresh = th * bkg.globalrms  # ensure the threshold is high enough wrt background
    objs = sep.extract(data, thresh)
    # calculate the Kron radius for each object, then we perform elliptical aperture photometry within that radius
    kronrad, krflag = sep.kron_radius(data, objs['x'], objs['y'], objs['a'], objs['b'], objs['theta'], ap)
    flux, fluxerr, flag = sep.sum_ellipse(data, objs['x'], objs['y'], objs['a'], objs['b'], objs['theta'],
                                          2.5 * kronrad, subpix=1)
    flag |= krflag  # combine flags into 'flag'

    r_min = 1.75  # minimum diameter = 3.5
    use_circle = kronrad * np.sqrt(objs['a'] * objs['b']) < r_min
    x = objs['x']
    y = objs['y']
    cflux, cfluxerr, cflag = sep.sum_circle(data, x[use_circle], y[use_circle],
                                            r_min, subpix=1)
    flux[use_circle] = cflux
    fluxerr[use_circle] = cfluxerr
    flag[use_circle] = cflag

    return objs

def detect_sources(data, threshold, elipse = False, *args, **kwargs):
    '''
    Detect the sources of the image.

    Parameters:
        data : ~numpy.ndarray~
            2D image data.
        threshold : float or ~numpy.ndarray~
            The threshold of the detection.
        elipse : bool
            Tell the program if you want the elipses parameters for each object.
        **kwargs will be passed integrally to the sep functions.

    Returns:
        x, y : ~numpy.ndarray~
            The positions of the detected sources.
        a, b, theta : ~numpy.ndarray~
            The parameters of the detected elipses.
    '''
    data = _fix_data(data)
    objs = sep.extract(data, threshold, **kwargs)

    if elipse:
        return objs['x'], objs['y'], objs['a'], objs['b'], objs['theta']
    else:
        return objs['x'], objs['y']

def extract_sources(data):
	data = data.byteswap(True).newbyteorder()
	bkg = sep.Background(data)
	back = bkg.back()
	data_woback = data-back
	thresh = 1.5 * bkg.globalrms
	objects = sep.extract(data_woback, thresh)
	return objects

    def run(self, step_name, mask=None):

        try:
            img = self.image.copy()
        except AttributeError:
            logger.error('You must setup image before running pipeline!')
            exit(1)

        logger.info('running ' + step_name)

        kws = self.step_kws[step_name].copy()

        sep_extract_kws = kws.pop('sep_extract_kws', {})
        sep_extract_kws = check_kwargs_defaults(sep_extract_kws,
                                                self.SEP_EXTRACT_DEFAULTS)

        sep_back_kws = kws.pop('sep_back_kws', {})
        sep_back_kws = check_kwargs_defaults(sep_back_kws,
                                             self.SEP_BACK_DEFAULTS)

        if mask is not None:
            logger.info('applying mask') 
            sep_extract_kws['mask'] = mask
            sep_back_kws['mask'] = mask

        # define gaussian kernel for detection
        # TODO: make filter shape optional
        num_fwhm = sep_extract_kws.pop('filter_num_fwhm', 1)
        kernel_fwhm = self.psf_fwhm * num_fwhm
        logger.info('smoothing with kernel with fwhm = {:.2f} arcsec'.\
                    format(kernel_fwhm * utils.pixscale))
        kern = Gaussian2DKernel(kernel_fwhm * gaussian_fwhm_to_sigma, 
                                mode='oversample')
        kern.normalize()
        sep_extract_kws['filter_kernel'] = kern.array

        # estimate background and subtract from image
        bkg = sep.Background(img, **sep_back_kws)
        img_sub = img - bkg

        # extract sources
        logger.info('detecting with a threshold of {} x background'.\
                    format(sep_extract_kws['thresh']))
        sources, segmap =  sep.extract(
            img_sub, err=bkg.rms(), segmentation_map=True, **sep_extract_kws)

        sources = Table(sources)
        sources = sources[sources['flux'] > 0]

        logger.info('found {} sources'.format(len(sources)))

        if kws['do_measure']:
            sources = self._measure(img, sources, mask)

        sources['seg_label'] = np.arange(1, len(sources) + 1)
        self.sources[step_name] = sources

        return sources, segmap

def createMask(data, thresh=100, title="Mask", plotMask=True):
    extraction = sep.extract(data, thresh, segmentation_map=True)
    mask = extraction[1]
    if plotMask:
        plt.figure()
        plt.imshow(mask, aspect="auto", interpolation="nearest", origin="lower")
        plt.title(title)
        plt.show()
    return mask

def test_long_error_msg():
    """Ensure that the error message is created successfully when
    there is an error detail."""

    # set extract pixstack to an insanely small value; this will trigger
    # a detailed error message when running sep.extract()
    old = sep.get_extract_pixstack()
    sep.set_extract_pixstack(5)

    data = np.ones((10, 10), dtype=np.float64)
    with pytest.raises(Exception) as excinfo:
        sep.extract(data, 0.1)
    msg = excinfo.value.args[0]
    assert type(msg) == str  # check that message is the native string type
    assert msg.startswith("internal pixel buffer full: The limit")

    # restore
    sep.set_extract_pixstack(old)

        def auto_fit(image_data):
            an_width = 3
            temp_image_data = image_data
            internal_image_data = temp_image_data.byteswap(True).newbyteorder()
            bkg = sep.Background(internal_image_data)
            thresh = 1.5 * bkg.globalback
            objects = sep.extract(internal_image_data, thresh)
            center_x = internal_image_data.shape[0]/2.0
            center_y = internal_image_data.shape[1]/2.0
            center = [center_x, center_y]
            smallest_dFoc = 1000000000
            radii = 21

            count = 0
            for i, j in zip(objects['x'], objects['y']):
                pos = [i, j]

                dFoc = math.sqrt(((pos[0]-center[0])**2) + ((pos[1]-center[1])**2))
                if abs(dFoc) < smallest_dFoc:
                    smallest_dFoc = dFoc
                    minx = objects['xmin'][count]
                    miny = objects['ymin'][count]
                    maxx = objects['xmax'][count]
                    maxy = objects['ymax'][count]
                    radii = abs(math.sqrt(((maxx-minx)**2)+((maxy-miny)**2)))
                else:
                    pass
                count += 1
            i = 0
            while i < 25:
                theta = 0
                area = 0
                while theta <= 2*pi:
                    area += (((i+.1)**2)/2) - ((i**2)/2)
                    theta += .001
                flux, fluxerr, flag = sep.sum_circann(internal_image_data, internal_image_data.shape[0]/2, internal_image_data.shape[1]/2, i, i+an_width)
                metric = (flux - bkg.globalback)/bkg.globalrms
                metric /= area
                if i > 1:
                    if metric < smallest_metric:
                        smallest_metric = metric
                        annulus = [i, i+an_width]
                else:
                    smallest_metric = metric
                    annulus = [i, i+an_width]
                i += 1
            i = 1
            an = [math.ceil(x*6) for x in annulus]
            if radii > an[0]:
                radii = an[0] - 0.5
            if radii > 30:
                radii = 30
            else:
                pass
            return {'Ap': radii, 'InAn': an[0], 'OutAn': an[1]}

def find_objects(obs, segmentation_map=False):
    import sep

    noise=np.sqrt(1.0/obs.weight[0,0])
    return sep.extract(
        obs.image,
        SEP_THRESH,
        segmentation_map=segmentation_map,
        err=noise,
        **SEP_PARS
    )

def find_stars_on_data(data, verbose = 0, useDS9 = False):
	bkg = sep.Background(data)
	bkg.subfrom(data)
	thres = 1.5 * bkg.globalrms
	if verbose > 1:
		print('global average background: {0:.2f} rms: {1:.3f} threshold: {2:.3f}'.format(bkg.globalback, bkg.globalrms, thres))
	objects = sep.extract(data, thres)
	# order by flux
	if len(objects) == 0:
		return []
	return sorted(objects, cmp=lambda x,y: cmp(y['flux'],x['flux']))

def extract_sources(data):
	try:
		bkg = sep.Background(data)
	except ValueError:
		data = data.byteswap().newbyteorder()
		bkg = sep.Background(data)
	back = bkg.back()
	data_woback = data-back
	thresh = 3.0 * bkg.globalrms
	objects = sep.extract(data_woback, thresh)
	return objects

def test_extract_segmentation_map():

    # Get some background-subtracted test data:
    data = np.copy(image_data)
    bkg = sep.Background(data, bw=64, bh=64, fw=3, fh=3)
    bkg.subfrom(data)

    objects, segmap = sep.extract(data, 1.5*bkg.globalrms,
                                  segmentation_map=True)

    assert type(segmap) is np.ndarray
    assert segmap.shape == data.shape
    for i in range(len(objects)):
        assert objects["npix"][i] == (segmap == i+1).sum()

    def ImageObjDetection(self,filter):
        print 'SEP object detection on %s......'%filter
        sigma = 3.0
        threshold = sigma * self.bkgRMS[filter]
        objs = sep.extract(self.dataList[filter], threshold, minarea=20)
        # add fields for mag and magerr
        desc = np.dtype([('RA','float64'),('DEC','float64')])
        newObjs = np.zeros(objs.shape, dtype = objs.dtype.descr + desc.descr)
        for name in objs.dtype.names:
            newObjs[name] = objs[name]
        # add RA and DEC to objs
        newObjs['RA'], newObjs['DEC'] = self.images[filter].xy_to_rd(objs['x'],objs['y'])

        return newObjs