How to use healpy - 10 common examples

# number of subpixels in each superpixel
    npix = np.array((nside_subpix // nside_superpix) ** 2, ndmin=1)
    x = np.arange(np.max(npix), dtype=int)
    idx = idx * npix

    if not np.all(npix[0] == npix):
        x = np.broadcast_to(x, idx.shape + x.shape)
        idx = idx[..., None] + x
        idx[x >= np.broadcast_to(npix[..., None], x.shape)] = INVALID_INDEX.int
    else:
        idx = idx[..., None] + x

    if not nest:
        m = idx == INVALID_INDEX.int
        idx[m] = 0
        idx = hp.nest2ring(nside_subpix[..., None], idx)
        idx[m] = INVALID_INDEX.int

    return idx

if npix > 12:
        # Determine which pixels comprise multi-pixel tiles.
        ipix = np.flatnonzero(
            (m[0::4] == m[1::4]) &
            (m[0::4] == m[2::4]) &
            (m[0::4] == m[3::4]))

        if len(ipix):
            ipix = (4 * ipix +
                np.expand_dims(np.arange(4, dtype=np.intp), 1)).T.ravel()

            nside = hp.npix2nside(npix)

            # Downsample.
            m_lores = hp.ud_grade(
                m, nside // 2, order_in='NESTED', order_out='NESTED')

            # Interpolate recursively.
            _interpolate_level(m_lores)

            # Record interpolated multi-pixel tiles.
            m[ipix] = hp.get_interp_val(
                m_lores, *hp.pix2ang(nside, ipix, nest=True), nest=True)

import healpy as hp
import matplotlib.pyplot as plt
import pymaster as nmt

# This script showcases the apodization routine included with pymaster
# and the three apodization modes supported.


# Read input binary mask
mask_raw = hp.read_map("mask.fits", verbose=False)

# The following function calls create apodized versions of the raw mask
# with an apodization scale of 2.5 degrees using three different methods

# Apodization scale in degrees
aposcale = 2.5

# C1 and C2: in these cases, pixels are multiplied by a factor f
#            (with 0<=f<=1) based on their distance to the nearest fully
#            masked pixel. The choices of f in each case are documented in
#            Section 3.4 of the C API documentation. All pixels separated
#            from any masked pixel by more than the apodization scale are
#            left untouched.
mask_C1 = nmt.mask_apodization(mask_raw, aposcale, apotype="C1")
mask_C2 = nmt.mask_apodization(mask_raw, aposcale, apotype="C2")

def setUp(self) :
        self.nside=64
        self.lmax=3*self.nside-1
        self.ntemp=5
        self.npix=int(hp.nside2npix(self.nside))
        self.msk=np.ones(self.npix)
        self.mps=np.zeros([3,self.npix])
        self.tmp=np.zeros([self.ntemp,3,self.npix])
        self.beam=np.ones(self.lmax+1)

        th,ph=hp.pix2ang(self.nside,np.arange(self.npix))
        sth=np.sin(th); cth=np.cos(th)
        self.mps[0]=np.sqrt(15./2./np.pi)*sth**2*np.cos(2*ph) #Re(Y_22)
        self.mps[1]=-np.sqrt(15./2./np.pi)*sth**2/4.          #_2Y^E_20 + _2Y^B_30
        self.mps[2]=-np.sqrt(105./2./np.pi)*cth*sth**2/2.
        for i in range(self.ntemp) :
            self.tmp[i][0]=np.sqrt(15./2./np.pi)*sth**2*np.cos(2*ph) #Re(Y_22)
            self.tmp[i][1]=-np.sqrt(15./2./np.pi)*sth**2/4.          #_2Y^E_20 + _2Y^B_30
            self.tmp[i][2]=-np.sqrt(105./2./np.pi)*cth*sth**2/2.

def setUp(self) :
        self.nside=256
        self.th0=np.pi/4
        self.msk=np.zeros(hp.nside2npix(self.nside),dtype=float)
        self.th,ph=hp.pix2ang(self.nside,np.arange(hp.nside2npix(self.nside),dtype=int))
        self.msk[self.th

print 'ini: mono=%s, dipole=%s'%(c,d)
mono,dipole = H.pixelfunc.fit_dipole(sig)
print 'fit: mono=%s, dipole=%s'%(mono,dipole)

resfit = dot(vec.T,dipole)+mono

diff = sig.copy()
diff[sig != H.UNSEEN] -= resfit[sig != H.UNSEEN]

H.mollview(diff)

# use remove_dipole
m2 = H.pixelfunc.remove_dipole(sig)
m3 = H.pixelfunc.remove_monopole(sig)

H.mollview(m2)
H.mollview(m3)

thb,phb = H.pix2ang(nside,arange(npix))
glat=90.-thb*180/pi

w=abs(glat)<10
sig[w] += 3.*randn(w.sum())+10

H.mollview(sig)

print 'ini: mono=%s, dipole=%s'%(c,d)
mono,dipole = H.pixelfunc.fit_dipole(sig)
print 'fit: mono=%s, dipole=%s'%(mono,dipole)

resfit = dot(vec.T,dipole)+mono

diff = sig.copy()
diff[sig != H.UNSEEN] -= resfit[sig != H.UNSEEN]

H.mollview(diff)

# use remove_dipole
m2 = H.pixelfunc.remove_dipole(sig)
m3 = H.pixelfunc.remove_monopole(sig)

H.mollview(m2)
H.mollview(m3)

#  GNU General Public License for more details.
# 
#  You should have received a copy of the GNU General Public License
#  along with Healpy; if not, write to the Free Software
#  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
# 
#  For more information about Healpy, see http://code.google.com/p/healpy
# 
import healpy as H
from numpy import *
from pylab import *

nside = 32
npix = H.nside2npix(nside)

vec = array(H.pix2vec(nside,arange(npix)))
d = 3*array([0.5,0.3,0.2])
#d /= norm(d)
c = 2.

sig = dot(vec.T,d)+c
sig += randn(npix)
sig[randint(sig.size,size=10000)] = H.UNSEEN

thb,phb = H.pix2ang(nside,arange(npix))
glat=90.-thb*180/pi

w=abs(glat)<10
sig[w] += 3.*randn(w.sum())+10

H.mollview(sig)

H.mollview(sig)

print 'ini: mono=%s, dipole=%s'%(c,d)
mono,dipole = H.pixelfunc.fit_dipole(sig)
print 'fit: mono=%s, dipole=%s'%(mono,dipole)

resfit = dot(vec.T,dipole)+mono

diff = sig.copy()
diff[sig != H.UNSEEN] -= resfit[sig != H.UNSEEN]

H.mollview(diff)

# use remove_dipole
m2 = H.pixelfunc.remove_dipole(sig)
m3 = H.pixelfunc.remove_monopole(sig)

H.mollview(m2)
H.mollview(m3)

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))