Python SpectralCube.spectral_slab方法代码示例

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import spectral_slab [as 别名]
def deblend_cube(region='OrionA',vmin=0.,vmax=20.):
  tex = fits.getdata('{0}/parameterMaps/{0}_Tex_DR1_rebase3_flag.fits'.format(region))
  mom0 = fits.getdata('{0}/{0}_NH3_11_DR1_rebase3_mom0_QA_trim.fits'.format(region))
  vlsr = fits.getdata('{0}/parameterMaps/{0}_Vlsr_DR1_rebase3_flag.fits'.format(region))
  sigv = fits.getdata('{0}/parameterMaps/{0}_Sigma_DR1_rebase3_flag.fits'.format(region))
  nnh3 = fits.getdata('{0}/parameterMaps/{0}_N_NH3_DR1_rebase3_flag.fits'.format(region))
  cube = SpectralCube.read('{0}/{0}_NH3_11_DR1_rebase3_trim.fits'.format(region))
  cube  = cube.with_spectral_unit(u.km/u.s,velocity_convention='radio')
  
  tpeak = mom0/(np.sqrt(2*np.pi)*sigv)

  vlsr[vlsr==0]=np.nan
  sigv[sigv==0]=np.nan
  
  deblend = np.zeros(cube.shape)
  hdr = cube.wcs.to_header()
  spaxis = cube.spectral_axis.value

  for plane in np.arange(deblend.shape[0]):
      deblend[plane,:,:] = tpeak*np.exp(-(spaxis[plane]-vlsr)**2/(2*sigv**2))

  newcube = SpectralCube(deblend,cube.wcs,header=hdr)
  slab = newcube.spectral_slab(vmin*u.km/u.s,vmax*u.km/u.s)
  slab.write('{0}/{0}_singlecomp.fits'.format(region),overwrite=True)  

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import spectral_slab [as 别名]
                        strides=(0,)+noise.strides)
noise_spcube = SpectralCube(data=noise_cube, wcs=cube_merge_high.wcs)

cube_merge_high_sm = SpectralCube.read(mpath('APEX_H2CO_merge_high_plait_all_smooth.fits'))
noise_sm = fits.getdata(mpath('APEX_H2CO_merge_high_plait_all_smooth_noise.fits'))
noise_cube_sm = as_strided(noise_sm, shape=cube_merge_high_sm.shape,
                           strides=(0,)+noise_sm.strides)
noise_spcube_sm = SpectralCube(data=noise_cube_sm, wcs=cube_merge_high_sm.wcs)

# Create a cutout of the cube covering the H2CO lines
# it's barely worth it; cuts off 10% of pixels
f1 = all_lines['H2CO_303_202']*u.GHz
f2 = all_lines['H2CO_321_220']*u.GHz
h2co_cube_merge_high = cube_merge_high.spectral_slab(f1*(1-(150*u.km/u.s/constants.c)),
                                                     f2*(1+(100*u.km/u.s/constants.c)))
h2co_noise_cube = noise_spcube.spectral_slab(f1*(1-(150*u.km/u.s/constants.c)),
                                           f2*(1+(100*u.km/u.s/constants.c)))

h2co_cube_merge_high_sm = cube_merge_high_sm.spectral_slab(f1*(1-(150*u.km/u.s/constants.c)),
                                                           f2*(1+(100*u.km/u.s/constants.c)))
h2co_noise_cube_sm = noise_spcube_sm.spectral_slab(f1*(1-(150*u.km/u.s/constants.c)),
                                                 f2*(1+(100*u.km/u.s/constants.c)))


# Pyspeckit cube made from spectralcube
pcube_merge_high = pyspeckit.Cube(cube=h2co_cube_merge_high._data,
                                  errorcube=h2co_noise_cube._data,
                                  header=h2co_cube_merge_high.header,
                                  xarr=h2co_cube_merge_high.spectral_axis,
                                 )
pcube_merge_high.xarr.refX = 218.22219
pcube_merge_high.xarr.refX_unit = 'GHz'

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import spectral_slab [as 别名]
                          header=high1e3dens_co_slab3_mom0.hdu.header)

# Purely detection-based thresholding
h2co11 = SpectralCube.read(h2co11subfn)
h2co22 = SpectralCube.read(h2co22subfn)

noisevrange = [-50,0]*u.km/u.s
h2co11noiseslab = h2co11.spectral_slab(*noisevrange)
h2co22noiseslab = h2co22.spectral_slab(*noisevrange)
h2co11noise = h2co11noiseslab.apply_numpy_function(np.std,axis=0)
h2co22noise = h2co22noiseslab.apply_numpy_function(np.std,axis=0)

# TODO: implement Cube (lazy?) Arithmetic: avoid filling data!
sn11 = SpectralCube(np.abs(h2co11.filled_data[:])/h2co11noise, wcs=h2co11.wcs)
sn22 = SpectralCube(np.abs(h2co22.filled_data[:])/h2co22noise, wcs=h2co22.wcs)
sn11_slab3 = sn11.spectral_slab(vrange1[0], vrange2[1])
sn22_slab3 = sn22.spectral_slab(vrange1[0], vrange2[1])

h2co11_total = np.array([cube13ss_slab3.with_mask(sn11_slab3>threshold).sum().value
                         for threshold in np.arange(1,6)])
h2co22_total = np.array([cube13ss_slab3.with_mask(sn22_slab3>threshold).sum().value
                         for threshold in np.arange(1,6)])
h2co11_fraction = h2co11_total / total_co_slab.value
h2co22_fraction = h2co22_total / total_co_slab.value
both_fraction = cube13ss_slab3.with_mask((sn11_slab3 > 2) &
                                         (sn22_slab3 > 2)).sum().value / total_co_slab.value

log.info("H2CO Fractions: 1-1: {0}".format(h2co11_fraction))
log.info("H2CO Fractions: 2-2: {0}".format(h2co22_fraction))
log.info("Both: {0}".format(both_fraction))

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import spectral_slab [as 别名]
    makefig(hdu5, 4, name='stdmasked_weighted_mean_{0}fit_density'.format(meastype), weight_hdu=weight.hdu)
    makefig(hdu6, 5, name='stdmasked_min_{0}fit_density'.format(meastype), weight_hdu=weight.hdu)
    makefig(hdu7, 6, name='mean_{0}fit_abundance'.format(meastype),
            vmin_lin=10**-14, vmax_lin=10**-6, vmin_log=-14, vmax_log=-6,
            labeltype='abundance', linscale=1, weight_hdu=weight.hdu)
    makefig(hdu8, 7, name='total_{0}fit_column'.format(meastype),
            vmin_lin=10**12, vmax_lin=10**16, vmin_log=12, vmax_log=16,
            labeltype='column', linscale=1, weight_hdu=weight.hdu)

    nextfig = 8

    for ii,vrange in enumerate((vrange1,vrange2)):
        vrange = vrange[0]*1.01, vrange[1]*0.99
        log.info("Integrating over {0}".format(vrange))

        wtdmeandens = np.log10(denscol.spectral_slab(*vrange).sum(axis=0) /
                               colcube_lin.spectral_slab(*vrange).sum(axis=0))
        hdu1 = fits.PrimaryHDU(data=wtdmeandens.value,
                                header=denscube.wcs.dropaxis(2).to_header())
        hdu1.writeto(paths.dpath("H2CO_ParameterFits_weighted_mean_{2}_density_v{0}to{1}.fits".format(vrange[0].value,
                                                                                                      vrange[1].value,
                                                                                                      meastype)),
                     clobber=True)

        masked_wtdmeans = np.log10(denscol.with_mask(goodmask).spectral_slab(*vrange).sum(axis=0) /
                                   colcube_lin.with_mask(goodmask).spectral_slab(*vrange).sum(axis=0))
        hdu2 = fits.PrimaryHDU(data=masked_wtdmeans.value,
                               header=denscube.wcs.dropaxis(2).to_header())
        hdu2.writeto(paths.dpath("H2CO_ParameterFits_weighted_mean_{2}_density_chi2masked_v{0}to{1}.fits".format(vrange[0].value,
                                                                                                                 vrange[1].value,
                                                                                                                 meastype)),