Python photutils.CircularAperture类代码示例

    def autocorr(self, ell_mask_scale=2, aperture_radius=5, annulus_width=4):
        # Compute 2D autocorrelation function

        try:
            masked_image = self.masked_image.copy() 
            if ell_mask_scale > 0:
                masked_image.mask |=  ~self.elliptical_mask(ell_mask_scale)
            masked_image = masked_image.filled(0.0)

            fft_imgae = np.fft.fft2(masked_image)
            acorr_image = np.fft.ifft2(fft_imgae * np.conjugate(fft_imgae)).real
            acorr_image = np.fft.ifftshift(acorr_image)

            ny, nx = masked_image.shape
            yy, xx = np.mgrid[:ny, :nx]

            circ = CircularAperture([nx // 2, ny // 2], r=aperture_radius)
            ann = CircularAnnulus([nx // 2, ny // 2],
                                  r_in=aperture_radius,
                                  r_out=aperture_radius + annulus_width)

            ann_mean = aperture_photometry(
                acorr_image, ann)['aperture_sum'][0] / ann.area()
            circ_mean = aperture_photometry(
                acorr_image, circ)['aperture_sum'][0] / circ.area()
        except:
            acorr_image = np.nan
            circ_mean = np.nan
            ann_mean = np.nan

        return acorr_image, circ_mean, ann_mean

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)

def measure_one_median_bg(image, center, aperRad, metric, nSig, apMethod='exact'):
    """Class methods are similar to regular functions.

    Note:
        Do not include the `self` parameter in the ``Args`` section.

    Args:
        param1: The first parameter.
        param2: The second parameter.

    Returns:
        True if successful, False otherwise.

    """
    
    aperture       = CircularAperture(center, aperRad)
    aperture       = aperture.to_mask(method=apMethod)[0]
    aperture       = aperture.to_image(image.shape).astype(bool)
    backgroundMask = ~aperture
    
    medFrame  = median(image[backgroundMask])
    madFrame  = std(image[backgroundMask])
    
    medianMask= abs(image - medFrame) < nSig*madFrame
    
    maskComb  = medianMask*backgroundMask
    
    return median(image[maskComb])

def measure_one_annular_bg(image, center, innerRad, outerRad, metric, apMethod='exact'):
    """Class methods are similar to regular functions.

    Note:
        Do not include the `self` parameter in the ``Args`` section.

    Args:
        param1: The first parameter.
        param2: The second parameter.

    Returns:
        True if successful, False otherwise.

    """
    
    innerAperture   = CircularAperture(center, innerRad)
    outerAperture   = CircularAperture(center, outerRad)
    
    inner_aper_mask = innerAperture.to_mask(method=apMethod)[0]
    inner_aper_mask = inner_aper_mask.to_image(image.shape).astype(bool)
    
    outer_aper_mask = outerAperture.to_mask(method=apMethod)[0]
    outer_aper_mask = outer_aper_mask.to_image(image.shape).astype(bool)
    
    backgroundMask = (~inner_aper_mask)*outer_aper_mask
    
    return metric(image[backgroundMask])

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()

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)

    def show_image(self,sourceRA,sourceDEC,refRA=None,refDEC=None, ref2RA=None, ref2DEC=None):
        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=10)
       
        #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='jet')
        apertures.plot(color='red',lw=5, alpha=1)
        
        if refRA != None:
           ref1aper  = CircularAperture((self.worldcoord.wcs_world2pix(refRA,refDEC,0)),     r=10)
           ref1aper.plot(color='red', lw=5, alpha=1)
        #apertures2.plot(color='white',lw=2.0,alpha=0.5)
        #largeaperture.plot(color='red',  lw=1.5, alpha=0.5)
        if ref2RA != None:
           ref2aper  = CircularAperture((self.worldcoord.wcs_world2pix(ref2RA,ref2DEC,0)),     r=10)
           ref2aper.plot(color='red', lw=5, alpha=1)

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()

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]

 def phot(self, image, objpos, aper):
     """
     Aperture photometry using Astropy's photutils.
     
     Parameters
     ----------
     image : numpy array
         2D image array
         
     objpos : list of tuple
         Object poistions as list of tuples
         
     aper : float
         Aperture radius in pixels
      
     Returns 
     -------
     phot_table : astropy table
          Output table with stellar photometry   
     """
     try:
         from astropy.table import hstack
         from photutils import aperture_photometry, CircularAnnulus, CircularAperture
     except ImportError:
         pass
 
     apertures = CircularAperture(objpos, r = aper) 
     annulus_apertures = CircularAnnulus(objpos, r_in = self.inner_radius, r_out = self.outer_radius)
     
     rawflux_table = aperture_photometry(image, apertures = apertures, method = self.method)
     bkgflux_table = aperture_photometry(image, apertures = annulus_apertures, method = self.method)
     phot_table = hstack([rawflux_table, bkgflux_table], table_names = ["raw", "bkg"])
     
     bkg = phot_table["aperture_sum_bkg"] / annulus_apertures.area()
     phot_table["msky"] = bkg
     phot_table["area"] = apertures.area()
     phot_table["nsky"] = annulus_apertures.area()
             
     bkg_sum = bkg * apertures.area()
     final_sum = phot_table["aperture_sum_raw"] - bkg_sum
     phot_table["flux"] = final_sum
     
     return phot_table

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()

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

def calc_bkg_rms(ap, image, src_ap_area, rpsrc, mask=None, min_ap=6):

    if isinstance(ap, CircularAnnulus):
        aback = bback =  ap.r_in + rpsrc
        ap_theta = 0
    elif isinstance(ap, EllipticalAnnulus):
        aback = ap.a_in + rpsrc
        bback = ap.b_in + rpsrc
        ap_theta = ap.theta
    
    ecirc = ellip_circumference(aback, bback)
    diam = 2*rpsrc
    
    # Estimate the number of background apertures that can fit around the source
    # aperture.
    naps = np.int(np.round(ecirc/diam))
    
    # Use a minimum of 6 apertures
    naps = np.max([naps, min_ap])
    #naps = 6
    
    theta_back = np.linspace(0, 2*np.pi, naps, endpoint=False)
    
    # Get the x, y positions of the background apertures
    x, y = ellip_point(ap.positions[0], aback, bback, ap_theta, theta_back)

    # Create the background apertures and calculate flux within each
    bkg_aps = CircularAperture(np.vstack([x,y]).T, rpsrc)
    flux_bkg = aperture_photometry(image, bkg_aps, mask=mask)
    flux_bkg = flux_bkg['aperture_sum']
    flux_bkg_adj = flux_bkg/bkg_aps.area() * src_ap_area
 	
	# Use sigma-clipping to determine the RMS of the background
	# Scale to the area of the source aperture
    me, md, sd = sigma_clipped_stats(flux_bkg_adj, sigma=3)


    bkg_rms = sd
    
    return bkg_rms, bkg_aps

    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)')

def do_phot(image,position,radius = 5, r_in=15., r_out=20.):
    
    aperture = CircularAperture(position,r=radius)

    bkg_aperture = CircularAnnulus(position,r_in=r_in,r_out=r_out)

    # perform the photometry; the default method is 'exact'
    phot = aperture_photometry(image, aperture)
    bkg = aperture_photometry(image, bkg_aperture)

    # calculate the mean background level (per pixel) in the annuli
    bkg_mean = bkg['aperture_sum'] / bkg_aperture.area()
    bkg_sum = bkg_mean * aperture.area()
   
    #look at ipython notebook; his may need editing; 'phot' in second line below may need brackets with 'flux_sum' inside
    #phot['bkg_sum'] = bkg_sum
    #phot['flux_sum'] = phot['flux'] - bkg_sum
    
    #these two lines from ipython notebook
    flux_bkgsub = phot['aperture_sum'] - bkg_sum
    phot['aperture_sum_bkgsub'] = flux_bkgsub
    
    return phot