Python SpectralCube.sum方法代码示例

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import sum [as 别名]

#.........这里部分代码省略.........
        update_low_hdr = \
            _update_beam_in_hdr(self.lowres.header, self.combined_beam)

        self.lowres_convolved = \
            SpectralCube(lowres_convolved*self.lowres.unit, self.lowres.wcs,
                         header=update_low_hdr)

    def flux_recovered(self, plot=True, filename=None, enable_interp=True,
                       interp_to='lower', diff_tol=1e-4*u.m/u.s):
        '''
        Check the amount of flux recovered in the high resolution image versus
        the low resolution data. If the spectral axes don't match, one is
        interpolated onto the other.

        Parameters
        ----------
        plot : bool, optional
            Enable plotting.
        filename : str, optional
            Give filename for the plot save file. When specified, the plot
            is automatically saved.
        interp_to : 'lower' or 'upper', optional
            If the spectral axes don't match, interpolated onto the same.
            The default 'lower' interpolates to the spectral axis with the
            lowest resolution.
        '''

        # Add up the total intensity in the cubes and compare
        # Assumes that there is some masking such that noise
        # doesn't dominate

        if self.highres_convolved is not None:
            high_channel_intensity = \
                self.highres_convolved.sum(axis=(1, 2))
        else:
            high_channel_intensity = self.highres.sum(axis=(1, 2))
            Warning("Should run convolve_to_common before. Using unconvolved"
                    " cube.")

        if self.lowres_convolved is not None:
            low_channel_intensity = self.lowres_convolved.sum(axis=(1, 2))
        else:
            low_channel_intensity = self.lowres.sum(axis=(1, 2))
            Warning("Should run convolve_to_common before. Using unconvolved"
                    " cube.")

        # If the spectral axes are not the same, try interpolating to match
        if enable_interp:

            better_vres_high = \
                np.abs(self.lowres.header['CDELT3']) > \
                np.abs(self.highres.header['CDELT3'])

            spec_diff = np.abs(np.abs(self.lowres.header['CDELT3']) -
                               np.abs(self.highres.header['CDELT3']))

            spec_unit = self.highres.spectral_axis.unit

            if spec_diff * spec_unit < diff_tol:
                warnings.warn("Channels are essentially the same. "
                              "Skipping interpolation. Lower diff_tol"
                              " to re-enable interpolation.")

            else:
                # Invert which to interpolate to
                if interp_to == "higher":

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import sum [as 别名]
    denscube = SpectralCube.read(paths.dpath("H2CO_ParameterFits_{0}dens.fits".format(meastype)))
    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)