Python CircularAperture.area方法代码示例

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import area [as 别名]
    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

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import area [as 别名]
 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

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import area [as 别名]
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

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import area [as 别名]
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

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import area [as 别名]
def photom_av(ima, pos, radius, r_in=False, r_out=False, mode='median'):
    '''
    Aperture photometry in an aperture located at pixel coordinates 
    pos = ( (x0, y0), (x1, y1), ... ) with a radius=radius.
    When r_in and r_out are given, background is estimated in CircularAnnulus and subtracted.
    
    mode refers to how the background is estimated within the circlar annulus.
    Can be 'median' or 'mean'
    
    Photometry is calculating by median averaging the pixels within the aperture and 
    multiplying by the number of pixels in the aperture (including fractions of pixels).

    '''
    # Setting up the mask 
    if hasattr(ima, 'mask'):
      if ima.mask.size == 1:
	mask = np.zeros(ima.shape, dtype=np.bool) | ima.mask
      else:
        mask = ima.mask.copy()
    else:
        mask = np.zeros(ima.shape, dtype=np.bool)
        
        
    ### Performing the actual photometry - identical for each method
    # Median averaging of flux in aperture
    # Setting up the aperture 
    apertures = CircularAperture(pos, r = radius) 
    ap_mask = apertures.to_mask(method='center')
    # Setting up arrays to store data
    nflx = len(ap_mask)
    flx = np.zeros(nflx, dtype=np.float)
    flux_max = np.zeros(nflx, dtype=np.float)
    flux_min = np.zeros(nflx, dtype=np.float)
    # Median averaging of flux
    for i, am in enumerate(ap_mask):
      fluxmask = ~mask & am.to_image(shape=mask.shape).astype(np.bool)
      flx[i] = np.median(ima[fluxmask])
      flux_max[i] = np.max(ima[fluxmask])
      flux_min[i] = np.min(ima[fluxmask])
      
      
      
      
    # Aperture photometry on mask to see how many masked pixels are in the 
    # aperture
    apm       = aperture_photometry(mask.astype(int), apertures)
    # Number of unmasked pixels in aperture
    ap_area   = Column(name = 'area_aper',
		       data=apertures.area() - apm['aperture_sum'].data)
    
    # Flux in aperture using median av flux and fractional no. pixels in aperture
    flux_init = flx*ap_area
    

    ### Two different modes for analysing the background
    if ( r_in and r_out and mode in ('mean', 'median') ):
      
      ### This stuff is the same regardless of method
      # Setting up the annulus
      anulus_apertures = CircularAnnulus(pos, r_in=r_in, r_out=r_out)
      # Performing annulus photometry on the mask
      bkgm = aperture_photometry(mask.astype(int), anulus_apertures)
      # Number of masked pixels in bkg
      mbkg_area = Column(name = 'bpix_bkg',
			 data=bkgm['aperture_sum'])  
      # Number of non-masked pixels in aperture and bkg        
      bkg_area  = Column(name = 'area_bkg',
			 data=anulus_apertures.area() - bkgm['aperture_sum'])
      
      
      ### This stuff is specific to the mean
      if mode == 'mean':
	# Perform the annulus photometry on the image
	bkg  = aperture_photometry(ima, anulus_apertures, mask=mask)
        # Average bkg where this divides by only number of NONMASKED pixels
        # as the aperture photometry ignores the masked pixels
        bkga = Column(name='background',
		      data=bkg['aperture_sum']/bkg_area)
        # Bkg subtracted flux
        flux = flux_init - bkga*ap_area
        # Adding that data
        ap.add_column(bkga)
        
        
      elif mode == 'median':
	# Number of pixels in the annulus, a different method
	aperture_mask = anulus_apertures.to_mask(method='center')
	nbkg = len(aperture_mask)
	
	# Background mask
	bkgm = np.zeros(nbkg, dtype=np.float)
	
	# Median averaging
	for i, am in enumerate(aperture_mask):
	  bmask = ~mask & am.to_image(shape=mask.shape).astype(np.bool)
	  bkgm[i] = np.median(ima[bmask])
		
	flux = flux_init - bkgm*ap_area
	bkgm = Column(name = 'background', data = bkgm)

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

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import area [as 别名]
def tso_aperture_photometry(datamodel, xcenter, ycenter, radius, radius_inner,
                            radius_outer):
    """
    Create a photometric catalog for NIRCam TSO imaging observations.

    Parameters
    ----------
    datamodel : `CubeModel`
        The input `CubeModel` of a NIRCam TSO imaging observation.

    xcenter, ycenter : float
        The ``x`` and ``y`` center of the aperture.

    radius : float
        The radius (in pixels) of the circular aperture.

    radius_inner, radius_outer : float
        The inner and outer radii (in pixels) of the circular-annulus
        aperture, used for local background estimation.

    Returns
    -------
    catalog : `~astropy.table.QTable`
        An astropy QTable (Quantity Table) containing the source
        photometry.
    """

    if not isinstance(datamodel, CubeModel):
        raise ValueError('The input data model must be a CubeModel.')

    aper1 = CircularAperture((xcenter, ycenter), r=radius)
    aper2 = CircularAnnulus((xcenter, ycenter), r_in=radius_inner,
                            r_out=radius_outer)

    nimg = datamodel.data.shape[0]
    aperture_sum = []
    aperture_sum_err = []
    annulus_sum = []
    annulus_sum_err = []

    for i in np.arange(nimg):
        tbl1 = aperture_photometry(datamodel.data[i, :, :], aper1,
                                   error=datamodel.err[i, :, :])
        tbl2 = aperture_photometry(datamodel.data[i, :, :], aper2,
                                   error=datamodel.err[i, :, :])

        aperture_sum.append(tbl1['aperture_sum'][0])
        aperture_sum_err.append(tbl1['aperture_sum_err'][0])
        annulus_sum.append(tbl2['aperture_sum'][0])
        annulus_sum_err.append(tbl2['aperture_sum_err'][0])

    # convert array of Quantities to Quantity arrays
    aperture_sum = u.Quantity(aperture_sum)
    aperture_sum_err = u.Quantity(aperture_sum_err)
    annulus_sum = u.Quantity(annulus_sum)
    annulus_sum_err = u.Quantity(annulus_sum_err)

    # construct metadata for output table
    meta = OrderedDict()
    meta['instrument'] = datamodel.meta.instrument.name
    meta['detector'] = datamodel.meta.instrument.detector
    meta['channel'] = datamodel.meta.instrument.channel
    meta['subarray'] = datamodel.meta.subarray.name
    meta['filter'] = datamodel.meta.instrument.filter
    meta['pupil'] = datamodel.meta.instrument.pupil

    meta['target_name'] = datamodel.meta.target.catalog_name
    meta['xcenter'] = xcenter
    meta['ycenter'] = ycenter
    ra_icrs, dec_icrs = datamodel.meta.wcs(xcenter, ycenter)
    meta['ra_icrs'] = ra_icrs
    meta['dec_icrs'] = dec_icrs

    info = ('Photometry measured in a circular aperture of r={0} pixels. '
            'Background calculated as the mean in a circular annulus with '
            'r_inner={1} pixels and r_outer={2} pixels.'
            .format(radius, radius_inner, radius_outer))
    meta['apertures'] = info

    tbl = QTable(meta=meta)

    dt = (datamodel.meta.exposure.group_time *
          (datamodel.meta.exposure.ngroups + 1))
    dt_arr = (np.arange(1, 1 + datamodel.meta.exposure.nints) *
              dt - (dt / 2.))
    int_dt = TimeDelta(dt_arr, format='sec')
    int_times = (Time(datamodel.meta.exposure.start_time, format='mjd') +
                 int_dt)
    tbl['MJD'] = int_times.mjd

    tbl['aperture_sum'] = aperture_sum
    tbl['aperture_sum_err'] = aperture_sum_err
    tbl['annulus_sum'] = annulus_sum
    tbl['annulus_sum_err'] = annulus_sum_err

    annulus_mean = annulus_sum / aper2.area()
    annulus_mean_err = annulus_sum_err / aper2.area()
    tbl['annulus_mean'] = annulus_mean
    tbl['annulus_mean_err'] = annulus_mean_err

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

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

#.........这里部分代码省略.........
    #                                target_centroid, plots=True,
    #                                min_flux=comparison_flux_threshold).T

    # Initialize some empty arrays to fill with data:
    times = np.zeros(len(image_paths))
    fluxes = np.zeros((len(image_paths), len(star_positions),
                       len(aperture_radii)))
    errors = np.zeros((len(image_paths), len(star_positions),
                       len(aperture_radii)))
    xcentroids = np.zeros((len(image_paths), len(star_positions)))
    ycentroids = np.zeros((len(image_paths), len(star_positions)))
    airmass = np.zeros(len(image_paths))
    airpress = np.zeros(len(image_paths))
    humidity = np.zeros(len(image_paths))
    telfocus = np.zeros(len(image_paths))
    psf_stddev = np.zeros(len(image_paths))

    medians = np.zeros(len(image_paths))

    with ProgressBar(len(image_paths)) as bar:
        for i in range(len(image_paths)):
            bar.update()

            # Subtract image by the dark frame, normalize by flat field
            imagedata = (fits.getdata(image_paths[i]) - master_dark) / master_flat

            # Collect information from the header
            imageheader = fits.getheader(image_paths[i])
            exposure_duration = imageheader['EXPTIME']
            times[i] = Time(imageheader['DATE-OBS'], format='isot', scale=imageheader['TIMESYS'].lower()).jd
            medians[i] = np.median(imagedata)
            airmass[i] = imageheader['AIRMASS']
            airpress[i] = imageheader['AIRPRESS']
            humidity[i] = imageheader['HUMIDITY']
            telfocus[i] = imageheader['TELFOCUS']

            # Initial guess for each stellar centroid informed by previous centroid
            for j in range(len(star_positions)):
                if i == 0:
                    init_x = star_positions[j][0]
                    init_y = star_positions[j][1]
                else:
                    init_x = ycentroids[i-1][j]
                    init_y = xcentroids[i-1][j]

                # Cut out a stamp of the full image centered on the star
                image_stamp = imagedata[int(init_y) - centroid_stamp_half_width:
                                        int(init_y) + centroid_stamp_half_width,
                                        int(init_x) - centroid_stamp_half_width:
                                        int(init_x) + centroid_stamp_half_width]

                # Measure stellar centroid with 2D gaussian fit
                x_stamp_centroid, y_stamp_centroid = centroid_com(image_stamp)
                y_centroid = x_stamp_centroid + init_x - centroid_stamp_half_width
                x_centroid = y_stamp_centroid + init_y - centroid_stamp_half_width

                xcentroids[i, j] = x_centroid
                ycentroids[i, j] = y_centroid

                # For the target star, measure PSF:
                if j == 0:
                    psf_model_init = models.Gaussian2D(amplitude=np.max(image_stamp),
                                                       x_mean=centroid_stamp_half_width,
                                                       y_mean=centroid_stamp_half_width,
                                                       x_stddev=psf_stddev_init,
                                                       y_stddev=psf_stddev_init)

                    fit_p = fitting.LevMarLSQFitter()
                    y, x = np.mgrid[:image_stamp.shape[0], :image_stamp.shape[1]]
                    best_psf_model = fit_p(psf_model_init, x, y, image_stamp -
                                           np.median(image_stamp))
                    psf_stddev[i] = 0.5*(best_psf_model.x_stddev.value +
                                          best_psf_model.y_stddev.value)

            positions = np.vstack([ycentroids[i, :], xcentroids[i, :]])

            for k, aperture_radius in enumerate(aperture_radii):
                target_apertures = CircularAperture(positions, aperture_radius)
                background_annuli = CircularAnnulus(positions,
                                                    r_in=aperture_radius +
                                                         aperture_annulus_radius,
                                                    r_out=aperture_radius +
                                                          2 * aperture_annulus_radius)
                flux_in_annuli = aperture_photometry(imagedata,
                                                     background_annuli)['aperture_sum'].data
                background = flux_in_annuli/background_annuli.area()
                flux = aperture_photometry(imagedata,
                                           target_apertures)['aperture_sum'].data
                background_subtracted_flux = (flux - background *
                                              target_apertures.area())

                fluxes[i, :, k] = background_subtracted_flux/exposure_duration
                errors[i, :, k] = np.sqrt(flux)

    ## Save some values
    results = PhotometryResults(times, fluxes, errors, xcentroids, ycentroids,
                                airmass, airpress, humidity, medians,
                                psf_stddev, aperture_radii)
    results.save(output_path)
    return results

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import area [as 别名]
def photom(ima, pos, radius, r_in=None, r_out=None, method='median'):
    '''
    Aperture photometry in an aperture located at pixel coordinates 
    pos = ( (x0, y0), (x1, y1), ... ) with a radius=radius.
    When r_in and r_out are given, background is estimated in CircularAnnulus and subtracted.
    
    method refers to how the background is estimated within the circlar annulus.
    Can be 'median' or 'mean' or 'mode'

    '''
    
    ima_local = np.ma.asanyarray(ima.copy())
    ima_local.fill_value = np.nan
    mask_ = ima_local.mask
    ima_  = ima_local.filled()
    
    ### Do photometry - identical for each method
    apertures = CircularAperture(pos, r = radius)
    ap        = aperture_photometry(ima_, apertures, 
                                    mask=mask_, method='exact')
    # Aperture photometry on mask to estimate # of masked pixels in aperture
    apm       = aperture_photometry(mask_.astype(int), apertures,
                                    method='exact')
    ap.add_columns( [apertures.area()-apm['aperture_sum'], apm['aperture_sum']],
                   names=['aperture_area', 'aperture_badpix'])
    ap.add_column(ap['aperture_sum'], index=3, name='Flux')

    if ( r_in == None or r_out == None or not method in ('mean', 'median', 'mode') ): 
      # Quit here if background correction is not requested
      return ap

    annulus_apertures = CircularAnnulus(pos, r_in=r_in, r_out=r_out)
    annulus_masks = annulus_apertures.to_mask(method='center')
    bg_values = []
    for annulus_mask in annulus_masks:
      bg_ima = annulus_mask.cutout(ima_)
      bg_mask = annulus_mask.cutout(mask_.astype(np.int))
      bg_ima = np.ma.array(bg_ima, mask= bg_mask.astype(np.bool) | ~annulus_mask.data.astype(np.bool))
      
      if method == 'mean': bg_val = bg_ima.mean()
      elif method == 'median': bg_val = np.ma.median(bg_ima)
      elif method == 'mode': 
        kernel = gaussian_kde(bg_ima.data[~bg_ima.mask], bw_method='silverman')
        mode = bg_ima.mean()
        std  = bg_ima.std()
        
        mode = minimize_scalar(lambda x: -kernel(x), bounds=(mode-3*std, mode+3*std),
                               method='bounded')
        bg_val=mode.x[0]
        
        if False:
          median = np.ma.median(bg_ima)
          h, b = np.histogram(bg_ima.data[~bg_ima.mask], bins=15, normed=True)
          bc = 0.5*(b[1:]+ b[:-1])
          plt.figure(33); plt.clf(); plt.ioff()
          fig, (ax0,ax1) = plt.subplots(ncols=2, nrows=1, num=33)
          ax0.plot(bc, h, 'x')
          x = np.linspace(bc.min(), bc.max(), 100)
          ax0.plot(x, kernel(x))
          ax0.vlines(mode.x, ax0.get_ylim()[0], ax0.get_ylim()[1])
          ax0.vlines(median, ax0.get_ylim()[0], ax0.get_ylim()[1])
          ax1.imshow(bg_ima)
          plt.show()
        
        
      bg_values.append(bg_val)
    ap.add_column(Column(data=bg_values, name = 'background'))  
    ap['Flux'] = ap['Flux'] - ap['aperture_area']*ap['background']
    return ap, bg_ima

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

#.........这里部分代码省略.........
        if method == "subpixel":
            log.debug("  subpixels = %s", str(subpixels))
    else:
        log.debug("  width = %s", str(width))
        log.debug("  height = %s", str(height))
        log.debug("  theta = %s degrees", str(extract_params['theta']))
        log.debug("  subtract_background = %s", str(subtract_background))
        log.debug("  method = %s", method)
        if method == "subpixel":
            log.debug("  subpixels = %s", str(subpixels))

    # Check for out of bounds.
    # The problem with having the background aperture extend beyond the
    # image is that the normalization would not account for the resulting
    # decrease in the area of the annulus, so the background subtraction
    # would be systematically low.
    outside = False
    f_nx = float(shape[2])
    f_ny = float(shape[1])
    if x_center < 0. or x_center >= f_nx - 1. or \
       y_center < 0. or y_center >= f_ny - 1.:
        outside = True
        log.error("Target location is outside the image.")
    if subtract_background and \
       (x_center - outer_bkg < -0.5 or x_center + outer_bkg > f_nx - 0.5 or
        y_center - outer_bkg < -0.5 or y_center + outer_bkg > f_ny - 0.5):
            outside = True
            log.error("Background region extends outside the image.")

    if outside:
        (ra, dec) = (0., 0.)
        wavelength = np.zeros(shape[0], dtype=np.float64)
        dq[:] = dqflags.pixel['DO_NOT_USE']
        return (ra, dec, wavelength, net, background, dq)       # all bad

    if hasattr(input_model.meta, 'wcs'):
        wcs = input_model.meta.wcs
    else:
        log.warning("WCS function not found in input.")
        wcs = None

    if wcs is not None:
        x_array = np.empty(shape[0], dtype=np.float64)
        x_array.fill(float(shape[2]) / 2.)
        y_array = np.empty(shape[0], dtype=np.float64)
        y_array.fill(float(shape[1]) / 2.)
        z_array = np.arange(shape[0], dtype=np.float64) # for wavelengths
        ra, dec, wavelength = wcs(x_array, y_array, z_array)
        nelem = len(wavelength)
        ra = ra[nelem // 2]
        dec = dec[nelem // 2]
    else:
        (ra, dec) = (0., 0.)
        wavelength = np.arange(1, shape[0] + 1, dtype=np.float64)

    position = (x_center, y_center)
    if source_type == 'point':
        aperture = CircularAperture(position, r=radius)
        if subtract_background:
            annulus = CircularAnnulus(position,
                                      r_in=inner_bkg, r_out=outer_bkg)
            normalization = aperture.area() / annulus.area()
    else:
        aperture = RectangularAperture(position, width, height, theta)
        # No background is computed for an extended source.

    for k in range(shape[0]):
        phot_table = aperture_photometry(data[k, :, :], aperture,
                                         method=method, subpixels=subpixels)
        net[k] = float(phot_table['aperture_sum'][0])
        if subtract_background:
            bkg_table = aperture_photometry(data[k, :, :], annulus,
                                            method=method, subpixels=subpixels)
            background[k] = float(bkg_table['aperture_sum'][0])
            net[k] = net[k] - background[k] * normalization

    # Check for NaNs in the wavelength array, flag them in the dq array,
    # and truncate the arrays if NaNs are found at endpoints (unless the
    # entire array is NaN).
    nan_mask = np.isnan(wavelength)
    n_nan = nan_mask.sum(dtype=np.intp)
    if n_nan > 0:
        log.warning("%d NaNs in wavelength array.", n_nan)
        dq[nan_mask] = np.bitwise_or(dq[nan_mask], dqflags.pixel['DO_NOT_USE'])
        not_nan = np.logical_not(nan_mask)
        flag = np.where(not_nan)
        if len(flag[0]) > 0:
            n_trimmed = flag[0][0] + nelem - (flag[0][-1] + 1)
            if n_trimmed > 0:
                log.info("Output arrays have been trimmed by %d elements",
                         n_trimmed)
                slc = slice(flag[0][0], flag[0][-1] + 1)
                wavelength = wavelength[slc]
                net = net[slc]
                background = background[slc]
                dq = dq[slc]
        else:
            dq |= dqflags.pixel['DO_NOT_USE']

    return (ra, dec, wavelength, net, background, dq)

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

#.........这里部分代码省略.........
					xpeaks_o[i]] - half_max_o)).argmin())
				nearest_above_y_e = ((np.abs(data_e[ypeaks_e[i]:-1,
					xpeaks_e[i]] - half_max_e)).argmin())
				nearest_below_y_e = ((np.abs(data_e[0:ypeaks_e[i],
					xpeaks_e[i]] - half_max_e)).argmin())
				fwhm.append((nearest_above_x_o + (xpeaks_o[i] -
					nearest_below_x_o)))
				fwhm.append((nearest_above_y_o + (ypeaks_o[i] -
					nearest_below_y_o)))
				fwhm.append((nearest_above_x_e + (xpeaks_e[i] -
					nearest_below_x_e)))
				fwhm.append((nearest_above_y_e + (ypeaks_e[i] -
					nearest_below_y_e)))
			
			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)
			

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import area [as 别名]
def addzb(fitsname, redo=False, fau_dir=None):
    telescopes = ['1','2','3','4']
#    night = (fitsname.split('/'))[3]
    night = (fitsname.split('/'))[-2]

    print 'beginning ' + fitsname


    # 2D spectrum
    h = pyfits.open(fitsname,mode='update')

    if 'BARYSRC4' in h[0].header.keys() and not redo:
        print fitsname + " already done"
        h.flush()
        h.close()
        return

    specstart = datetime.datetime.strptime(h[0].header['DATE-OBS'],"%Y-%m-%dT%H:%M:%S.%f")
    specmid = specstart + datetime.timedelta(seconds = h[0].header['EXPTIME']/2.0)
    specend = specstart + datetime.timedelta(seconds = h[0].header['EXPTIME'])

    t0 = datetime.datetime(2000,1,1)
    t0jd = 2451544.5

    aperture_radius = 3.398 # fiber radius in pixels
    annulus_inner = 2.0*aperture_radius
    annulus_outer = 3.0*aperture_radius

    for telescope in telescopes:
        print 'beginning telescope ' + telescope + ' on ' + fitsname

        # get the barycentric redshift for each time
        ra = h[0].header['TARGRA' + telescope]
        dec = h[0].header['TARGDEC' + telescope]
        try: pmra = h[0].header['PMRA' + telescope]
        except: pmra = 0.0
        try: pmdec = h[0].header['PMDEC' + telescope]
        except: pmdec = 0.0
        try: parallax = h[0].header['PARLAX' + telescope]
        except: parallax = 0.0
        try: rv = h[0].header['RV' + telescope]
        except: rv = 0.0

        objname = h[0].header['OBJECT' + telescope]
        if fau_dir is None:
            faupath = '/Data/t' + telescope + '/' + night + '/' + night + '.T' + telescope + '.FAU.' + objname + '.????.fits'
        else:
            faupath = fau_dir + '/t' + telescope + '/' + night + '/' + night + '.T' + telescope + '.FAU.' + objname + '.????.fits'
        guideimages = glob.glob(faupath)

#        if telescope == '2' and "HD62613" in fitsname: ipdb.set_trace()

        times = []
        fluxes = np.array([])

        for guideimage in guideimages:
            try:
                fauimage = pyfits.open(guideimage)
            except:
                print "corrupt file for " + guideimage
                continue

            # midtime of the guide image (UTC)
            midtime = datetime.datetime.strptime(fauimage[0].header['DATE-OBS'],"%Y-%m-%dT%H:%M:%S") +\
                datetime.timedelta(seconds=fauimage[0].header['EXPTIME']/2.0)
            
            # convert to Julian date
            midjd = t0jd + (midtime-t0).total_seconds()/86400.0

            # only look at images during the spectrum
            if midtime < specstart or midtime > specend: continue

            # find the fiber position
            try:
                fiber_x = fauimage[0].header['XFIBER' + telescope]
                fiber_y = fauimage[0].header['YFIBER' + telescope]
            except:
                print "keywords missing for " + guideimage
                continue

            # do aperture photometry
            positions = [(fiber_x,fiber_y)]
            apertures = CircularAperture(positions,r=aperture_radius)
            annulus_apertures = CircularAnnulus(positions, r_in=annulus_inner, r_out=annulus_outer)

            # calculate the background-subtracted flux at the fiber position
            rawflux_table = aperture_photometry(fauimage[0].data, apertures)
            bkgflux_table = aperture_photometry(fauimage[0].data, annulus_apertures)
            bkg_mean = bkgflux_table['aperture_sum'].sum() / annulus_apertures.area()
            bkg_sum = bkg_mean * apertures.area()
            flux = rawflux_table['aperture_sum'].sum() - bkg_sum
            
            # append to the time and flux arrays
            times.append(midjd)
            fluxes = np.append(fluxes,flux)

        if len(times) == 0:
            print "No guider images for " + fitsname + " on Telescope " + telescope + "; assuming mid time"
            if 'daytimeSky' in fitsname: 
                h[0].header['BARYCOR' + telescope] = ('UNKNOWN','Barycentric redshift')
#.........这里部分代码省略.........

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import area [as 别名]
    error in photometry is assumed as Poissonian -- sqrt(N)
    but smooth background error not included for now
    """
for eachFile in imFiles:
    print('Prcoessing %s ...' %(eachFile))
    
    try: 
        hdu = pyfits.open(eachFile)
        im = hdu[0].data
    except: print('File could not be opened: %s' %(eachFile))

    # create aperture and annulus objects
    apertures = CircularAperture(coords, r=aper_size)
    annulus_apertures = CircularAnnulus(coords, r_in=annulus, r_out=dannulus)
    
    npix_src, npix_bkg = apertures.area(),  annulus_apertures.area()
    
    # calculate the object and annulus flux
    data_error = np.sqrt(im)
    rawflux_table = aperture_photometry(im, apertures, error=data_error)
    bkgflux_table = aperture_photometry(im, annulus_apertures, error=data_error)
    # stack the two tables into one
    phot_table = hstack([rawflux_table, bkgflux_table], table_names=['raw', 'bkg'])

    # calculate bkg mean & normalize by area
    bkg_mean = phot_table['aperture_sum_bkg'] / npix_bkg

    # final photometric counts -- in ADUs
    Nsrc = phot_table['aperture_sum_raw'] - bkg_mean * npix_src
    # change counts to flux
    flux = gain * Nsrc / hdu[0].header['EXPTIME']

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

#.........这里部分代码省略.........
    for i in range(len(fnlist)):
        image = FPImage(fnlist[i])
        skyavg[i], skysig[i], _skyvar = image.skybackground()
        image.close()

    # Identify the stars in each image
    xlists = []
    ylists = []
    print "Identifying stars in each image..."
    for i in range(len(fnlist)):
        xlists.append([])
        ylists.append([])
        image = FPImage(fnlist[i])
        axcen = image.axcen
        aycen = image.aycen
        arad = image.arad
        sources = daofind(image.inty-skyavg[i],
                          fwhm=fwhm[i],
                          threshold=thresh*skysig[i]).as_array()
        for j in range(len(sources)):
            # If the source is not near the center or edge
            centermask = ((sources[j][1]-axcen)**2 +
                          (sources[j][2]-aycen)**2 > (0.05*arad)**2)
            edgemask = ((sources[j][1]-axcen)**2 +
                        (sources[j][2]-aycen)**2 < (0.95*arad)**2)
            if np.logical_and(centermask, edgemask):
                xlists[i].append(sources[j][1])
                ylists[i].append(sources[j][2])
        image.close()

    # Match objects between fields
    print "Matching objects between images..."
    xcoo = []
    ycoo = []
    for i in range(len(xlists[0])):
        # For each object in the first image
        accept = True
        for j in range(1, len(fnlist)):
            # For each other image
            dist2 = ((np.array(xlists[j])-xlists[0][i])**2 +
                     (np.array(ylists[j])-ylists[0][i])**2)
            if (min(dist2) > tolerance**2):
                accept = False
                break
        if accept:
            # We found an object at that position in every image
            xcoo.append(xlists[0][i])
            ycoo.append(ylists[0][i])

    # Create coordinate arrays for the photometry and shifting
    x = np.zeros((len(fnlist), len(xcoo)))
    y = np.zeros_like(x)
    for i in range(len(xcoo)):
        # For every object found in the first image
        for j in range(len(fnlist)):
            # Find that object in every image
            dist2 = ((np.array(xlists[j])-xcoo[i])**2 +
                     (np.array(ylists[j])-ycoo[i])**2)
            index = np.argmin(dist2)
            x[j, i] = xlists[j][index]
            y[j, i] = ylists[j][index]

    # Do aperture photometry on the matched objects
    print "Performing photometry on matched stars..."
    counts = np.zeros_like(x)
    dcounts = np.zeros_like(x)
    for i in range(len(fnlist)):
        image = FPImage(fnlist[i])
        apertures = CircularAperture((x[i], y[i]), r=2*fwhm[i])
        annuli = CircularAnnulus((x[i], y[i]), r_in=3*fwhm[i], r_out=4*fwhm[i])
        phot_table = aperture_photometry(image.inty,
                                         apertures, error=np.sqrt(image.vari))
        sky_phot_table = aperture_photometry(image.inty, annuli,
                                             error=np.sqrt(image.vari))
        counts[i] = phot_table["aperture_sum"] / apertures.area()
        counts[i] -= sky_phot_table["aperture_sum"] / annuli.area()
        counts[i] *= apertures.area()
        dcounts[i] = phot_table["aperture_sum_err"] / apertures.area()
        image.close()

    # Calculate the shifts and normalizations
    norm, dnorm = calc_norm(counts, dcounts)
    for i in range(x.shape[1]):
        x[:, i] = -(x[:, i] - x[0, i])
        y[:, i] = -(y[:, i] - y[0, i])
    xshifts = np.average(x, axis=1)
    yshifts = np.average(y, axis=1)

    # Normalize the images and put shifts in the image headers
    for i in range(len(fnlist)):
        image = FPImage(fnlist[i], update=True)
        image.phottog = "True"
        image.dnorm = dnorm[i]
        image.inty /= norm[i]
        image.vari = image.vari/norm[i]**2
        image.xshift = xshifts[i]
        image.yshift = yshifts[i]
        image.close()

    return

# 需要导入模块: from photutils import CircularAperture [as 别名]
# 或者: from photutils.CircularAperture import area [as 别名]
        # overplotted apertures, smaller should be better.
        circ_apertures = CircularAperture(positions, r=arb.radius[band])
        annulus_apertures = CircularAnnulus(
                positions, r_in=arb.annuli_r[band][0], r_out=arb.annuli_r[band][1])
        apers = [circ_apertures, annulus_apertures]
        phot_table = aperture_photometry(image, apers)

        #################################
        # Local background subtraction. #
        #################################

        ## Method 1: Follow
        ## http://photutils.readthedocs.io/en/stable/photutils (...)
        ##       /aperture.html#local-background-subtraction
        bkg_mean = phot_table['aperture_sum_1'] / annulus_apertures.area()
        bkg_sum = bkg_mean * circ_apertures.area()
        final_sum = phot_table['aperture_sum_0'] - bkg_sum
        phot_table['residual_aperture_sum'] = final_sum

        ## Method 2: Follow
        ## https://github.com/astropy/photutils/pull/453, sigclipping away stars
        ## in the annulus.
        #ann_masks = annulus_apertures.to_mask(method='center')
        #ann_masked_data = [am.apply(image) for am in ann_masks]

        ## Sigma clip stars in annular aperture.
        #pre_sc_median = [np.nanmedian(amd[amd != 0])
        #                    for amd in ann_masked_data]

        #pre_sc_std = [
        #        (np.nanmedian(np.abs(amd[amd != 0] - pre_sc_median[ix])))*1.483

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

#.........这里部分代码省略.........
    annulus_sum_err = np.array(annulus_sum_err)

    # construct metadata for output table
    meta = OrderedDict()
    meta['instrument'] = datamodel.meta.instrument.name
    meta['detector'] = datamodel.meta.instrument.detector
    meta['channel'] = datamodel.meta.instrument.channel
    meta['subarray'] = datamodel.meta.subarray.name
    meta['filter'] = datamodel.meta.instrument.filter
    meta['pupil'] = datamodel.meta.instrument.pupil
    meta['target_name'] = datamodel.meta.target.catalog_name
    meta['xcenter'] = xcenter
    meta['ycenter'] = ycenter
    ra_icrs, dec_icrs = datamodel.meta.wcs(xcenter, ycenter)
    meta['ra_icrs'] = ra_icrs
    meta['dec_icrs'] = dec_icrs
    meta['apertures'] = info

    # initialize the output table
    tbl = QTable(meta=meta)

    if hasattr(datamodel, 'int_times') and datamodel.int_times is not None:
        nrows = len(datamodel.int_times)
    else:
        nrows = 0
    if nrows == 0:
        log.warning("There is no INT_TIMES table in the input file.")

    if nrows > 0:
        shape = datamodel.data.shape
        if len(shape) == 2:
            num_integ = 1
        else:                                   # len(shape) == 3
            num_integ = shape[0]
        int_start = datamodel.meta.exposure.integration_start
        if int_start is None:
            int_start = 1
            log.warning("INTSTART not found; assuming a value of %d",
                        int_start)

        # Columns of integration numbers & times of integration from the
        # INT_TIMES table.
        int_num = datamodel.int_times['integration_number']
        mid_utc = datamodel.int_times['int_mid_MJD_UTC']
        offset = int_start - int_num[0]                 # both are one-indexed
        if offset < 0:
            log.warning("Range of integration numbers in science data extends "
                        "outside the range in INT_TIMES table.")
            log.warning("Can't use INT_TIMES table.")
            del int_num, mid_utc
            nrows = 0                   # flag as bad
        else:
            log.debug("Times are from the INT_TIMES table.")
            time_arr = mid_utc[offset: offset + num_integ]
            int_times = Time(time_arr, format='mjd', scale='utc')
    else:
        log.debug("Times were computed from EXPSTART and TGROUP.")

        dt = (datamodel.meta.exposure.group_time *
              (datamodel.meta.exposure.ngroups + 1))
        dt_arr = (np.arange(1, 1 + datamodel.meta.exposure.nints) *
                  dt - (dt / 2.))
        int_dt = TimeDelta(dt_arr, format='sec')
        int_times = (Time(datamodel.meta.exposure.start_time, format='mjd') +
                     int_dt)

    tbl['MJD'] = int_times.mjd

    tbl['aperture_sum'] = aperture_sum
    tbl['aperture_sum_err'] = aperture_sum_err

    if not sub64p_wlp8:
        tbl['annulus_sum'] = annulus_sum
        tbl['annulus_sum_err'] = annulus_sum_err

        annulus_mean = annulus_sum / bkg_aper.area()
        annulus_mean_err = annulus_sum_err / bkg_aper.area()
        tbl['annulus_mean'] = annulus_mean
        tbl['annulus_mean_err'] = annulus_mean_err

        aperture_bkg = annulus_mean * phot_aper.area()
        aperture_bkg_err = annulus_mean_err * phot_aper.area()
        tbl['aperture_bkg'] = aperture_bkg
        tbl['aperture_bkg_err'] = aperture_bkg_err

        net_aperture_sum = aperture_sum - aperture_bkg
        net_aperture_sum_err = np.sqrt(aperture_sum_err ** 2 +
                                       aperture_bkg_err ** 2)
        tbl['net_aperture_sum'] = net_aperture_sum
        tbl['net_aperture_sum_err'] = net_aperture_sum_err
    else:
        colnames = ['annulus_sum', 'annulus_sum_err', 'annulus_mean',
                    'annulus_mean_err', 'aperture_bkg', 'aperture_bkg_err']
        for col in colnames:
            tbl[col] = np.full(nimg, np.nan)

        tbl['net_aperture_sum'] = aperture_sum
        tbl['net_aperture_sum_err'] = aperture_sum_err

    return tbl