How to use pixell - 10 common examples

# Get the pixel area. We assume a rectangular pixelization, so this is just
    # a function of y
    ipixsize = 4 * np.pi / (12 * nside ** 2)
    opixsize = get_pixsize_rect(shape, wcs)
    nblock = (shape[-2] + rstep - 1) // rstep
    for bi in range(comm.rank, nblock, comm.size):
        if bi % comm.size != comm.rank:
            continue
        i = bi * rstep
        rdec = dec[i : i + rstep]
        opos = np.zeros((2, len(rdec), len(ra)))
        opos[0] = rdec[:, None]
        opos[1] = ra[None, :]
        if rot:
            # This is unreasonably slow
            ipos = coordinates.transform("equ", "gal", opos[::-1], pol=True)
        else:
            ipos = opos[::-1]
        pix[i : i + rstep, :] = hp.ang2pix(nside, np.pi / 2 - ipos[1], ipos[0])
        del ipos, opos
    for i in range(0, shape[-2], rstep):
        pix[i : i + rstep] = utils.allreduce(pix[i : i + rstep], comm)
    omap = enmap.zeros((1,) + shape, wcs, dtype)
    imap = np.array(hmap).astype(dtype)
    imap = imap[None]
    if do_mask:
        bad = hp.mask_bad(imap)
        bad |= imap <= 0
        imap[bad] = 0
        del bad
    # Read off the nearest neighbor values
    omap[:] = imap[:, pix]

def test_pospix(self):
        # Posmap separable and non-separable on CAR
        for res in [6,12,24]:
            shape,wcs = enmap.fullsky_geometry(res=np.deg2rad(res/60.),proj='car')
            posmap1 = enmap.posmap(shape,wcs)
            posmap2 = enmap.posmap(shape,wcs,separable=True)
            assert np.all(np.isclose(posmap1,posmap2))

        # Pixmap plain
        pres = 0.5
        shape,wcs = enmap.geometry(pos=(0,0),shape=(30,30),res=pres*u.degree,proj='plain')
        yp,xp = enmap.pixshapemap(shape,wcs)
        assert np.all(np.isclose(yp,pres*u.degree))
        assert np.all(np.isclose(xp,pres*u.degree))
        yp,xp = enmap.pixshape(shape,wcs)
        parea = enmap.pixsize(shape,wcs)
        assert np.isclose(parea,(pres*u.degree)**2)
        assert np.isclose(yp,pres*u.degree)
        assert np.isclose(xp,pres*u.degree)
        pmap = enmap.pixsizemap(shape,wcs)
        assert np.all(np.isclose(pmap,(pres*u.degree)**2))

import sys
sys.path.append('../../tests')
import test_pixell as ptests
import os
import numpy as np
from pixell import enmap

version = sys.argv[1]

obs_pos,grad,raw_pos = ptests.get_offset_result(1.,np.float64)
enmap.write_map("MM_offset_obs_pos_%s.fits" % version,obs_pos)
enmap.write_map("MM_offset_grad_%s.fits"  % version,grad)
enmap.write_map("MM_offset_raw_pos_%s.fits"  % version,raw_pos)

lensed,unlensed = ptests.get_lens_result(1.,400,np.float64)
enmap.write_map("MM_lensed_%s.fits"  % version,lensed)
enmap.write_map("MM_unlensed_%s.fits"  % version,unlensed)

# Flat-sky
        px = 2.0
        N = 300
        shape,iwcs = enmap.geometry(pos=(0,0),res=np.deg2rad(px/60.),shape=(300,300))
        shape = (3,) + shape
        a = enmap.zeros(shape,iwcs)
        a = a[...,::-1]
        wcs = a.wcs

        seed = 100
        imap = enmap.rand_map(shape,wcs,ps_cmb,seed=seed)
        kmap = enmap.map2harm(imap.copy())
        rmap = enmap.harm2map(kmap,spin=0) # reference map

        imap = imap[...,::-1]
        kmap = enmap.map2harm(imap.copy())
        rmap2 = enmap.harm2map(kmap,spin=0)[...,::-1] # comparison map
        
        assert np.all(np.isclose(rmap[0],rmap2[0]))
        assert np.all(np.isclose(rmap[1],rmap2[1],atol=1e0))
        assert np.all(np.isclose(rmap[2],rmap2[2],atol=1e0))

box = np.deg2rad(np.array(e['box_deg']))
                    cutout = enmap.read_map(filename,box=box)
                    cutout_internal = imap.submap(box=box)
                elif e['type']=='postage':
                    dec_deg,ra_deg = e['center_deg']
                    width_arcmin = e['width_arcmin']
                    res_arcmin = e['res_arcmin']
                    file_proxy = enmap.read_map(filename,delayed=True)
                    coords = np.array([dec_deg,ra_deg]) * utils.degree
                    cutout = reproject.thumbnails(file_proxy,coords,
                                                  r=width_arcmin/2.*utils.arcmin,
                                                  res=res_arcmin*utils.arcmin,
                                                  proj='tan')
                    cutout_internal = reproject.thumbnails(imap,coords,
                                                           r = width_arcmin/2.*utils.arcmin,
                                                           res=res_arcmin*utils.arcmin,
                                                           proj='tan')
                check_equality(cutout,cutout_internal)
                pixels = get_reference_pixels(cutout.shape)
                results[g][s]['refpixels'] = get_pixel_values(cutout,pixels)
                results[g][s]['meansquare'] = get_meansquare(cutout)

    os.remove(filename)
    return results,config['result_name']

imap2 = imap.copy()[...,::-1]
        oalm = curvedsky.map2alm(imap2.copy(),lmax=lmax)
        rmap2 = curvedsky.alm2map(oalm,enmap.empty(shape,wcs),spin=0)

        assert np.all(np.isclose(rmap[0],rmap2[0]))
        assert np.all(np.isclose(rmap[1],rmap2[1]))
        assert np.all(np.isclose(rmap[2],rmap2[2]))
        

        # Flat-sky
        px = 2.0
        N = 300
        shape,iwcs = enmap.geometry(pos=(0,0),res=np.deg2rad(px/60.),shape=(300,300))
        shape = (3,) + shape
        a = enmap.zeros(shape,iwcs)
        a = a[...,::-1]
        wcs = a.wcs

        seed = 100
        imap = enmap.rand_map(shape,wcs,ps_cmb,seed=seed)
        kmap = enmap.map2harm(imap.copy())
        rmap = enmap.harm2map(kmap,spin=0) # reference map

        imap = imap[...,::-1]
        kmap = enmap.map2harm(imap.copy())
        rmap2 = enmap.harm2map(kmap,spin=0)[...,::-1] # comparison map
        
        assert np.all(np.isclose(rmap[0],rmap2[0]))
        assert np.all(np.isclose(rmap[1],rmap2[1],atol=1e0))
        assert np.all(np.isclose(rmap[2],rmap2[2],atol=1e0))

We compare these maps.
        """
        ells,cltt,clee,clbb,clte = np.loadtxt(DATA_PREFIX+"cosmo2017_10K_acc3_lensedCls.dat",unpack=True)
        ps_cmb = np.zeros((3,3,ells.size))
        ps_cmb[0,0] = cltt
        ps_cmb[1,1] = clee
        ps_cmb[2,2] = clbb
        ps_cmb[1,0] = clte
        ps_cmb[0,1] = clte
        np.random.seed(100)

        # Curved-sky is fine
        lmax = 1000
        alm = curvedsky.rand_alm_healpy(ps_cmb,lmax=lmax)
        shape,iwcs = enmap.fullsky_geometry(res=np.deg2rad(10./60.))
        wcs = enmap.empty(shape,iwcs)[...,::-1].wcs
        shape = (3,) + shape
        imap = curvedsky.alm2map(alm,enmap.empty(shape,wcs))
        oalm = curvedsky.map2alm(imap.copy(),lmax=lmax)
        rmap = curvedsky.alm2map(oalm,enmap.empty(shape,wcs),spin=0)

        imap2 = imap.copy()[...,::-1]
        oalm = curvedsky.map2alm(imap2.copy(),lmax=lmax)
        rmap2 = curvedsky.alm2map(oalm,enmap.empty(shape,wcs),spin=0)

        assert np.all(np.isclose(rmap[0],rmap2[0]))
        assert np.all(np.isclose(rmap[1],rmap2[1]))
        assert np.all(np.isclose(rmap[2],rmap2[2]))
        

        # Flat-sky
        px = 2.0

ps_cmb[0,0] = cltt
        ps_cmb[1,1] = clee
        ps_cmb[2,2] = clbb
        ps_cmb[1,0] = clte
        ps_cmb[0,1] = clte
        np.random.seed(100)

        # Curved-sky is fine
        lmax = 1000
        alm = curvedsky.rand_alm_healpy(ps_cmb,lmax=lmax)
        shape,iwcs = enmap.fullsky_geometry(res=np.deg2rad(10./60.))
        wcs = enmap.empty(shape,iwcs)[...,::-1].wcs
        shape = (3,) + shape
        imap = curvedsky.alm2map(alm,enmap.empty(shape,wcs))
        oalm = curvedsky.map2alm(imap.copy(),lmax=lmax)
        rmap = curvedsky.alm2map(oalm,enmap.empty(shape,wcs),spin=0)

        imap2 = imap.copy()[...,::-1]
        oalm = curvedsky.map2alm(imap2.copy(),lmax=lmax)
        rmap2 = curvedsky.alm2map(oalm,enmap.empty(shape,wcs),spin=0)

        assert np.all(np.isclose(rmap[0],rmap2[0]))
        assert np.all(np.isclose(rmap[1],rmap2[1]))
        assert np.all(np.isclose(rmap[2],rmap2[2]))
        

        # Flat-sky
        px = 2.0
        N = 300
        shape,iwcs = enmap.geometry(pos=(0,0),res=np.deg2rad(px/60.),shape=(300,300))
        shape = (3,) + shape
        a = enmap.zeros(shape,iwcs)

parea = enmap.pixsize(shape,wcs)
        assert np.isclose(parea,(pres*u.degree)**2)
        assert np.isclose(yp,pres*u.degree)
        assert np.isclose(xp,pres*u.degree)
        pmap = enmap.pixsizemap(shape,wcs)
        assert np.all(np.isclose(pmap,(pres*u.degree)**2))

        # Pixmap CAR
        pres = 0.1
        dec_cut = 89.5 # pixsizemap is not accurate near the poles currently
        shape,wcs = enmap.band_geometry(dec_cut=dec_cut*u.degree,res=pres*u.degree,proj='car')
        # Current slow and general but inaccurate near the poles implementation
        pmap = enmap.pixsizemap(shape,wcs)
        # Fast CAR-specific pixsizemap implementation
        dra, ddec = wcs.wcs.cdelt*u.degree
        dec = enmap.posmap([shape[-2],1],wcs)[0,:,0]
        area = np.abs(dra*(np.sin(np.minimum(np.pi/2.,dec+ddec/2))-np.sin(np.maximum(-np.pi/2.,dec-ddec/2))))
        Nx = shape[-1]
        pmap2 = enmap.ndmap(area[...,None].repeat(Nx,axis=-1),wcs)
        assert np.all(np.isclose(pmap,pmap2))

def get_offset_result(res=1.,dtype=np.float64,seed=1):
    shape,wcs  = enmap.fullsky_geometry(res=np.deg2rad(res))
    shape = (3,) + shape
    obs_pos = enmap.posmap(shape, wcs)
    np.random.seed(seed)
    grad = enmap.enmap(np.random.random(shape),wcs)*1e-3
    raw_pos = enmap.samewcs(lensing.offset_by_grad(obs_pos, grad, pol=shape[-3]>1, geodesic=True), obs_pos)
    return obs_pos,grad,raw_pos