How to use the pixell.enmap function in pixell

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

def cutout(imap, width=None, ra=None, dec=None, pad=1, corner=False,
           res=None, npix=None, return_slice=False,sindex=None):
    if type(imap) == str:
        shape, wcs = enmap.read_map_geometry(imap)
    else:
        shape, wcs = imap.shape, imap.wcs
    Ny, Nx = shape[-2:]
    def fround(x):
        return int(np.round(x))
    iy, ix = enmap.sky2pix(shape, wcs, coords=(dec, ra), corner=corner)
    if res is None:
        res = np.min(enmap.extent(shape, wcs) / shape[-2:])
    if npix is None:
        npix = int(width / res)
    if fround(iy - npix / 2) < pad or fround(ix - npix / 2) < pad or \
       fround(iy + npix / 2) > (Ny - pad) or \
       fround(ix + npix / 2) > (Nx - pad):
        return None
    if sindex is None:
        s = np.s_[...,fround(iy - npix / 2. + 0.5):fround(iy + npix / 2. + 0.5),

def prepare_map_field(map, args, crange=None, printer=noprint):
	if crange is None:
		with printer.time("ranges", 3):
			crange = get_color_range(map, args)
	if map.ndim == 2:
		map = map[None]
	if args.autocrop_each:
		map = enmap.autocrop(map)
	with printer.time("colorize", 3):
		color = map_to_color(map, crange, args)
	return map, color

def rotate_map(imap, shape_target=None, wcs_target=None, shape_source=None,
               wcs_source=None, pix_target=None, **kwargs):
    if pix_target is None:
        pix_target = get_rotated_pixels(
            shape_source, wcs_source, shape_target, wcs_target)
    else:
        assert (shape_target is None) and (
            wcs_target is None), "Both pix_target and shape_target, \
            wcs_target must not be specified."
    rotmap = enmap.at(imap, pix_target[:2], unit="pix", **kwargs)
    return rotmap

bsize = int(shape[-2]/(nblock+0.5))
	# Allocate output maps
	if "p" in output: phi_map   = enmap.empty(shape[-2:], wcs, dtype=dtype)
	if "k" in output:
		kappa_map = enmap.empty(shape[-2:], wcs, dtype=dtype)
		kappa_alm = phi_to_kappa(phi_alm,phi_ainfo=phi_ainfo)
		for i1 in range(0, shape[-2], bsize):
			curvedsky.alm2map(kappa_alm, kappa_map[...,i1:i1+bsize,:])
		del kappa_alm
	if "a" in output: grad_map  = enmap.empty((2,)+shape[-2:], wcs, dtype=dtype)
	if "u" in output: cmb_raw   = enmap.empty(shape, wcs, dtype=dtype)
	if "l" in output: cmb_obs   = enmap.empty(shape, wcs, dtype=dtype)
	# Then loop over dec bands
	for i1 in range(0, shape[-2], bsize):
		i2 = min(i1+bsize, shape[-2])
		lshape, lwcs = enmap.slice_geometry(shape, wcs, (slice(i1,i2),slice(None)))
		if "p" in output:
			if verbose: print("Computing phi map")
			curvedsky.alm2map(phi_alm, phi_map[...,i1:i2,:])
		if verbose: print("Computing grad map")
		if "a" in output: grad = grad_map[...,i1:i2,:]
		else: grad = enmap.zeros((2,)+lshape[-2:], lwcs, dtype=dtype)
		curvedsky.alm2map(phi_alm, grad, deriv=True)
		if "l" not in output: continue
		if verbose: print("Computing observed coordinates")
		obs_pos = enmap.posmap(lshape, lwcs)
		if verbose: print("Computing alpha map")
		raw_pos = enmap.samewcs(offset_by_grad(obs_pos, grad, pol=shape[-3]>1, geodesic=geodesic), obs_pos)
		del obs_pos, grad
		if "u" in output:
			if verbose: print("Computing unlensed map")
			curvedsky.alm2map(cmb_alm, cmb_raw[...,i1:i2,:], spin=spin)