How to use the pixell.enmap.empty function in pixell

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)

# would require us to copy out multiple slices.
	nslice = (nx+nphi-1)//nphi
	islice, oslice = [], []
	def negnone(x): return x if x >= 0 else None
	for i in range(nslice):
		# i1:i2 is the range of pixels in the original map to use
		i1, i2 = i*nphi, min((i+1)*nphi,nx)
		islice.append((Ellipsis, slice(i1,i2)))
		# yslice and xslice give the range of pixels in our temporary map to use.
		# This is 0:(i2-i1) if we're not flipping, but if we flip we count from
		# the opposite direction: nx-1:nx-1-(i2-i1):-1
		yslice = slice(-1,None,-1)  if flipy else slice(None)
		xslice = slice(nx-1,negnone(nx-1-(i2-i1)),-1) if flipx else slice(0,i2-i1)
		oslice.append((Ellipsis,yslice,xslice))
	if verbose: print("Allocating shape %s dtype %s intermediate map" % (str(oshape),np.dtype(map.dtype).char))
	return enmap.empty(oshape, owcs, dtype=map.dtype), islice, oslice

if delta_theta is None: bsize = shape[-2]
	else:
		bsize = utils.nint(abs(delta_theta/utils.degree/wcs.wcs.cdelt[1]))
		# Adjust bsize so we don't get any tiny blocks at the end
		nblock= shape[-2]//bsize
		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")

def displace_map(imap, pix, order=3, mode="spline", border="cyclic", trans=False, deriv=False):
	"""Displace map m[{pre},ny,nx] by pix[2,ny,nx], where pix indicates the location
	in the input map each output pixel should get its value from (float). The output
	is [{pre},ny,nx]."""
	if not deriv: omap = imap.copy()
	else:         omap = enmap.empty((2,)+imap.shape, imap.wcs, imap.dtype)
	if not trans:
		interpol.map_coordinates(imap, pix, omap, order=order, mode=mode, border=border, trans=trans, deriv=deriv)
	else:
		interpol.map_coordinates(omap, pix, imap, order=order, mode=mode, border=border, trans=trans, deriv=deriv)
	return omap

def lens_map_curved(shape, wcs, phi_alm, cmb_alm, phi_ainfo=None, maplmax=None, dtype=np.float64, oversample=2.0, spin=[0,2], output="l", geodesic=True, verbose=False, delta_theta=None):
	from . import curvedsky, sharp
	# Restrict to target number of components
	oshape  = shape[-3:]
	if len(oshape) == 2: shape = (1,)+tuple(shape)
	if delta_theta is None: bsize = shape[-2]
	else:
		bsize = utils.nint(abs(delta_theta/utils.degree/wcs.wcs.cdelt[1]))
		# Adjust bsize so we don't get any tiny blocks at the end
		nblock= shape[-2]//bsize
		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,:])