Python pyregion.open函数代码示例

def test_region():

    ref_region_name = "test01_img.reg"

    region_list = ["test01_fk5_sexagecimal.reg",
                   "test01_gal.reg",
                   "test01_ds9_physical.reg",
                   "test01_fk5_degree.reg",
                   "test01_mixed.reg",
                   "test01_ciao.reg",
                   "test01_ciao_physical.reg",
                   ]

    header = demo_header()

    ref_region = pyregion.open(join(rootdir,ref_region_name)).as_imagecoord(header)

    for reg_name in region_list:
        r = pyregion.open(join(rootdir,reg_name)).as_imagecoord(header)
        for reg0, reg in zip(ref_region, r):
            if reg.name == "rotbox":
                reg.name = "box"

            assert reg0.name == reg.name
            if reg0.name in ["ellipse", "box"]:
                assert np.allclose(reg0.coord_list[:-1], reg.coord_list[:-1],
                                   atol=0.01)
                a0 = reg0.coord_list[-1]
                a1 = fix_lon(reg.coord_list[-1], 0)
                assert np.allclose([a0], [a1], atol=0.02)
            else:
                assert np.allclose(reg0.coord_list, reg.coord_list,
                                   atol=0.01)
            assert reg0.exclude == reg.exclude

def subcubes_from_ds9(cube, region_file='../nro_maps/SouthShells.reg', pad_factor=1., shape='exact'):
    """
    Extracts subcubes using the ds9 region file.
    
    Parameters
    ----------
    cube : SpectralCube, str
        The cube to be chopped. Must be type spectral_cube.SpectralCube or str filename.
    region_file : str
        Path to a ds9 region file.
    pad_factor : float, optional
        Expand the subcube around the region by this factor.
    shape : {'square', 'exact'}
        The shape of the subcube returned. 'square' returns the
        smallest square subcube that contains the region.
        'exact' returns only the pixels contained within the region.
    
    Returns
    -------
    subcubes: list of SpectralCube of SpectralCube
    """
    from spectral_cube import SpectralCube
    import pyregion

    try:
        #If cube is a str filename, read a SpectralCube.
        cube = SpectralCube.read(cube)
    except ValueError:
        pass

    if shape == 'square':
        import astropy.units as u
        subcube_list = []
        region_list = pyregion.open(region_file)
        for region in region_list:
            half_width = region.coord_list[2] * pad_factor * u.deg
            ra_center = region.coord_list[0] * u.deg
            dec_center = region.coord_list[1] * u.deg
            ra_range = [ra_center - half_width, ra_center + half_width]
            dec_range = [dec_center - half_width, dec_center + half_width]
            #print(ra_range, dec_range)
            subcube_list.append(cube.subcube(ra_range[1], ra_range[0], dec_range[0], dec_range[1]))
    if shape == 'exact':
        region_list = pyregion.open(region_file)
        subcube_list = []
        for region in region_list:
            
            if pad_factor != 1.:
                new_string = '{};{}({},{},{}")'.format(region.coord_format, region.name,
                                region.coord_list[0], region.coord_list[1],
                                region.coord_list[2]*3600.*pad_factor)
                region = pyregion.parse(new_string)[0]
                
            subcube_list.append(cube.subcube_from_ds9region(pyregion.ShapeList([region])))
    if len(subcube_list) == 1:
        return subcube_list[0]
    else:
        return subcube_list

    def __init__(self, x, y, regfile, header=None):
        
	""" Parse the ds9 region file here. Also save the number 
	of regions and the x y dimensions of the data"""
	if(header):
            self.reg = pyreg.open(regfile).as_imagecoord(header=header)
        else:
            self.reg = pyreg.open(regfile)

	self.filter = self.reg.get_filter()
	self.x = x
	self.y = y
	self.nreg = len(self.reg)

def radioflux(files,fgr,bgr=None,individual=False,action='Flux',fluxerr=0,nsigma=0,verbose=False):
    """Determine the flux in a region file for a set of files. This is the
    default action for the code called on the command line, but
    may be useful to other code as well.

    Keyword arguments:
    files -- list of files (mandatory)
    fdr -- foreground region name (mandatory)
    bgr -- background region name (optional)
    individual -- separate region into individual sub-regions
    action -- what to do once fluxes are measured: allows a user-defined action
              which must be a drop-in replacement for printflux
    fluxerr -- flux error in % for spidxmap
    nsigma -- keep only pixels above these sigma level in ALL maps (bgr must be specified)
    """
    action = {'flux':printflux, 'mean':printmean, 'spidx':printspidx}[action]

    rms = [] # radio maps
    for filename in files:
        rms.append(radiomap(filename,verbose=verbose))

    # if using the sigma all the images must have the same size
    if nsigma > 0: assert all(rms[i].d[0].size == rms[0].d[0].size for i in range(len(rms)))
    # initial mask
    mask = (np.zeros_like(rms[0].d) == 0)

    bgs = [] #1d list: [ radiomap ]
    for rm in rms:
        if bgr:
            bg_ir=pyregion.open(bgr).as_imagecoord(rm.headers[0])
            bg=applyregion(rm,bg_ir)
            bgs.append(bg.rms)
            # likely brakes with channelled images
            if nsigma > 0: mask = np.logical_and(mask, np.array(rm.d) > (np.array(bg.rms)*nsigma) )
        else:
            bgs.append(None)

    fgs = [] # 2d list: [ radiomap x forground_region]
    for i, rm in enumerate(rms):
        fg_ir=pyregion.open(fgr).as_imagecoord(rm.headers[0])
        if individual:
            fgs.append([])
            for fg_ir_split in fg_ir:
                fg=pyregion.ShapeList([fg_ir_split])
                fgs[-1].append(applyregion(rm,fg,offsource=bgs[i],mask=mask))
        else:
            fgs.append([applyregion(rm,fg_ir,offsource=bgs[i],mask=mask)])

    # cycle before on regions and than on rm
    fgs = np.array(fgs).swapaxes(0,1)
    action(fgs, fluxerr)

def photometry(dither=1):
    # read in region file
    if dither == 2:
        image,hdr = pyfits.getdata('L1544_J_dither2.fits', header=True)
        regpix = pyregion.open('L1544_pix_dither2.reg')
    else:
        image,hdr = pyfits.getdata('L1544_J_dither1.fits', header=True)
        regpix = pyregion.open('L1544_pix_dither1.reg')
        nx1 = hdr['NAXIS1']
        nx2 = hdr['NAXIS2']
#        calc_maglim(image, nx1, nx2)

    regsky = pyregion.open('L1544_sky.reg')
    
    names = ('star', 'ra', 'dec', 'xc', 'yc', 'flux', 'eflux', 
             'mag', 'emag')
    dtypes = ('int', 'S11', 'S11', 'float', 'float', 'float', 'float',
              'float', 'float')
    t = Table(data=None, names=names, dtype=dtypes)

    # centroid sources
    cenx = np.zeros(len(regpix), dtype=float)
    ceny = np.zeros(len(regpix), dtype=float)
    for i,star in enumerate(regpix):
        id = star.comment
        ra = regsky[i].coord_list[0]
        dec = regsky[i].coord_list[1]

#        x_init = star.coord_list[0]
#        y_init = star.coord_list[1]
        cenx = star.coord_list[0]
        ceny = star.coord_list[1]
#        # centroid by fitting a 2D Gaussian
#        cenx,ceny = gauss_centroid(image, x_init, y_init)
        
        t.add_row([id, ra, dec, cenx, ceny, 0.0, 0.0, 0.0, 0.0])
    
    # perform aperture photometry
    # 30 pixels in radius and sky annuli at 30-50 pixels
    table = apphot(image, t, 15, [20,30], return_mag=True)

    formats = {'xc':'%8.3f', 'xc':'%8.3f', 'flux':'%7.2f', \
               'eflux':'%6.2f', 'mag':'%6.3f', 'emag':'%5.3f'}

    if dither == 2:
        table.write('flux_dither2.txt', format='ascii.tab', formats=formats)
    else:
        table.write('flux_dither1.txt', format='ascii.tab', formats=formats)

def read_mask(maskfile,fitsfile,wolly=False,return_shape=False):
    '''Parse mask file and return masked region. If 'wolly' is True, return two masks.'''
    import pyregion
    if isinstance(fitsfile,str):
        hdu = fits.open(fitsfile)
    else:
        hdu = fitsfile
    reg = pyregion.open(maskfile)

    # total mask
    mask = ~pyregion.get_mask(reg,hdu[0])
    
    if wolly:
        # return three masks: (top, bottom, total)
        if len(reg) != 2:
            raise AttributeError("Wolly masks must have two regions")

        # make total mask
        reg = [pyregion.ShapeList([r]) for r in reg]
        pmask = [~pyregion.get_mask(r,hdu[0]) for r in reg]
        pmask.append(mask)
        mask = pmask

        if return_shape:
            shapes = [mask_shape(m) for m in mask]
            mask = tuple(pmask + shapes)
            #mask is (top,bottom,total,shape_top,shape_bot,shape_total)
    else:
        if return_shape:
            mask = (mask,mask_shape(mask))
    return mask

def getregionboxcenter(regionfile):
    """
    Extract box center of a DS9 box region. 
    Input is regionfile Return NDPPP compatible string for phasecenter shifting
    """
    r = pyregion.open(regionfile)
    
    if len(r[:]) > 1:
      print 'Only one region can be specified, your file contains', len(r[:])
      sys.exit() 
    
    if r[0].name != 'box':
      print 'Only box region supported'
      sys.exit()
    
    ra  = r[0].coord_list[0]
    dec = r[0].coord_list[1]
    boxsizex = r[0].coord_list[2]
    boxsizey = r[0].coord_list[3]
    angle = r[0].coord_list[4]
    if boxsizex != boxsizey:
      print 'Only a square box region supported, you have these sizes:', boxsizex, boxsizey
      sys.exit()
    if np.abs(angle) > 1:
      print 'Only normally oriented sqaure boxes are supported, your region is oriented under angle:', angle
      sys.exit()   
    
    regioncenter =  ('{:12.8f}'.format(ra) + 'deg,' + '{:12.8f}'.format(dec) + 'deg').replace(' ', '')
    return regioncenter

def ds9(regionfile, header, zorder=3, **kwargs):
    """
    Wrapper to return a PatchCollection given a ds9 region file
    and a fits header.

    zorder - defaults to 3 so that regions are on top of contours
    """

    # read region file
    rr = pyregion.open(regionfile)

    # convert coordinates to image coordinates
    rrim = rr.as_imagecoord(header)

    # pyregion and aplpy both correct for the FITS standard origin=1,1
    # need to avoid double-correcting
    for r in rrim:
        r.coord_list[0] += 1
        r.coord_list[1] += 1

    if 'text_offset' in kwargs:
        text_offset = kwargs['text_offset']
        del kwargs['text_offset']
    else:
        text_offset = 5.0

    # grab the shapes to overplot
    pp, aa = rrim.get_mpl_patches_texts(text_offset=text_offset)

    PC = ArtistCollection(pp, **kwargs) # preserves line style (dashed)
    TC = ArtistCollection(aa, **kwargs)
    PC.set_zorder(zorder)
    TC.set_zorder(zorder)

    return PC, TC

def _cut_stamps_regions(source_file, region_file, source_ext=0, file_prefix='',
                        selected_labels=None, **kwargs):
    """
    Cuts stamps from a FITS image based on a DS9 region file. See docstring for
    cut_stamps for argument explanations
    """
    # Open region file first, so parse errors will fail before trying to open
    # a potentially large FITS file
    regions = pyregion.open(region_file)

    # Open cutting image
    hdu_list = fits.open(source_file, mode='readonly', ignore_missing_end=True)

    # Translate to image coordinates, since we're working in pixel space
    regions = regions.as_imagecoord(hdu_list[source_ext].header)

    # If a region has a text label, use that as the default file name
    # reg.attr[1] is a dictionary of extended region attributes
    names = [reg.attr[1].get('text', str(num))
             for num, reg in enumerate(regions)]
    if selected_labels is None:
        selected_labels = set(names)

    all_stamp_names = []
    for name in set(names) & selected_labels:
        regs_for_name = [reg for reg, nm in zip(regions, names) if nm == name]
        stamp = _stamp_from_regions(hdu_list, regs_for_name,
                                    extension=source_ext, **kwargs)
        filename = file_prefix + name + '.fits'
        stamp.writeto(filename, clobber=True)
        all_stamp_names += [filename]
    hdu_list.close()

    return all_stamp_names

def gal2cel(regfile):
    """
    Converts a region file from galactic to celestial coordinates including
    position angle reference from the center of the box (right now only works
    on box regions)

    Requires pyregion with the ShapeList.write() function implemented...
    not clear if that exists in 1.0
    """
    reg = pyregion.open(regfile)

    for R in reg:
        if R.name == "box":
            x, y, dx, dy, angle = R.coord_list

            # posn = coords.Position([x,y],system='galactic')
            # ra,dec = posn.j2000()
            posn = coordinates.Galactic(x * u.deg, y * u.deg)
            ra, dec = posn.fk5.ra.deg, posn.fk5.dec.deg

            newang = posang.posang(x - dx, y, x + dx, y, system="galactic")

            coord_list = [ra, dec, dx, dy, angle - newang - 90]

            R.coord_format = "fk5"
            R.coord_list = coord_list
            R.params = coord_list

    reg.write(regfile[:-4] + "_fk5.reg")

def make_reg_masks(regfile,shape):
    r = pyregion.open(regfile)

    if len(r) != 2:
        raise Exception('Exactly two box regions required')

    paths = []
    for reg in r:
        paths.append(get_region_box(reg.coord_list))

    #Always have A be the top half
    if paths[0].get_extents().ymax > paths[1].get_extents().ymax:
        pathA = paths[0]
        pathB = paths[1]
    else:
        pathA = paths[1]
        pathB = paths[2]


    print 'Building skymasks'
    maskA = np.array([True if pathA.contains_point([x,y]) else False for x,y in np.ndindex(shape)])
    maskA = maskA.reshape(shape).T
    
    maskB = np.array([True if pathB.contains_point([x,y]) else False for x,y in np.ndindex(shape)])
    maskB = maskB.reshape(shape).T

    return (~maskA, ~maskB)

def show_region(fig, region_list):
    h = demo_header()

    n = len(region_list)
    nx = int(math.ceil(n ** .5))
    ny = int(math.ceil(1. * n / nx))

    nrows_ncols = (ny, nx)

    grid = [plt.subplot(ny, nx, i + 1) for i in range(n)]

    for ax, reg_name in zip(grid, region_list):
        ax.set_aspect(1)

        r = pyregion.open(reg_name).as_imagecoord(h)

        patch_list, text_list = r.get_mpl_patches_texts()
        for p in patch_list:
            ax.add_patch(p)
        for t in text_list:
            ax.add_artist(t)

        if plt.rcParams["text.usetex"]:
            reg_name = reg_name.replace("_", r"\_")
        ax.set_title(reg_name, size=10)
        for t in ax.get_xticklabels() + ax.get_yticklabels():
            t.set_visible(False)

    return grid

def mask_satellite_trail(flc='jcdu36e5q_flc.backup.fits', mask_files=['jcdu36e5q_mask_1.reg', 'jcdu36e5q_mask_2.reg']):
    """
    Use Pyregion to mask satellite trails with a region file. 
    
    Note: save the region file in "fk5" sky coordinates
    """
    import astropy.io.fits as pyfits
    import pyregion
    
    ### open the FITS image for modifying
    im = pyfits.open(flc) 
    
    ### Make the mask for both "DQ" extensions for both ACS/WFC chips
    for extension in [1,2]:
        ### Get a mask returning 1 *within* the region and 0 elsewhere
        im['SCI',extension].header.remove('CPDIS1')
        im['SCI',extension].header.remove('CPDIS2')
        reg = pyregion.open(mask_files[extension-1]).as_imagecoord(im['SCI', extension].header)
        mask = reg.get_mask(im['SCI', extension])
        ### Apply the mask to the DQ extension
        im['DQ', extension].data |= 4096*mask
    
    imf = pyfits.open(flc, mode='update')
    for extension in [1,2]:
        imf['DQ', extension].data = im['DQ',extension].data
        
    ### Write the file
    imf.flush()

def area_within_reg(reg_name,base_fits=m31):
    '''
    returns deg^s 
    deg^2 to arcsec^2: areas*3600.*3600.
    deg^2 to arcmin^2: areas*3600.
    '''
    # finds area of concentric reg shapes
    # use GenUtils.area for the area of a polygon.
    tot_areas,areas = [],[]
    proj = wcs.Projection(f_header)
    r = pyregion.open(reg_name).as_imagecoord(header=f_header)
    patch_list, artist_list = r.get_mpl_patches_texts()
    for p in patch_list[::-1]:
        verts_xy = p.get_xy()
        pixels = (verts_xy[:,0],verts_xy[:,1])
        ra_v,dec_v = proj.toworld(pixels)
        verts = np.column_stack((ra_v,dec_v))
        tot_areas.append(area(verts))
        
    for i in range(len(tot_areas)):
        if i == len(tot_areas)-1:
            areas.append(tot_areas[i])
        else:
            areas.append(tot_areas[i]-tot_areas[i+1])

    # deg^2 to arcsec^2: areas*3600.*3600.
    # deg^2 to arcmin^2: areas*3600.
    return np.array(areas)

def regtoapo(regfile,outfile,filtercolor=None):
    """
    Take an input ds9 .reg file and create an output file
    formatted to work with APO's TUI

    filtercolor - if specified, will ignore any regions of this color
    """

    reg = pyregion.open(regfile)
    outf = open(outfile,'w')

    for r in reg:
        if r.attr[1]['color'] == filtercolor or not r.attr[1].has_key('text'):
            continue
        if r.name == 'box':
            x,y,dx,dy,posang = r.coord_list
            radec = position_region(r).hmsdms()
            # ds9's position angles are 90 degrees offset from APO's
            if posang+90 > 360: posang -= 360
            print >>outf,'%32s %26s rotangle=%f' % ( '"%s"' % r.attr[1]['text'],radec,posang+90)
        else:
            radec = position_region(r).hmsdms()
            print >>outf,'%32s %26s' % ('"%s"' % r.attr[1]['text'],radec)

    outf.close()