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

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

interpolate: if True, bilinear interpolation using 4 nearest neighbours
        is done.

    """
    import healpy as hp
    from astropy.coordinates import SkyCoord
    import astropy.units as u
    eq_coords = ['fk5', 'j2000', 'equatorial']
    gal_coords = ['galactic']
    imap = enmap.zeros(shape, wcs)
    Ny, Nx = shape
    pixmap = enmap.pixmap(shape, wcs)
    y = pixmap[0, ...].T.ravel()
    x = pixmap[1, ...].T.ravel()
    del pixmap
    posmap = enmap.posmap(shape, wcs)
    if rot is not None:
        s1, s2 = rot.split(",")
        opos = coordinates.transform(s2,s1, posmap[::-1], pol=None)
        posmap[...] = opos[1::-1]
    th = np.rad2deg(posmap[1, ...].T.ravel())
    ph = np.rad2deg(posmap[0, ...].T.ravel())
    del posmap
    if interpolate:
        imap[y, x] = hp.get_interp_val(
            hp_map, th, ph, lonlat=True)
    else:
        ind = hp.ang2pix(hp.get_nside(hp_map),
                         th, ph, lonlat=True)
        del th
        del ph
        imap[:] = 0.

In general
		 aberrator = Aberrator(shape, wcs, ...)
		 omap = aberrator.boost(imap)
		is equivalent to
		 omap = boost_map(imap, ...)
		However, Aberrator will be more efficient when multiple maps with the
		same geometry all need to be boosted the same way, as much of the calculation
		can be precomputed in the constructor and reused for each map."""
		# Save parameters for later
		self.shape, self.wcs  = shape, wcs                                          # geometry
		self.dir,   self.beta, self.pol, self.recenter = dir, beta, pol, recenter   # boost
		self.boundary, self.order = boundary, order                                 # interpolation
		self.T0, self.freq, self.dipole = T0, freq, dipole                          # modulation
		self.map_unit, self.modulation = map_unit, modulation                       # modulation
		# Precompute displacement
		opos = enmap.posmap(shape, wcs)
		ipos, A = calc_boost(opos[::-1], dir, -beta, pol=pol, recenter=recenter)
		self.A = 1/A
		self.ipix = enmap.ndmap(enmap.sky2pix(shape, wcs, ipos[1::-1]), wcs)
		if pol:
			self.cos = np.cos(2*ipos[2])
			self.sin = np.sin(2*ipos[2])
	def aberrate(self, imap):

from one CAR geometry to another CAR geometry.
    """
    # what are the center coordinates of each geometries
    if center_source is None:
        center_source = enmap.pix2sky(
            shape_source, wcs_source,
            (shape_source[0] / 2., shape_source[1] / 2.))
    if center_target is None:
        center_target = enmap.pix2sky(
            shape_target, wcs_target,
            (shape_target[0] / 2., shape_target[1] / 2.))
    decs, ras = center_source
    dect, rat = center_target
    # what are the angle coordinates of each pixel in the target geometry
    if pos_target is None:
        pos_target = enmap.posmap(shape_target, wcs_target)
    #del pos_target
    # recenter the angle coordinates of the target from the target center
    # to the source center
    if inverse:
        transfun = lambda x: coordinates.decenter(x, (rat, dect, ras, decs))
    else:
        transfun = lambda x: coordinates.recenter(x, (rat, dect, ras, decs))
    res = coordinates.transform_meta(transfun, pos_target[1::-1], fields=["ang"])
    pix_new = enmap.sky2pix(shape_source, wcs_source, res.ocoord[1::-1])
    pix_new = np.concatenate((pix_new,res.ang[None]))
    return pix_new

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)
		if verbose: print("Computing lensed map")
		cmb_obs[...,i1:i2,:] = curvedsky.alm2map_pos(cmb_alm, raw_pos[:2], oversample=oversample, spin=spin)
		if raw_pos.shape[0] > 2 and np.any(raw_pos[2]):
			if verbose: print("Rotating polarization")
			cmb_obs[...,i1:i2,:] = enmap.rotate_pol(cmb_obs[...,i1:i2,:], raw_pos[2])
		del raw_pos

	del cmb_alm, phi_alm
	# Output in same order as specified in output argument
	res = []

# FIXME: Specifying a geometry manually is broken - see usage of r in neighborhood_pixboxes below
	# Handle arbitrary coords shape
	coords = np.asarray(coords)
	ishape = coords.shape[:-1]
	coords = coords.reshape(-1, coords.shape[-1])
	# If the output geometry was not given explicitly, then build one
	if oshape is None:
		if res is None: res = min(np.abs(imap.wcs.wcs.cdelt))*utils.degree/2
		oshape, owcs = enmap.thumbnail_geometry(r=r, res=res, proj=proj)
	# Check if we should be doing polarization rotation
	pol_compat = imap.ndim >= 3 and imap.shape[-3] == 3
	if pol is None: pol = pol_compat
	if pol and not pol_compat: raise ValueError("Polarization rotation requested, but can't interpret map shape %s as IQU map" % (str(imap.shape)))
	nsrc = len(coords)
	if verbose: print("Extracting %d %dx%d thumbnails from %s map" % (nsrc, oshape[-2], oshape[-1], str(imap.shape)))
	opos = enmap.posmap(oshape, owcs)
	# Get the pixel area around each of the coordinates
	rtot     = r + apod
	apod_pix = utils.nint(apod/(np.min(np.abs(imap.wcs.wcs.cdelt))*utils.degree))
	pixboxes = enmap.neighborhood_pixboxes(imap.shape, imap.wcs, coords, rtot)
	# Define our output maps, which we will fill below
	omaps = enmap.zeros((nsrc,)+imap.shape[:-2]+oshape, owcs, imap.dtype)
	for si, pixbox in enumerate(pixboxes):
		if oversample > 1:
			# Make the pixbox fft-friendly
			for i in range(2):
				pixbox[1,i] = pixbox[0,i] + fft.fft_len(pixbox[1,i]-pixbox[0,i], direction="above", factors=[2,3,5])
		ithumb = imap.extract_pixbox(pixbox)
		if extensive: ithumb /= ithumb.pixsizemap()
		ithumb = ithumb.apod(apod_pix, fill="median")
		if filter is not None: ithumb = filter(ithumb)
		if verbose: