Python CircularAperture.plot方法代码示例

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
def plot_peaks(box, x_peaks, y_peaks, radius=None, title=None, vmin=None, vmax=None):

    """
    This function plots the data with peaks marked ...
    :param box:
    :param x_peaks:
    :param y_peaks:
    :return:
    """

    # Determine the maximum value in the box and the minium value for plotting
    if vmin is None: vmin = max(np.nanmin(box), 0.)
    if vmax is None: vmax = 0.5 * (np.nanmax(box) + vmin)

    # Set the normalization
    norm = ImageNormalize(stretch=SqrtStretch())

    # Make the plot
    plt.figure(figsize=(8,2.5))
    plt.imshow(box, origin='lower', norm=norm, interpolation='nearest', vmin=vmin, vmax=vmax, cmap="viridis")

    if radius is None: plt.plot(x_peaks, y_peaks, ls='none', color='white', marker='+', ms=40, lw=10, mew=4)
    else:

        positions = (x_peaks, y_peaks)
        apertures = CircularAperture(positions, r=radius)
        apertures.plot(color='green', lw=1.5, alpha=0.5)

    plt.xlim(0, box.shape[1]-1)
    plt.ylim(0, box.shape[0]-1)

    if title is not None: plt.title(title)

    plt.show()

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
def plot_light_curves(diff_cube, unique_extracted_objects):
    
    frame_data = [i for i in range(len(diff_cube))]
    colors = [(random.uniform(0.5, 1),random.uniform(0.5, 1),random.uniform(0.5,1)) for i in range(len(unique_extracted_objects))]


    plt.figure(figsize=(20,10))


    for i, extracted_obj in enumerate(unique_extracted_objects):
        ap_data=[]
        plt.figure(i, figsize=(10, 12))
        for frame in diff_cube:
            diff_cube_test = frame.copy()
            flux, fluxerr, flag = sep.sum_ellipse(diff_cube_test, x=extracted_obj[0], y=extracted_obj[1], a=extracted_obj[2], b=extracted_obj[3], theta=extracted_obj[4])
            #flux /= diff_cube_test.sum()
            ap_data.append(flux)

  
        plt.ylim((0,800))
        plt.plot(frame_data, ap_data, '-o', color=colors[i],linewidth=5.0, )

    plt.show()

    plt.figure(2, figsize=(10, 12))

    plt.imshow(diff_cube[1], cmap='gray', vmin=1, vmax=12)
    plt.colorbar()

    for i, extracted_obj in enumerate(unique_extracted_objects):
        positions = (extracted_obj[0], extracted_obj[1])
        apertures = CircularAperture(positions, r=5.)
        apertures.plot(color=colors[i], linewidth=10.0, lw=2.5, alpha=0.5)

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
def plot(sources, data, path):
    positions = (sources['xcentroid'], sources['ycentroid'])
    apertures = CircularAperture(positions, r=4.)
    norm = ImageNormalize(stretch=SqrtStretch())
    plt.imshow(data, cmap='Greys', origin='lower', norm=norm)
    apertures.plot(color='blue', lw=1.5, alpha=0.5)

    plt.savefig(path)

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
def plot_light_curves(diff_cube, unique_extracted_objects):
    # The diff_cube has to be byte swapped BEFORE being sent as parameter (diff_cube.byteswap(True).newbyteorder()), otherwise the method is not goint to work. Unique_extracted_objects only work for elliptic-shapped apertures

    # We get the number of frames from the cube
    frame_data = [i for i in range(len(diff_cube))]
    # Random colours array
    colours = [
        (random.uniform(0.5, 1), random.uniform(0.5, 1), random.uniform(0.5, 1))
        for i in range(len(unique_extracted_objects))
    ]

    maxVal = 0
    minVal = float("inf")

    plt.figure(2, figsize=(10, 12))

    # Bonus: Show the image with the sources on the same colour than the plots.
    if len(diff_cube) == 1:
        plt.imshow(diff_cube[0], cmap="gray", vmin=1, vmax=12)
    else:
        plt.imshow(diff_cube[1], cmap="gray", vmin=1, vmax=12)
    plt.colorbar()
    for i, extracted_obj in enumerate(unique_extracted_objects):
        positions = (extracted_obj[0], extracted_obj[1])
        apertures = CircularAperture(positions, r=5.0)
        apertures.plot(color=colours[i], linewidth=10.0, lw=2.5, alpha=0.5)
    # For every object we are going to calculate the aperture
    plt.figure(1, figsize=(20, 12))
    for i, extracted_obj in enumerate(unique_extracted_objects):
        ap_data = []
        # The standard size of each independent figure
        # plt.figure(i, figsize=(10, 12))
        # For every frame...
        for frame in diff_cube:
            diff_cube_test = frame.copy()
            # The parameters passed in order are x, y, a, b and theta
            flux, fluxerr, flag = sep.sum_ellipse(
                diff_cube_test,
                x=extracted_obj[0],
                y=extracted_obj[1],
                a=extracted_obj[2],
                b=extracted_obj[3],
                theta=extracted_obj[4],
            )

            ap_data.append(flux)
        maxVal = np.maximum(maxVal, np.max(ap_data))
        minVal = np.minimum(minVal, np.min(ap_data))
        # Hard-coded value!!! ALERT!!!

        # Plot every curve as a dotted line with the points visible
        plt.plot(frame_data, ap_data, "-o", color=colours[i], linewidth=5.0)
    plt.ylim((minVal * 1.1, maxVal * 0.9))
    # Voila
    plt.show()

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
    def show_image(self,sourceRA,sourceDEC,refRA,refDEC):
        print('RED circle is the cepheid. WHITE circle is the reference object(s).')
        print('Add more reference stars by defining ref1aper = blahblahbelow, and ref1aper.plot(etc...)')
        #aper_annulus = CircularAnnulus((sourceRA, sourceDEC), r_in=6., r_out = 8.)

        apertures = CircularAperture((self.worldcoord.wcs_world2pix(sourceRA,sourceDEC,0)), r=6)
        ref1aper  = CircularAperture((self.worldcoord.wcs_world2pix(refRA,refDEC,0)),     r=6)
        #ref2aper  = CircularAperture((worldcoord.wcs_world2pix(ref2RA,ref2DEC,0)),     r=7)
        #ref3aper  = CircularAperture((worldcoord.wcs_world2pix(ref3RA,ref3DEC,0)),     r=3.5)
        #darkaper  = CircularAperture((worldcoord.wcs_world2pix(darkRA,darkDEC,0)),     r=3.5)
        
        
        fig = plt.figure()
        fig.add_subplot(111, projection = self.worldcoord)
        plt.imshow(self.stardata,origin='lower', cmap='Spectral')
        apertures.plot(color='red',lw=1.5, alpha=0.5)
        ref1aper.plot(color='white', lw=2.5, alpha=0.5)

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
def starbright(fnstar,fnflat,istar,axs,fg):
    #%% load data
    data = meanstack(fnstar,100)[0]
    #%% flat field
    flatnorm = readflat(fnflat,fnstar)
    data = (data/flatnorm).round().astype(data.dtype)
    #%% background
    mean, median, std = sigma_clipped_stats(data, sigma=3.0)

    rfact=data.shape[0]//40
    cfact=data.shape[1]//40
    bg = Background(data,(rfact,cfact),interp_order=1, sigclip_sigma=3)
# http://docs.astropy.org/en/stable/units/#module-astropy.units
    #dataphot = (data - bg.background)*u.ph/(1e-4*u.m**2 * u.s * u.sr)
 #   data = (data-0.97*data.min()/bg.background.min()*bg.background) * u.ph/(u.cm**2 * u.s * u.sr)
    data = data* u.ph/(u.cm**2 * u.s * u.sr)
    #%% source extraction
    sources = daofind(data, fwhm=3.0, threshold=5*std)
    #%% star identification and quantification
    XY = column_stack((sources['xcentroid'], sources['ycentroid']))
    apertures = CircularAperture(XY, r=4.)
    norm = ImageNormalize(stretch=SqrtStretch())

    flux = apertures.do_photometry(data,effective_gain=camgain)[0]
#%% plots
    fg.suptitle('{}'.format(fnflat.parent),fontsize='x-large')

    hi = axs[-3].imshow(flatnorm,interpolation='none',origin='lower')
    fg.colorbar(hi,ax=axs[-3])
    axs[-3].set_title('flatfield {}'.format(fnflat.name))

    hi = axs[-2].imshow(bg.background,interpolation='none',origin='lower')
    fg.colorbar(hi,ax=axs[-2])
    axs[-2].set_title('background {}'.format(fnstar.name))

    hi = axs[-1].imshow(data.value,
                    cmap='Greys', origin='lower', norm=norm,interpolation='none')
    fg.colorbar(hi,ax=axs[-1])
    for i,xy in enumerate(XY):
        axs[-1].text(xy[0],xy[1], str(i),ha='center',va='center',fontsize=16,color='w')
    apertures.plot(ax=axs[-1], color='blue', lw=1.5, alpha=0.5)
    axs[-1].set_title('star {}'.format(fnstar.name))

    return flux[istar]

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
def process_file(inpath, file_name, t_constant, sigma, fwhm, r, kernel_size, outpath, plot):
    print "Processing " + file_name
    hdulist = fits.open(inpath + file_name)
    image = hdulist[0].data

    if isinstance(sigma, list):
        threshold = calc_sigma(image, sigma[0], sigma[1]) * t_constant
    else:
        threshold = t_constant*sigma

    median_out = signal.medfilt(image,kernel_size)
    median_sub = np.subtract(image,median_out)
    sources = daofind(median_sub, threshold, fwhm)

    sources_2 = np.array(sources["id", "xcentroid", "ycentroid", "sharpness", "roundness1", "roundness2", "npix", "sky", "peak", "flux", "mag"])
    print_line= (file_name+","+str(sources_2))
    base_name = os.path.splitext(file_name)[0]
    file = open(outpath + base_name + ".out", "a")
    file.write(print_line)
    file.close()

    positions = (sources['xcentroid'], sources['ycentroid'])
#    print positions
    apertures = CircularAperture(positions, r)
    phot_table = aperture_photometry(median_sub, apertures)
    phot_table_2 = np.array(phot_table["aperture_sum", "xcenter", "ycenter"])
    print_line= (","+str(phot_table_2)+"\n")
    file = open(outpath + base_name + ".out", "a")
    file.write(print_line)
    file.write("\n")
    file.close()

    hdulist[0].data = median_sub
    file = open(outpath + base_name + ".fits", "w")
    hdulist.writeto(file)
    file.close()

    if plot:
        median_sub[median_sub<=0]=0.0001
        plt.imshow(median_sub, cmap='gray', origin='lower')
        apertures.plot(color='blue', lw=1.5, alpha=0.5)
        plt.show()

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
def find_stars(image, plot = False, fwhm = 20.0, threshold=3.):

    from astropy.stats import sigma_clipped_stats
    mean, median, std = sigma_clipped_stats(image, sigma=3.0)
    from photutils import daofind
    sources = daofind(image - median, fwhm=fwhm, threshold=threshold*std)
    
   # stars already found accurately, vet_sources will be implemented when working properly
   # vet_sources(10.0,10.0)
        
    if plot == True:
       # from astropy.visualization import SqrtStretch
       # from astropy.visualization.mpl_normalize import ImageNormalize
        positions = (sources['xcentroid'], sources['ycentroid'])
        apertures = CircularAperture(positions, r=4.)
        #norm = ImageNormalize(stretch=SqrtStretch())
        #plt.imshow(image, cmap='Greys', origin='lower', norm=norm)
        qi.display_image(image)
        apertures.plot(color='blue', lw=1.5, alpha=0.5)
        
    return sources

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
    def plot(self,scale='log'):
        apertures = CircularAperture([self.locx,self.locy], r=self.r)

        z = self.copy()
        z -= np.nanmedian(z)
        if scale=='log':
            z = np.log10(self)
            z = ma.masked_invalid(z)
            z.mask = z.mask | (z < 0)
            z.fill_value = 0
            z = z.filled()

        imshow2(z)

        if self.pixels is not None:
            for i,pos in enumerate(self.pixels):
                r,c = pos
                plt.text(c,r,i,va='center',ha='center',color='Orange')

        apertures.plot(color='Lime',lw=1.5,alpha=0.5)
        plt.xlabel('Column (pixels)')
        plt.ylabel('Row (pixels)')

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
def try_dao(fitsfile, outfile='out.png'):
   
    hdulist = fits.open(fitsfile)
    hdu = hdulist[0]

    if len(hdu.data.shape) == 3:
        data = hdu.data[0]
    else:
        data = hdu.data

    mean, median, std = sigma_clipped_stats(data[800:900, 800:900], sigma=3.0, iters=5)
    print(mean, median, std)

    #data = hdu.data
    sources = daofind(data - median, fwhm=3.0, threshold=5.*std)
    print 'Found %i sources' % len(sources)

    positions = (sources['xcentroid'], sources['ycentroid'])
    apertures = CircularAperture(positions, r=4.)
    norm = ImageNormalize(stretch=SqrtStretch(), vmin=2000, vmax=3000)
    plt.imshow(data, cmap='Greys', origin='lower', norm=norm)
    plt.title('%i Sources from a single all-sky frame' % len(sources))
    apertures.plot(color='blue', lw=1.5, alpha=0.5)
    plt.savefig(filename=outfile)

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
hdulist = fits.open(inpath+file_name)
image = hdulist[0].data
#image = image.astype(float) - np.median(image)
from photutils import daofind
from astropy.stats import mad_std
bkg_sigma = mad_std(image)
sources = daofind(image, fwhm, threshold*bkg_sigma)
#print_line= (file_name+","+str(sources_2)+"\n")
sources_2 = np.array(sources["id", "xcentroid", "ycentroid", "sharpness", "roundness1", "roundness2", "npix", "sky", "peak", "flux", "mag"])
print_line= (file_name+","+str(sources_2))
file= open(outpath, "a")
file.write(print_line)
file.close()

from photutils import aperture_photometry, CircularAperture
positions = (sources['xcentroid'], sources['ycentroid'])
apertures = CircularAperture(positions, r)
phot_table = aperture_photometry(image, apertures)
phot_table_2 = np.array(phot_table["aperture_sum", "xcenter", "ycenter"])
print_line= (","+str(phot_table_2)+"\n")
file= open(outpath, "a")
file.write(print_line)
file.close()

import matplotlib.pylab as plt
im2 = image
im2[im2<=0]=0.0001
plt.imshow(im2, cmap='gray', origin='lower')
apertures.plot(color='blue', lw=1.5, alpha=0.5)
plt.show()

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
def init_centroids(first_image_path, master_flat, master_dark, target_centroid,
                   max_number_stars=10, min_flux=0.2, plots=False):

    first_image = np.median([(fits.getdata(path) - master_dark)/master_flat
                             for path in first_image_path], axis=0)

    tophat_kernel = Tophat2DKernel(5)
    convolution = convolve_fft(first_image, tophat_kernel, fftn=fft2, ifftn=ifft2)

    convolution -= np.median(convolution)

    mad = mad_std(convolution)

    convolution[convolution < -5*mad] = 0.0

    from skimage.filters import threshold_yen
    from skimage.measure import label, regionprops

    thresh = threshold_yen(convolution)/4 # Use /4 for planet c, /2 for planet b
    #thresh = threshold_otsu(convolution)/15

    masked = np.ones_like(convolution)
    masked[convolution <= thresh] = 0

    label_image = label(masked)

    plt.figure()
    plt.imshow(label_image, origin='lower', cmap=plt.cm.viridis)
    plt.show()

    # regions = regionprops(label_image, convolution)
    regions = regionprops(label_image, first_image)

    # reject regions near to edge of detector
    buffer_pixels = 50
    regions = [region for region in regions
               if ((region.weighted_centroid[0] > buffer_pixels and
                   region.weighted_centroid[0] < label_image.shape[0] - buffer_pixels)
               and (region.weighted_centroid[1] > buffer_pixels and
                    region.weighted_centroid[1] < label_image.shape[1] - buffer_pixels))]

    #centroids = [region.weighted_centroid for region in regions]
    #intensities = [region.mean_intensity for region in regions]

    target_intensity = regions[0].mean_intensity
    target_diameter = regions[0].equivalent_diameter
    #  and region.equivalent_diameter > 0.8 * target_diameter
    centroids = [region.weighted_centroid for region in regions
                 if min_flux * target_intensity < region.mean_intensity]
    # intensities = [region.mean_intensity for region in regions
    #                if min_flux * target_intensity < region.mean_intensity]
#    centroids = np.array(centroids)[np.argsort(intensities)[::-1]]

    distances = [np.sqrt((target_centroid[0] - d[0])**2 +
                         (target_centroid[1] - d[1])**2) for d in centroids]

    centroids = np.array(centroids)[np.argsort(distances)]

    positions = np.vstack([[y for x, y in centroids], [x for x, y in centroids]])

    if plots:
        apertures = CircularAperture(positions, r=12.)
        apertures.plot(color='r', lw=2, alpha=1)
        plt.imshow(first_image, vmin=np.percentile(first_image, 0.01),
                   vmax=np.percentile(first_image, 99.9), cmap=plt.cm.viridis,
                   origin='lower')
        plt.scatter(positions[0, 0], positions[1, 0], s=150, marker='x')

        plt.show()
    return positions

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]

#.........这里部分代码省略.........
			
			fwhm = np.mean(fwhm)
		
		# Stack both ord and exord sources together
		tot_sources = vstack([sources_o,sources_e])
				
		# Store the ordinary and extraordinary beam source images and
		# create apertures for aperture photometry 
		positions = np.swapaxes(np.array((tot_sources['xcentroid'],
			tot_sources['ycentroid']),dtype='float'),0,1)
		aperture = CircularAperture(positions, r=0.5*apermul*fwhm)
		phot_table = aperture_photometry(image_data,aperture)   
					  
		# Set up arrays of ord and exord source parameters
		s_id = np.zeros([len(np.array(phot_table['id']))])
		xp = np.zeros([len(s_id)])
		yp = np.zeros([len(s_id)])
		fluxbgs = np.zeros([len(s_id)])
		mean_bg = np.zeros([len(s_id)])
		bg_err = np.zeros([len(s_id)])
		s_area = []
		ann_area = []
		
		for i in range(0,len(np.array(phot_table['id'])),1):
			s_id[i] = np.array(phot_table['id'][i])
			xpos = np.array(phot_table['xcenter'][i])
			ypos = np.array(phot_table['ycenter'][i])
			xp[i] = xpos
			yp[i] = ypos
			s_area.append(np.pi*(0.5*apermul*fwhm)**2)
			j = i%2				
			fluxbgs[i] = (phot_table['aperture_sum'][i] -
				aperture.area()*glob_bgm[j])
			mean_bg[i] = glob_bgm[j]
			bg_err[i] = glob_bgerr[j]
			ann_area.append(80*80)			
		
		# Create and save the image in z scale and overplot the ordinary and
		# extraordinary apertures and local background annuli if applicable
		fig = plt.figure()
		zscale = ZScaleInterval(image_data)
		norm = ImageNormalize(stretch=SqrtStretch(),interval=zscale)
		image = plt.imshow(image_data,cmap='gray',origin='lower',norm=norm)
		bg_annulus_o = RectangularAnnulus((843,159),w_in=0,w_out=80,h_out=80,
			theta=0)
		bg_annulus_e = RectangularAnnulus((843,69),w_in=0,w_out=80,h_out=80,
			theta=0)
		bg_annulus_o.plot(color='skyblue',lw=1.5,alpha=0.5)
		bg_annulus_e.plot(color='lightgreen',lw=1.5,alpha=0.5)
		
		for i in range(0,len(np.array(phot_table['id'])),1):
			aperture = CircularAperture((xp[i],yp[i]),r=0.5*apermul*fwhm)
			
			if i < int(len(np.array(phot_table['id']))/2):
				aperture.plot(color='blue',lw=1.5,alpha=0.5)
		
			else:
				aperture.plot(color='green',lw=1.5,alpha=0.5)
			
		plt.xlim(760,920)
		plt.ylim(20,210)
		plt.title(label[k])
		image_fn = folder_path + angle[k] + '_image.png'
		fig.savefig(image_fn)

		# Write ordinary and extraordinary beams to file following the 
		# convention angleXXX_ord.txt and angleXXX_exord.txt
		orig_stdout = sys.stdout
		ord_result_file= folder_path + 'angle' + angle[k] + '_ord.txt'
		ordresultf = open(ord_result_file, 'w')
		sys.stdout = ordresultf
		
		print("# id, xpix, ypix, fluxbgs, sourcearea, meanbg, bgerr, bgarea") 
		for i in range(0,int(len(np.array(phot_table['id']))/2),1):
			print(i+1,xp[i],yp[i],fluxbgs[i],s_area[i],mean_bg[i],bg_err[i],
				ann_area[i])
		sys.stdout = orig_stdout
		ordresultf.close()

		orig_stdout = sys.stdout
		exord_result_file = folder_path + 'angle' + angle[k] + '_exord.txt'
		exordresultf = open(exord_result_file, 'w')
		sys.stdout = exordresultf
		
		print("# id, xpix, ypix, fluxbgs, sourcearea, meanbg, bgerr, bgarea")
		for i in range(int(len(np.array(phot_table['id']))/2),len(np.array
			(phot_table['id'])),1):
			print(i+1-int(len(np.array(phot_table['id']))/2),xp[i],yp[i],
				fluxbgs[i],s_area[i],mean_bg[i],bg_err[i],ann_area[i])  
		sys.stdout = orig_stdout
		exordresultf.close()
		
		# Save the number of sources in each beam to a list
		numsource.append(int(len(np.array(phot_table['id']))/2))
	
	# Print number of sources per half-wave plate image
	for i in range(0,len(numsource),1):
		print("No of sources detected at",ang_dec[i],"degrees:",numsource[i])
	
	return 0

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
  def extension(self, extension_idx, threshold='', FWHM=3.0, sigma=3.0, snr=50., plot=False):
    '''
    A method to run aperatue photometry routines on an individual extension and save the results to the exposure class
    
    Parameters
    ----------
    extension_idx: int
      Index of the extension
    threshold: float (optional)
      The absolute image value above which to select sources
    FWHM: float
      The full width at half maximum
    sigma: float
      Number of standard deviations to use for background estimation
    snr: float
      The signal-to-noise ratio to use in the threshold detection
    plot: bool
      Plot the field with identified sources circled      

    Returns
    -------
    source_list: table
      A source list for the image

    '''

    # Define the data array
    data = self.hdulist[extension_idx].data.astype(np.float)
    
    # Extract the header and create a WCS object
    hdr = self.hdulist[extension_idx].header
    wcs = WCS(hdr)

    # Estimate the background and background noise
    mean, median, std = sigma_clipped_stats(data, sigma=sigma, iters=5)

    # Calculate the detection threshold and FWHM if not provided
    if not threshold: threshold = np.mean(detect_threshold(data, snr=snr))
    
    # Print the parameters being used
    for p,v in zip(['mean','median','std','threshold','FWHM'],[mean,median,std,threshold,FWHM]): print '{!s:10}: {:.3f}'.format(p,v)

    # Subtract background and generate sources list of all detections
    sources = daofind(data-median, threshold, FWHM)
    
    # Map RA and Dec to pixels
    positions = (sources['xcentroid'], sources['ycentroid'])
    skycoords = pixel_to_skycoord(*positions, wcs=wcs)
    
    # Calculate magnitudes at given source positions
    apertures = CircularAperture(positions, r=2.)
    photometry_table = aperture_photometry(data, apertures)
    
    # 'skycoords' IRCS object is problematic for stacking tables so for now we'll just add the ra and dec
    # photometry_table['sky_center'] = skycoords
    photometry_table['ra'], photometry_table['dec'] = skycoords.ra, skycoords.dec
    
    # Update data in the exposure object
    self.source_table = vstack([self.source_table,photometry_table], join_type='inner')  
    
    # Plot the sources
    if plot:
      norm = ImageNormalize(stretch=SqrtStretch())
      plt.imshow(data, cmap='Greys', origin='lower', norm=norm)
      apertures.plot(color='blue', lw=1.5, alpha=0.5)
    
    print '{!s:10}: {}'.format('sources',len(sources))

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import plot [as 别名]
groups = daogroup(intab, crit_separation=2.0*sigma_psf*gaussian_sigma_to_fwhm)

plt.subplot(1, 2, 1)
plt.imshow(image, origin='lower', interpolation='nearest')
plt.title('Simulated data')
plt.xlabel('x-position (pixel units)')
plt.ylabel('y-position (pixel units)')

plt.subplot(1, 2, 2)
for i in range(len(groups)):
    for j in range(len(groups[i]['id'])):
        # show ids
        # plt.text(groups[i]['x_0'][j], groups[i]['y_0'][j],
        #          str(groups[i]['id'][j]))
        aperture = CircularAperture((groups[i]['x_0'][j],
                                     groups[i]['y_0'][j]),
                                     r=sigma_psf*gaussian_sigma_to_fwhm)
        aperture.plot(lw=1.5, alpha=0.5)
tab, residual_image = nstar(image-bkg, groups, (5,5),
                            fitting.LevMarLSQFitter(),
                            IntegratedGaussianPRF, sigma=2.0)
tab.sort('id')
tab.write('starlist_estimated.html')

plt.imshow(residual_image, origin='lower', interpolation='nearest')
plt.title('Residual')
plt.xlabel('x-position (pixel units)')
plt.ylabel('y-position (pixel units)')
plt.show()