Использование healpy.anafast на картах гамма-лучей

Question

У меня есть карты гамма-излучения (изображение с поверхностной яркостью) в формате Fits, а также.hpx, выводимые конвертером Aladin.

Я хочу вычислить угловой спектр мощности. Как создать файл для чтенияhealpy.anafast? Кажется, я неправильно понимаю формат данных (TypeErrors).

Одним из изображений гамма-излучения, которое я пробовал, был Fermi Galactic Diffuse. Файл является общедоступной Галактической диффузной картой LAT с именем gll_iem_v02_P6_V11_DIFFUSE.fit на:

http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html

Я вставил код ниже, как я его использую, но это по сути сценарий под названием plot_wmap_power_spectra на астромле:

    """
    WMAP power spectrum analysis with HealPy
    ----------------------------------------

    This demonstrates how to plot and take a power spectrum of the WMAP data
    using healpy, the python wrapper for healpix.  Healpy is available for
    download at the `github site <https://github.com/healpy/healpy>`_
    """
    # Author: Jake VanderPlas <[e-mail]>
    # License: BSD
    #   The figure is an example from astroML: see http://astroML.github.com
    import numpy as np
    from matplotlib import pyplot as plt

    # warning: due to a bug in healpy, importing it before pylab can cause
    #  a segmentation fault in some circumstances.
    import pylab
    import healpy as hp
    ###
    from astroML.datasets import fetch_wmap_temperatures
    ###

    #------------------------------------------------------------
    # Fetch the data
    ###
    wmap_unmasked = fetch_wmap_temperatures(masked=False)

    #PredictedSurfaceFluxFromModelMap =           np.arange(hp.read_map('PredictedSurfaceFluxFromModelMap.hpx[1]'))
    PredictedSurfaceFluxFromModelMap =   hp.read_map('gll_iem_v02_p6_V11_DIFFUSE.fit',dtype=np.float,verbose=True)
    #PredictedSurfaceFluxFromModelMap = hp.read_map('all.fits',dtype=np.float,verbose=True)
    #cl_out =  hp.read_cl('PredictedSurfaceFluxFromModelMap.hpx',dtype=np.float)#,verbose=True)
    wmap_masked = fetch_wmap_temperatures(masked=True)
    ###
    white_noise = np.ma.asarray(np.random.normal(0, 0.062, wmap_masked.shape))
    len(cl_out)

    #------------------------------------------------------------
    # plot the unmasked map
    fig = plt.figure(1)
    #hp.mollview(wmap_unmasked, min=-1, max=1, title='Unmasked map',
    #            fig=1, unit=r'$\Delta$T (mK)')
    ########----------------
    ##hp.mollview(PredictedSurfaceFluxFromModelMap, min=-1, max=1, title='Unmasked map',
    ##            fig=1, unit=r'$\Delta$T (mK)')
    ########----------------
    #------------------------------------------------------------
    # plot the masked map
    #  filled() fills the masked regions with a null value.
    ########----------------
    #fig = plt.figure(2)
    #hp.mollview(wmap_masked.filled(), title='Masked map',
    #            fig=2, unit=r'$\Delta$T (mK)')

    ########----------------
    #------------------------------------------------------------
    # compute and plot the power spectrum
    ########----------------
    #cl = hp.anafast(wmap_masked.filled(), lmax=1024)
    cl = hp.anafast(PredictedSurfaceFluxFromModelMap, lmax=1024)
    #cl = cl_out
    ########----------------
    ell = np.arange(len(cl))

    cl_white = hp.anafast(white_noise, lmax=1024)

    fig = plt.figure(3)
    ax = fig.add_subplot(111)
    ax.scatter(ell, ell * (ell + 1) * cl,
       s=4, c='black', lw=0,
       label='data')
     ax.scatter(ell, ell * (ell + 1) * cl_white,
       s=4, c='gray', lw=0,
       label='white noise')

     ax.set_xlabel(r'$\ell$')
     ax.set_ylabel(r'$\ell(\ell+1)C_\ell$')
     ax.set_title('Angular Power (not mask corrected)')
     ax.legend(loc='upper right')
     ax.grid()
     ax.set_xlim(0, 1100)

     plt.show()

Answer 1

Я загрузил вашу карту также в Figshare, где, вероятно, будет доступно в будущем.

Если у вас есть карта в формате HEALPix, ее легко прочитать с помощью healpy:

import healpy as hp
m = hp.ma(hp.read_map("gll_iem_v02_p6_V11_DIFFUSE.hpx"))

Маска NaN пикселей:

m.mask = np.isnan(m)

Сюжет:

hp.mollview(m, min=-1e-5, max=1e-5, xsize=2000)
title("gll_iem_v02_p6_V11_DIFFUSE")

Вычислить и построить график угловой мощности:

plt.loglog(hp.anafast(m))

Смотрите также блокнот IPython: http://nbviewer.ipython.org/7553252