Python SpectralCube.with_mask方法代码示例

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
    def __init__(self, cube, wcs=None, mask=None, sigma=None, empty_channel=0,
                 keep_threshold_mask=True, distance=None, galaxy_props={}):
        super(BubbleFinder, self).__init__()

        if not isinstance(cube, SpectralCube):
            if wcs is None:
                raise TypeError("When cube is not a SpectralCube, wcs must be"
                                " given.")
            cube = SpectralCube(cube, wcs)
            if mask is not None:
                cube = cube.with_mask(mask)

        self.cube = cube

        self.empty_channel = empty_channel

        if sigma is None:
            self.estimate_sigma()
        else:
            self.sigma = sigma

        self.keep_threshold_mask = keep_threshold_mask
        self._mask = None
        self.distance = distance
        self.galaxy_props = galaxy_props

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
def load_and_reduce(filename, add_noise=False, rms_noise=0.001,
                    nsig=3):
    '''
    Load the cube in and derive the property arrays.
    '''

    if add_noise:
        if rms_noise is None:
            raise TypeError("Must specify value of rms noise.")

        cube, hdr = getdata(filename, header=True)

        from scipy.stats import norm
        cube += norm.rvs(0.0, rms_noise, cube.shape)

        sc = SpectralCube(data=cube, wcs=WCS(hdr))

        mask = LazyMask(np.isfinite, sc)
        sc = sc.with_mask(mask)

    else:
        sc = filename

    reduc = Mask_and_Moments(sc, scale=rms_noise)
    reduc.make_mask(mask=reduc.cube > nsig * reduc.scale)
    reduc.make_moments()
    reduc.make_moment_errors()

    return reduc.to_dict()

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
def reduce_and_save(filename, add_noise=False, rms_noise=0.001,
                    output_path="", cube_output=None,
                    nsig=3, slicewise_noise=True):
    '''
    Load the cube in and derive the property arrays.
    '''

    if add_noise:
        if rms_noise is None:
            raise TypeError("Must specify value of rms noise.")

        cube, hdr = getdata(filename, header=True)

        # Optionally scale noise by 1/10th of the 98th percentile in the cube
        if rms_noise == 'scaled':
            rms_noise = 0.1*np.percentile(cube[np.isfinite(cube)], 98)

        from scipy.stats import norm
        if not slicewise_noise:
            cube += norm.rvs(0.0, rms_noise, cube.shape)
        else:
            spec_shape = cube.shape[0]
            slice_shape = cube.shape[1:]
            for i in range(spec_shape):
                cube[i, :, :] += norm.rvs(0.0, rms_noise, slice_shape)

        sc = SpectralCube(data=cube, wcs=WCS(hdr))

        mask = LazyMask(np.isfinite, sc)
        sc = sc.with_mask(mask)

    else:
        sc = filename

    reduc = Mask_and_Moments(sc, scale=rms_noise)
    reduc.make_mask(mask=reduc.cube > nsig * reduc.scale)

    reduc.make_moments()
    reduc.make_moment_errors()

    # Remove .fits from filename
    save_name = filename.split("/")[-1][:-4]

    reduc.to_fits(output_path+save_name)

    # Save the noisy cube too
    if add_noise:
        if cube_output is None:
            reduc.cube.hdu.writeto(output_path+save_name)
        else:
            reduc.cube.hdu.writeto(cube_output+save_name)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
def warp_ellipse_to_circle(cube, a, b, pa, stop_if_huge=True):
    '''
    Warp a SpectralCube such that the given ellipse is a circle int the
    warped frame.

    Since you should **NOT** be doing this with a large cube, we're going
    to assume that the given cube is a subcube centered in the middle of the
    cube.

    This requires a rotation, then scaling. The equivalent matrix is:
    [b cos PA    b sin PA]
    [-a sin PA   a cos PA ].

    '''

    if cube._is_huge:
        if stop_if_huge:
            raise Warning("The cube has the huge flag enabled. Disable "
                          "'stop_if_huge' if you would like to continue "
                          "anyways with the warp.")
        else:
            warn("The cube has the huge flag enabled. This may use a lot "
                 "of memory!")

    # Let NaNs be 0
    data = cube.with_fill_value(0.0).filled_data[:].value

    warped_array = []

    for i in range(cube.shape[0]):
        warped_array.append(nd.zoom(nd.rotate(data[i], np.rad2deg(-pa)),
                                    (1, a / b)))

    warped_array = np.array(warped_array)

    # We want to mask outside of the original bounds
    mask = np.ones(data.shape[1:])
    warp_mask = \
        np.isclose(nd.zoom(nd.rotate(mask, np.rad2deg(-pa)),
                           (1, a / b)), 1)

    # There's probably a clever way to transform the WCS, but all the
    # solutions appear to need pyast/starlink. The output of the wrap should
    # give a radius of b and the spectral dimension is unaffected.
    # Also this is hidden and users won't be able to use this weird cube
    # directly
    warped_cube = SpectralCube(warped_array * cube.unit, cube.wcs)
    warped_cube = warped_cube.with_mask(warp_mask)

    return warped_cube

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
def load_and_reduce(filename, add_noise=False, rms_noise=0.001,
                    nsig=3, slicewise_noise=True):
    '''
    Load the cube in and derive the property arrays.
    '''

    if add_noise:
        if rms_noise is None:
            raise TypeError("Must specify value of rms noise.")

        cube, hdr = getdata(filename, header=True)

        # Optionally scale noise by 1/10th of the 98th percentile in the cube
        if rms_noise == 'scaled':
            rms_noise = 0.1*np.percentile(cube[np.isfinite(cube)], 98)

        from scipy.stats import norm
        if not slicewise_noise:
            cube += norm.rvs(0.0, rms_noise, cube.shape)
        else:
            spec_shape = cube.shape[0]
            slice_shape = cube.shape[1:]
            for i in range(spec_shape):
                cube[i, :, :] += norm.rvs(0.0, rms_noise, slice_shape)

        sc = SpectralCube(data=cube, wcs=WCS(hdr))

        mask = LazyMask(np.isfinite, sc)
        sc = sc.with_mask(mask)

    else:
        sc = filename

    reduc = Mask_and_Moments(sc, scale=rms_noise)
    reduc.make_mask(mask=reduc.cube > nsig * reduc.scale)
    reduc.make_moments()
    reduc.make_moment_errors()

    return reduc.to_dict()

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
def calc_structure_fcn(catalog,
                       bootiter=0,doPlot=False):
    
    cat = catalog
    if 'sf_offset' not in cat.keys():
        keylist = ['sf_offset','sf_index','sf_ngood',
                   'sf_index_err','sf_offset_err', 'vca_index',
                   'vca_index_err']
        for thiskey in keylist:
            c = Column(np.zeros(len(cat))+np.nan,name=thiskey)
            cat.add_column(c)
        
    current_open_file = ''
    for cloud in cat:
        root = cloud['orig_file'].split('_')
        orig_file = 'COHRS_{0}_{1}.fits'.format(root[-2], root[-1])
        asgn_file = cloud['orig_file'] + '_fasgn.fits'
        if os.path.isfile(datadir+'COHRS_tiles/'+orig_file):
            if current_open_file != datadir+'COHRS_tiles/'+orig_file:
                hdu = fits.open(datadir+'COHRS_tiles/'+orig_file,memmap=False)
                cat.write('cohrs_structurefunc.fits',overwrite=True)
                w = wcs.WCS(hdu[0].header)
                hdr2 = w.to_header()
                hdr2['BMAJ'] = 15./3600
                hdr2['BMIN'] = 15./3600
                hdr2['BPA'] = 0.
                co = SpectralCube(hdu[0].data,w,header=hdr2)
                hdu = fits.open(datadir+'ASGN/'+asgn_file,memmap=False)
                w = wcs.WCS(hdu[0].header)
                hdr2 = w.to_header()
                hdr2['BMAJ'] = 15./3600
                hdr2['BMIN'] = 15./3600
                hdr2['BPA'] = 0.
                asgn = SpectralCube(hdu[0].data,w,header=hdr2)

#                masked_co = co.with_mask(asgn>0*u.dimensionless_unscaled)
#                moment = masked_co.moment(0)
                
                current_open_file = datadir+'COHRS_tiles/'+orig_file

            print(cloud['cloud_id'])
            mask = (asgn == cloud['cloud_id']*u.dimensionless_unscaled)
            subcube = co.with_mask(mask)
            subcube = subcube.subcube_from_mask(mask)
#            subcube.moment0().quicklook()
#            plt.show()
            if subcube.shape[0] > 10:
                # zeros = np.zeros_like(subcube) # np.random.randn(*subcube.shape)*0.5
                # concatdata = np.vstack([subcube.filled_data[:],zeros])
                # hdr2 = subcube.wcs.to_header()
                # hdr2['BMAJ'] = 15./3600
                # hdr2['BMIN'] = 15./3600
                # hdr2['BPA'] = 0.
                # newcube = SpectralCube(concatdata,subcube.wcs,header=hdr2)
                # if True:
                try:
                    pcaobj = pca.PCA(subcube,
                                     distance=cloud['bgps_distance_pc']*u.pc)
                    vcaobj = vca.VCA(subcube,
                                     distance=cloud['bgps_distance_pc']*u.pc)
                    vcaobj.run(verbose=False, high_cut = 0.3/u.pix)
                    pcaobj.run(verbose=True,mean_sub=True,
                               min_eigval=0.9,eigen_cut_method='proportion')
                    cloud['sf_index']= pcaobj.index
                    cloud['sf_offset'] = pcaobj.intercept.value
                    cloud['sf_ngood'] = np.min([np.isfinite(pcaobj.spatial_width).sum(),
                                                np.isfinite(pcaobj.spectral_width).sum()])
                    
                    cloud['sf_index_err'] = (pcaobj.index_error_range[1]-pcaobj.index_error_range[0])*0.5
                    cloud['sf_offset_err'] = ((pcaobj.intercept_error_range[1]-pcaobj.intercept_error_range[0])*0.5).value
                    cloud['vca_index'] = vcaobj.slope
                    cloud['vca_index_err'] = vcaobj.slope_err
                    print('{0} +/- {1}'.format(cloud['sf_index'],
                                               cloud['sf_index_err']),
                          cloud['sf_offset'])
                    # import pdb; pdb.set_trace()
                # else:
                except:
                    pass
        
    cat.write('output_catalog2.fits',overwrite=True)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
dataset1 = np.load(path1)

cube1 = np.empty((500, 32, 32))

count = 0
for posn, kept in zip(*dataset1["channels"]):
    posn = int(posn)
    if kept:
        cube1[posn, :, :] = dataset1["cube"][count, :, :]
        count += 1
    else:
        cube1[posn, :, :] = ra.normal(0.005, 0.005, (32, 32))

sc1 = SpectralCube(data=cube1, wcs=WCS(header))
mask = LazyMask(np.isfinite, sc1)
sc1 = sc1.with_mask(mask)
# Set the scale for the purposes of the tests
props1 = Moments(sc1, scale=0.003031065017916262 * u.Unit(""))
# props1.make_mask(mask=mask)
props1.make_moments()
props1.make_moment_errors()

dataset1 = props1.to_dict()

moment0_hdu1 = fits.PrimaryHDU(dataset1["moment0"][0],
                               header=dataset1["moment0"][1])

moment0_proj = Projection.from_hdu(moment0_hdu1)

##############################################################################

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
#    88.63185 GHz (0.48257208, 0.0)
#    88.63394 GHz (0.3509686, 0.0)
#    88.6316024 GHz (0.49761587, 0.0)

from FITS_tools import regrid_cube
hnc_regrid = regrid_cube(cubehnc.filled_data[:], cubehnc.header, cubehcn.header)
cubehnc = SpectralCube(data=hnc_regrid, wcs=cubehcn.wcs)

#cubehnc._wcs = cubehcn.wcs
#cubehnc.mask._wcs = cubehcn.wcs

mask  = (cubehcn > 0.5)
mask2 = (cubehcn > 0.5)
#mask2._wcs = cubehnc.wcs
hcn_flat = cubehcn.with_mask(mask).flattened()
hnc_flat = cubehnc.with_mask(mask2).flattened()

mask3 = cubehnc > 0.5
mask4 = cubehnc > 0.5
#mask4._wcs = cubehcn.wcs
hnc_flat2 = cubehnc.with_mask(mask3).flattened()
hcn_flat2 = cubehcn.with_mask(mask4).flattened()

import pylab as pl
pl.clf()
pl.plot(hnc_flat, hcn_flat, ',', alpha=0.5)
pl.plot(hnc_flat2, hcn_flat2, ',', alpha=0.5)
pl.plot([0,5],[0,5],'k--')
pl.axis([0,5,0,5])

dd = cubehcn.to_ds9()

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
def calc_irlum(catalog = 'cohrs_finalcatalog_clumpflux_medaxis.fits',
               refresh=False):
    ctrr = 0
    cat = Table.read(catalog)
    current_open_file = ''
    if 'ir_luminosity' not in cat.keys():
        keylist = ['ir_luminosity','ir_flux','ir_flux_short',
                'ir_lum_short','bg_flux','bg_lum']
        for thiskey in keylist:
            c = Column(np.zeros(len(cat))+np.nan,name=thiskey)
            cat.add_column(c)

    for cloud in cat:
        if np.isnan(cloud['ir_luminosity']) or refresh:
            orig_file = cloud['orig_file']+'.fits'
            asgn_file = cloud['orig_file']+'_fasgn.fits'
            higal_file = orig_file.replace('cohrs','higal_xmatch')
            if os.path.isfile(datadir+'COHRS/'+orig_file) and \
                    os.path.isfile(datadir+'HIGAL_MATCHED/'+higal_file):
                if current_open_file != datadir+'COHRS/'+orig_file:
                    hdu = fits.open(datadir+'COHRS/'+orig_file,memmap=False)
                    w = wcs.WCS(hdu[0].header)
                    hdr2 = w.to_header()
                    hdr2['BMAJ'] = 15./3600
                    hdr2['BMIN'] = 15./3600
                    hdr2['BPA'] = 0.
                    co = SpectralCube(hdu[0].data,w,header=hdr2)
                    irfull= fits.open(datadir+'HIGAL_MATCHED/'+higal_file,memmap=False)
                    irlong = fits.open(datadir+'HIGAL_MATCHED2/'+higal_file,memmap=False)
                    irmap = (irfull[0].data-irlong[0].data)
                    irmap2 = irfull[0].data
                    asgn = fits.getdata(datadir+'ASSIGNMENTS/'+asgn_file,memmap=False)
                    masked_co = co.with_mask(asgn>0*u.dimensionless_unscaled)
                    moment = masked_co.moment(0)
                    current_open_file = datadir+'COHRS/'+orig_file
                    cat.write('output_catalog.fits',overwrite=True)
#                mask = np.zeros(co.shape,dtype=np.bool)
#                mask[asgn == cloud['cloud_id']]=True
                cloud_cube = co.with_mask(asgn == cloud['cloud_id'])
                cloud_moment = cloud_cube.moment(0)
                cloud_cube = 0.0
                fraction = (cloud_moment.value/moment.value)
                planemask = skm.binary_closing(fraction > 0,selem=skm.disk(3))
                fraction = np.nanmean(fraction)
                rind = (skm.binary_dilation(planemask,selem=skm.disk(6)) ^\
                        skm.binary_dilation(planemask,selem=skm.disk(3)))*\
                    np.isfinite(irmap2)
                if np.any(rind):
                    rindvals = irmap2[rind]
                    clipval = 4*np.percentile(rindvals,15.87)-\
                              3*(np.percentile(rindvals,2.28))
                    rind *= irmap2 <= clipval
                    yv,xv = np.where(rind)
                    x0,y0 = np.median(xv),np.median(yv)
                    dataz = np.c_[np.ones(xv.size), 
                                  xv-x0, 
                                  yv-y0]
                    try:

                        lsqcoeffs,_,_,_ = np.linalg.lstsq(dataz,irmap2[rind])
                        outputs = lsq(myplane,np.r_[lsqcoeffs,0],
                                     args=(xv-x0,
                                           yv-y0,
                                           irmap2[rind]),
                                     loss = 'soft_l1')
                        coeffs = outputs.x
                        yhit,xhit = np.where(planemask)
                        bg = coeffs[0]+coeffs[1]*(xhit-x0)+\
                             coeffs[2]*(yhit-y0)+coeffs[3]*(yhit-y0)*(xhit-x0)

                        # I am sitcking a 6e11 in here as the frequency of 
                        # the 500 microns
                        bgavg = np.sum(fraction*bg)/6e11
                        bglum = bgavg*cloud['distance']**2*\
                                3.086e18**2*np.pi*4/3.84e33

                    except ValueError:
                        import pdb; pdb.set_trace()


                    ir_flux = np.nansum(fraction*(irmap[planemask]))/6e11
                    ir_lum = ir_flux * cloud['distance']**2*\
                             3.086e18**2*np.pi*4/3.84e33
                    ir_flux2 = np.nansum(fraction*irmap2[planemask])/6e11
                    ir_lum2 = ir_flux2 * cloud['distance']**2*\
                             3.086e18**2*np.pi*4/3.84e33
                    cloud['ir_flux'] = ir_flux2
                    cloud['ir_luminosity'] = ir_lum2
                    cloud['ir_flux_short'] = ir_flux
                    cloud['ir_lum_short'] = ir_lum
                    cloud['bg_flux'] = bgavg
                    cloud['bg_lum'] = bglum
                    print(cloud['cloud_id'],ir_flux,ir_lum,ir_lum2,bglum)
                    cat.write('output_catalog.fits',overwrite=True)
                    ctrr+=1
                    if ctrr==20:
                        return(cat)
    #                if cloud['volume_pc2_kms']>1e2:
    #                    import pdb; pdb.set_trace()
    
#.........这里部分代码省略.........

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
    denscube = denscube.with_mask(okmask)
    colcube  = SpectralCube.read(paths.dpath("H2CO_ParameterFits_{0}col.fits".format(meastype)))
    colcube = colcube.with_mask(okmask)
    goodmask_std = stdcube < 0.5

    flathead = denscube.wcs.dropaxis(2).to_header()

    denscube_lin = SpectralCube(10**denscube.filled_data[:].value,
                                              wcs=denscube.wcs,
                                              mask=okmask)
    colcube_lin = SpectralCube(10**colcube.filled_data[:].value,
                                             wcs=denscube.wcs, mask=okmask)


    totalcol = colcube_lin.sum(axis=0)
    totalcolgood = colcube_lin.with_mask(goodmask).sum(axis=0)
    totalcolgoodstd = colcube_lin.with_mask(goodmask_std).sum(axis=0)

    denscol = SpectralCube(denscube_lin.filled_data[:] * colcube_lin.filled_data[:], wcs=denscube.wcs,
                                         mask=okmask)
    wtdmeandens = np.log10(denscol.sum(axis=0) / totalcol)

    mindens_std = denscube.with_mask(goodmask_std).min(axis=0)
    mindens_chi2 = denscube.with_mask(goodmask).min(axis=0)

    hdu1 = fits.PrimaryHDU(data=wtdmeandens.value,
                           header=flathead)
    hdu1.writeto(paths.dpath("H2CO_ParameterFits_weighted_mean_{0}_density.fits".format(meastype)), clobber=True)

    masked_wtdmeans = np.log10(denscol.with_mask(goodmask).sum(axis=0) / totalcolgood)
    hdu2 = fits.PrimaryHDU(data=masked_wtdmeans.value,

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import with_mask [as 别名]
add_noise = sys.argv[3]
if add_noise == "True" or add_noise == "T":
    add_noise = True
else:
    add_noise = False


data1, hdr1 = fits.getdata(fits1, header=True)
data2, hdr2 = fits.getdata(fits2, header=True)

if add_noise:
    data1 += np.random.normal(0.0, 0.788 / 10, data1.shape)
    data2 += np.random.normal(0.0, 0.788 / 10, data2.shape)

cube1 = SpectralCube(data=data1, wcs=WCS(hdr1))
cube1 = cube1.with_mask(LazyMask(np.isfinite, cube1))
cube2 = SpectralCube(data=data2, wcs=WCS(hdr2))
cube2 = cube2.with_mask(LazyMask(np.isfinite, cube2))

set1 = Mask_and_Moments(cube1)
mask = cube1 > 3*set1.scale
set1.make_mask(mask=mask)
set1.make_moments()
set1.make_moment_errors()
dataset1 = set1.to_dict()

set2 = Mask_and_Moments(cube2)
mask = cube2 > 3*set2.scale
set2.make_mask(mask=mask)
set2.make_moments()
set2.make_moment_errors()