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

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

cylindrical grid and then interpolates to the actual pixels. This is the
	general way of doing things, but not the fastest. Computing pos and
	interpolating takes a significant amount of time."""
	alm_full = np.atleast_2d(alm)
	if ainfo is None: ainfo = sharp.alm_info(nalm=alm_full.shape[-1])
	ashape, ncomp = alm_full.shape[:-2], alm_full.shape[-2]
	if deriv:
		# If we're computing derivatives, spin isn't allowed.
		# alm must be either [ntrans,nelem] or [nelem],
		# and the output will be [ntrans,2,ny,nx] or [2,ny,nx]
		ashape = ashape + (ncomp,)
		ncomp = 2
	tmap   = make_projectable_map_by_pos(pos, ainfo.lmax, ashape+(ncomp,), oversample, alm.real.dtype)
	alm2map_cyl(alm, tmap, ainfo=ainfo, spin=spin, deriv=deriv, direct=True, verbose=verbose)
	# Project down on our final pixels. This will result in a slight smoothing
	res = enmap.samewcs(tmap.at(pos[:2], mode="wrap"), pos)
	# Remove any extra dimensions we added
	if alm.ndim == alm_full.ndim-1: res = res[0]
	return res

def lens_map_flat(cmb_map, phi_map):
	raw_pix  = cmb_map.pixmap() + enmap.grad_pix(phi_map)
	# And extract the interpolated values. Because of a bug in map_pixels with
	# mode="wrap", we must handle wrapping ourselves.
	npad = int(np.ceil(max(np.max(-raw_pix),np.max(raw_pix-np.array(cmb_map.shape[-2:])[:,None,None]))))
	pmap = enmap.pad(cmb_map, npad, wrap=True)
	return enmap.samewcs(utils.interpol(pmap, raw_pix+npad, order=4, mode="wrap"), cmb_map)

def aberrate(self, imap):
		"""Apply the aberration part of the doppler boost to the map imap"""
		omap = enmap.samewcs(imap.at(self.ipix, unit="pix", mode=self.boundary, order=self.order), imap)
		if self.pol and imap.ndim > 2:
			omap1 = omap[...,1,:,:].copy()
			omap[...,1,:,:] =  self.cos*omap1 + self.sin*omap[...,2,:,:]
			omap[...,2,:,:] = -self.sin*omap1 + self.cos*omap[...,2,:,:]
		return omap
	def modulate(self, imap):

# 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 = []
	for c in output:
		if   c == "l": res.append(cmb_obs.reshape(oshape))