Python SpectralCube.write方法代码示例

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import write [as 别名]
def measure_dendrogram_properties(dend=None, cube303=cube303,
                                  cube321=cube321, cube13co=cube13co,
                                  cube18co=cube18co, noise_cube=noise_cube,
                                  sncube=sncube,
                                  suffix="",
                                  last_index=None,
                                  plot_some=True,
                                  line='303',
                                  write=True):

    assert (cube321.shape == cube303.shape == noise_cube.shape ==
            cube13co.shape == cube18co.shape == sncube.shape)
    assert sncube.wcs is cube303.wcs is sncube.mask._wcs

    metadata = {}
    metadata['data_unit'] = u.K
    metadata['spatial_scale'] =  7.2 * u.arcsec
    metadata['beam_major'] =  30 * u.arcsec
    metadata['beam_minor'] =  30 * u.arcsec
    metadata['wavelength'] =  218.22219*u.GHz
    metadata['velocity_scale'] = u.km/u.s
    metadata['wcs'] = cube303.wcs

    keys = [
            'density_chi2',
            'expected_density',
            'dmin1sig_chi2',
            'dmax1sig_chi2',
            'column_chi2',
            'expected_column',
            'cmin1sig_chi2',
            'cmax1sig_chi2',
            'temperature_chi2',
            'expected_temperature',
            'tmin1sig_chi2',
            'tmax1sig_chi2',
            'eratio321303',
            'ratio321303',
            'logh2column',
            'elogh2column',
            'logabundance',
            'elogabundance',
           ]
    obs_keys = [
            'Stot303',
            'Smin303',
            'Smax303',
            'Stot321',
            'Smean303',
            'Smean321',
            'npix',
            'e303',
            'e321',
            'r321303',
            'er321303',
            '13cosum',
            'c18osum',
            '13comean',
            'c18omean',
            's_ntotal',
            'index',
            'is_leaf',
            'parent',
            'root',
            'lon',
            'lat',
            'vcen',
            'higaldusttem',
            'reff',
            'dustmass',
            'dustmindens',
            'bad',
            #'tkin_turb',
    ]
    columns = {k:[] for k in (keys+obs_keys)}

    log.debug("Initializing dendrogram temperature fitting loop")

    # FORCE wcs to match
    # (technically should reproject here)
    cube13co._wcs = cube18co._wcs = cube303.wcs
    cube13co.mask._wcs = cube18co.mask._wcs = cube303.wcs

    if line == '303':
        maincube = cube303
    elif line == '321':
        maincube = cube321
    else:
        raise ValueError("Unrecognized line: {0}".format(line))

    # Prepare an array to hold the fitted temperatures
    tcubedata = np.empty(maincube.shape, dtype='float32')
    tcubedata[:] = np.nan
    tcubeleafdata = np.empty(maincube.shape, dtype='float32')
    tcubeleafdata[:] = np.nan


    nbad = 0

    catalog = ppv_catalog(dend, metadata)
#.........这里部分代码省略.........

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import write [as 别名]
                                                               column_flat, utline303,
                                                               utline321, unoise)):
    if tcube[z,y,x] == 0:
        logh2column = np.log10(col)+22

        mf.set_constraints(ratio303321=rat, eratio303321=erat,
                           #ratio321322=ratio2, eratio321322=eratio2,
                           logh2column=logh2column, elogh2column=elogh2column,
                           logabundance=logabundance, elogabundance=elogabundance,
                           taline303=ta303.value, etaline303=err,
                           taline321=ta321.value, etaline321=err,
                           linewidth=linewidth)
        row_data = mf.get_parconstraints()
        tcube[z,y,x] = row_data['temperature_chi2']
        row_data['ratio303321'] = rat
        row_data['eratio303321'] = erat

        if ii % 100 == 0 or ii < 50:
            log.info("T: [{tmin1sig_chi2:7.2f},{temperature_chi2:7.2f},{tmax1sig_chi2:7.2f}]  R={ratio303321:6.2f}+/-{eratio303321:6.2f}".format(**row_data))
        else:
            pb.update(ii)
        tcube.flush()
    else:
        pb.update(ii)

tcube[tcube==0] = np.nan
tCube = SpectralCube(tcube, cube303.wcs, mask=BooleanArrayMask(np.isfinite(tcube), wcs=cube303.wcs))
tCube.write(hpath('chi2_temperature_cube.fits'), overwrite=True)

print()

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import write [as 别名]
def rebaseline(filename, blorder=3,
               baselineRegion=[slice(0, 262, 1), slice(-512, 0, 1)],
               windowFunction=None, blankBaseline=False,
               flagSpike=True, v0=None, **kwargs):
    """
    Rebaseline a data cube using robust regression of Legendre polynomials.

    Parameters
    ----------
    filename : string
        FITS filename of the data cube
    blorder : int
        Order of the polynomial to fit to the data
    baselineRegion : list
        List of slices defining the default region of the spectrum, in
        channels, to be used for the baseline fitting.
    windowFunction : function
        Name of function to be used that will accept spectrum data, and
        velocity axis and will return a binary mask of the channels to be used
        in the  baseline fitting.  Extra **kwargs are passed to windowFunction
        to do with as it must.
    blankBaseline : boolean
        Blank the baseline region on a per-spectrum basis

    Returns
    -------
    Nothing.  A new FITS file is written out with the suffix 'rebaseN' where N
    is the baseline order

    """
    cube = SpectralCube.read(filename)
    originalUnit = cube.spectral_axis.unit
    cube = cube.with_spectral_unit(u.km / u.s, velocity_convention='radio')
    spaxis = cube.spectral_axis.to(u.km / u.s).value

    goodposition = np.isfinite(cube.apply_numpy_function(np.max, axis=0))
    y, x = np.where(goodposition)
    outcube = np.zeros(cube.shape) * np.nan
    RegionName = (filename.split('_'))[0]

    if hasattr(windowFunction, '__call__'):
        catalog = catalogs.GenerateRegions()

    nuindex = np.arange(cube.shape[0])
    runmin = nuindex[-1]
    runmax = nuindex[0]

    for thisy, thisx in console.ProgressBar(zip(y, x)):
        spectrum = cube[:, thisy, thisx].value

        if v0 is not None:
            baselineIndex = windowFunction(spectrum, spaxis,
                                           v0=v0, **kwargs)
        elif hasattr(windowFunction, '__call__'):
            _, Dec, RA = cube.world[0, thisy, thisx]
            # This determines a v0 appropriate for the region
            v0 = VlsrByCoord(RA.value, Dec.value, RegionName,
                             regionCatalog=catalog)
            baselineIndex = windowFunction(spectrum, spaxis,
                                           v0=v0, **kwargs)
        else:
            baselineIndex = np.zeros_like(spectrum,dtype=np.bool)
            for ss in baselineRegion:
                baselineIndex[ss] = True

        runmin = np.min([nuindex[baselineIndex].min(), runmin])
        runmax = np.max([nuindex[baselineIndex].max(), runmax])

        # Use channel-to-channel difference as the noise value.
        if flagSpike:
            jumps = (spectrum - np.roll(spectrum, -1))
            noise = mad1d(jumps) * 2**(-0.5)
            baselineIndex *= (np.abs(jumps) < 5 * noise)
            noise = mad1d((spectrum -
                           np.roll(spectrum, -2))[baselineIndex]) * 2**(-0.5)
        else:
            noise = mad1d((spectrum -
                           np.roll(spectrum, -2))[baselineIndex]) * 2**(-0.5)

        if blankBaseline:
            spectrum = robustBaseline(spectrum, baselineIndex,
                                      blorder=blorder,
                                      noiserms=noise)
            spectrum[baselineIndex] = np.nan
            outcube[:, thisy, thisx] = spectrum
        else:
            outcube[:, thisy, thisx] = robustBaseline(spectrum, baselineIndex,
                                                      blorder=blorder,
                                                      noiserms=noise)

    outsc = SpectralCube(outcube, cube.wcs, header=cube.header)
    outsc = outsc[runmin:runmax, :, :]  # cut beyond baseline edges
    # Return to original spectral unit
    outsc = outsc.with_spectral_unit(originalUnit)
    outsc.write(filename.replace('.fits', '_rebase{0}.fits'.format(blorder)),
                overwrite=True)

# 需要导入模块: from spectral_cube import SpectralCube [as 别名]
# 或者: from spectral_cube.SpectralCube import write [as 别名]
        rcubedata[structure.get_mask()] = r321303
        tcubedata[structure.get_mask()] = pwtem(np.array([r321303]))

        pb.update(ii+1)

    # Note that there are overlaps in the catalog, which means that ORDER MATTERS
    # in the above loop.  I haven't yet checked whether large scale overwrites
    # small or vice-versa; it may be that both views of the data are interesting.
    tcube = SpectralCube(data=tcubedata, wcs=cubeA.wcs,
                         mask=cubeA.mask, meta={'unit':'K'},
                         header=cubeA.header,
                        )

    outpath = 'TemperatureCube_DendrogramObjects{0}_Piecewise.fits'.format(sm)
    tcube.write(hpath(outpath), overwrite=True)

    rcube = SpectralCube(data=rcubedata, wcs=cubeA.wcs,
                         mask=cubeA.mask, meta={'unit':'K'},
                         header=cubeA.header,
                        )

    outpath = 'RatioCube_DendrogramObjects{0}.fits'.format(sm)
    rcube.write(hpath(outpath), overwrite=True)

    max_temcube = tcube.max(axis=0)
    max_temcube.hdu.writeto(hpath('TemperatureCube_DendrogramObjects{0}_Piecewise_max.fits'.format(sm)), clobber=True)
    max_rcube = rcube.max(axis=0)
    max_rcube.hdu.writeto(hpath('RatioCube_DendrogramObjects{0}_Piecewise_max.fits'.format(sm)), clobber=True)

    mean_temcube = tcube.mean(axis=0)