How to use the healpy.alm2map function in healpy

To help you get started, we’ve selected a few healpy examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github damonge / CoLoRe / test / study_pk.py View on Github external
tbdata = hdulist[1].data
    
    pix=hp.ang2pix(nside,(90-tbdata['DEC'])*np.pi/180,tbdata['RA']*np.pi/180)
    n=np.bincount(pix,minlength=npix)
    e1=np.bincount(pix,minlength=npix,weights=tbdata['E1'])
    e2=np.bincount(pix,minlength=npix,weights=tbdata['E2'])
    nmap+=n
    e1map+=e1
    e2map+=e2
    ifile+=1

ndens=(np.sum(nmap)+0.0)/(4*np.pi)
mp_e1=e1map/nmap; mp_e1[nmap<=0]=0
mp_e2=e2map/nmap; mp_e2[nmap<=0]=0
mp_d=(nmap+0.0)/np.mean(nmap+0.0)-1
mp_db,mp_E,mp_B=hp.alm2map(hp.map2alm(np.array([mp_d,mp_e1,mp_e2]),pol=True),pol=False,nside=nside)

lt,cls_dd=np.loadtxt('test/outlj_cl_dd.txt',unpack=True);
lt,clt_dl=np.loadtxt('test/outlj_cl_d1l2.txt',unpack=True);
lt,clt_ll=np.loadtxt('test/outlj_cl_ll.txt',unpack=True);
lt,clt_kd=np.loadtxt('test/outlj_cl_dc.txt',unpack=True);
lt,clt_kk=np.loadtxt('test/outlj_cl_cc.txt',unpack=True);
lt,clt_id=np.loadtxt('test/outlj_cl_di.txt',unpack=True);
lt,clt_ii=np.loadtxt('test/outlj_cl_ii.txt',unpack=True);
cln_dd=np.ones_like(lt)/ndens
clt_dd=cls_dd+cln_dd
cld_dd,cld_ee,cld_bb,cld_de,cld_eb,cld_db=hp.anafast(np.array([mp_d,mp_e1,mp_e2]),pol=True);
ld=np.arange(len(cld_dd));

#Analyze kappa
mp_k=hp.read_map("test/out_kappa_z000.fits")
cld_kk=hp.anafast(mp_k); ld=np.arange(len(cld_kk))
github astro-informatics / s2let / src / main / python / pys2let_test_directional.py View on Github external
# Call pys2let and compute directional wavelet transform.
# Returns an array of MW maps
# The way to access them is described below.
print 'Running analysis_lm2wav'
f_wav, f_scal = analysis_lm2wav(f_lm, B, L, J_min, N, spin, upsample)
print 'Done'


print 'size f_scal f_wav', f_scal.size, f_wav.size, f_wav.size / f_scal.size

# Uses synthesis to reconstruct the input alms.
print 'Running synthesis_wav2lm'
f_lm_rec = synthesis_wav2lm(f_wav, f_scal, B, L, J_min, N, spin, upsample)
print 'Done'
f_rec = hp.alm2map(f_lm_rec, nside=nside, lmax=L-1)
#plot to compare the input/output maps
#hp.mollview(f)
#hp.mollview(f_rec)

# Convert to MW sampling from spherical harmonics
f_mw = alm2map_mw(f_lm, L, spin)
f_lm_rec = map2alm_mw(f_mw, L, spin) # To check that the analysis routine works
f_mw_rec = alm2map_mw(f_lm, L, spin)


# Home made plotting routine! inputs : function f (1D array of MW signal), bandlimit L, plot axis ax, and title
def myplot(f, L, ax, title=''):
    thetas, phis = mw_sampling(L) # Compute the theta and phi values of the MW equiangular grid.
    ntheta = len(thetas)
    nphi = len(phis)
    arr = f.reshape((ntheta, nphi)) # Convert the input MW 1D array into 2D array with rows for theta and columns for phi. As simple as that!
github astro-informatics / s2let / src / main / python / pys2let_test_axisym_hpx.py View on Github external
import matplotlib.pyplot as plt

nside = 128
L = 128
J_min = 1
B = 3
J = pys2let_j_max(B, L, J_min)

# The filename of some random healpix map
fname = '/Users/bl/Dropbox/Wavelets/s2let/data/somecmbsimu_hpx_128.fits'

# Read healpix map and compute alms. 
# f_lm has size L*(L+1)/2
f_ini = hp.read_map(fname) # Initial map
f_lm = hp.map2alm(f_ini, lmax=L-1) # Its alms
f = hp.alm2map(f_lm, nside=nside, lmax=L-1) # Band limited version

hp.mollview(f)

# Call pys2let and compute wavelet transform. Returns the harmonic coefficients of the wavelets.
# f_scal_lm has size L*(L+1)/2
# f_wav_lm has size L*(L+1)/2 by J-J_min+1
f_wav_lm, f_scal_lm = analysis_axisym_lm_wav(f_lm, B, L, J_min)

# Reconstruct healpix maps on the sphere and plot them
f_scal = hp.alm2map(f_scal_lm, nside=nside, lmax=L-1)
hp.mollview(f_scal)
f_wav = np.empty([12*nside*nside, J-J_min+1])
for j in range(J-J_min+1):
	flm = f_wav_lm[:,j].ravel()
	f_wav[:,j] = hp.alm2map(flm, nside=nside, lmax=L-1)
	hp.mollview(f_wav[:,j])
github astro-informatics / s2let / src / main / python / pys2let_test_directional_manual.py View on Github external
J_min = 0
L_transitions = np.array([64])
L = 256
hybrid_scal_l, hybrid_wav_l, hybrid_scal_bandlimit, hybrid_wav_bandlimits, J, L_bounds = construct_hybrid_tiling(L, J_min, Bs, L_transitions)
hybrid_wav_l = hybrid_wav_l.T.ravel()
res = verify_tiling(L, hybrid_scal_l, hybrid_wav_l, hybrid_scal_bandlimit, hybrid_wav_bandlimits)
print 'Is the tiling OK and giving an invertible wavelet transform?', res


# The filename of some arbitrary healpix map
fname = '/Users/bl/Dropbox/Astrodata/SDSS/Fields/Planck_EBV_256rQ.fits'
#fname = '/Users/bl/Dropbox/Wavelets/s2let/data/somecmbsimu_hpx_128.fits'

f_ini = hp.read_map(fname) # Initial map
f_lm = hp.map2alm(f_ini, lmax=L-1) # Its alms
f = hp.alm2map(f_lm, nside=nside, lmax=L-1) # Band limited version

print 'Running analysis_lm2wav'
f_wav, f_scal = analysis_lm2wav_manualtiling(f_lm, L, N, spin, hybrid_scal_l, hybrid_wav_l, hybrid_scal_bandlimit, hybrid_wav_bandlimits)

# Uses synthesis to reconstruct the input alms.
print 'Running synthesis_wav2lm'
f_lm_rec = synthesis_wav2lm_manualtiling(f_wav, f_scal, L, N, spin, hybrid_scal_l, hybrid_wav_l, hybrid_scal_bandlimit, hybrid_wav_bandlimits)
f_rec = hp.alm2map(f_lm_rec, nside=nside, lmax=L-1)

# Convert to MW sampling from spherical harmonics
f_mw = alm2map_mw(lm_hp2lm(f_lm, L), L, spin)
#f_lm_rec = map2alm_mw(f_mw, L, spin) # To check that the analysis routine works
#f_mw_rec = alm2map_mw(f_lm_rec, L, spin)
f_mw_rec = alm2map_mw(lm_hp2lm(f_lm_rec, L), L, spin)
github lscsoft / lalsuite-archive / lalinference / python / lalinference / bayestar / postprocess.py View on Github external
m_rotated : np.ndarray
        The rotated HEALPix array.
    """
    npix = len(m)
    nside = hp.npix2nside(npix)

    theta = 0.5 * np.pi - dec
    phi = ra

    if method == 'fft':
        if nest:
            m = hp.reorder(m, n2r=True)
        alm = hp.map2alm(m)
        hp.rotate_alm(alm, -phi, -theta, 0.0)
        ret = hp.alm2map(alm, nside, verbose=False)
        if nest:
            ret = hp.reorder(ret, r2n=True)
    elif method == 'direct':
        R = hp.Rotator(
            rot=np.asarray([0, theta, -phi]),
            deg=False, inv=False, eulertype='Y')
        theta, phi = hp.pix2ang(nside, np.arange(npix), nest=nest)
        ipix = hp.ang2pix(nside, *R(theta, phi), nest=nest)
        ret = m[ipix]
    else:
        raise ValueError('Unrecognized method: {0}'.format(method))

    return ret
github lscsoft / lalsuite-archive / lalinference / python / lalinference / bayestar / postprocess.py View on Github external
m_rotated : np.ndarray
        The rotated HEALPix array.
    """
    npix = len(m)
    nside = hp.npix2nside(npix)

    theta = 0.5 * np.pi - dec
    phi = ra

    if method == 'fft':
        if nest:
            m = hp.reorder(m, n2r=True)
        alm = hp.map2alm(m)
        hp.rotate_alm(alm, -phi, -theta, 0.0)
        ret = hp.alm2map(alm, nside, verbose=False)
        if nest:
            ret = hp.reorder(ret, r2n=True)
    elif method == 'direct':
        R = hp.Rotator(
            rot=np.asarray([0, theta, -phi]),
            deg=False, inv=False, eulertype='Y')
        theta, phi = hp.pix2ang(nside, np.arange(npix), nest=nest)
        ipix = hp.ang2pix(nside, *R(theta, phi), nest=nest)
        ret = m[ipix]
    else:
        raise ValueError('Unrecognized method: {0}'.format(method))

    return ret
github simonsobs / pixell / pixell / reproject.py View on Github external
Args:
        imap: ndmap of shape (Ny,Nx)
        lmax: integer specifying maximum multipole of map
        nside: integer specifying nside of healpix map

    Returns:
        retmap: (Npix,) healpix map as array

    """
    import healpy as hp
    alm = curvedsky.map2alm(imap, lmax=lmax, spin=0)
    if alm.ndim > 1:
        assert alm.shape[0] == 1
        alm = alm[0]
    retmap = hp.alm2map(alm.astype(np.complex128), nside, lmax=lmax)
    return retmap
github nanograv / enterprise / enterprise / signals / anis_coefficients.py View on Github external
def mapFromClm_fast(clm, nside):
    """
    Given an array of C_{lm} values, produce a pixel-power-map (non-Nested) for
    healpix pixelation with nside

    @param clm:     Array of C_{lm} values (inc. 0,0 element)
    @param nside:   Nside of the healpix pixelation

    return:     Healpix pixels

    Use Healpix spherical harmonics for computational efficiency
    """
    maxl = int(np.sqrt(len(clm))) - 1
    alm = almFromClm(clm)

    h = hp.alm2map(alm, nside, maxl, verbose=False)

    return h
github HERA-Team / aipy / aipy / healpix.py View on Github external
sigma : float, scalar, optional
          The sigma of the Gaussian used to smooth the map (applied on alm)
          [in radians]
        pol : bool, optional
          If True, assumes input alms are TEB. Output will be TQU maps.
          (input must be 1 or 3 alms)
          If False, apply spin 0 harmonic transform to each alm.
          (input can be any number of alms)
          If there is only one input alm, it has no effect. Default: True.

        Returns
        -------
        maps : array or list of arrays
          A Healpix map in RING scheme at nside or a list of T,Q,U maps (if
          polarized input)'''
        return healpy.alm2map(self.get_data(), nside, 
            lmax=self._lmax, mmax=self._mmax,
            pixwin=pixwin, fwhm=fwhm, sigma=sigma, pol=pol, verbose=verbose)
    def from_map(self, data, iter=3, pol=True, use_weights=False, gal_cut=0):