How to use the pixell.utils function in pixell

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']

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

def hor2cel(coord, time, site, copy=True):
	from enlib import pyfsla
	from enlib import iers
	coord  = np.array(coord, copy=copy)
	trepr  = time[len(time)/2]
	info   = iers.lookup(trepr)
	ao = pyfsla.sla_aoppa(trepr, info.dUT, site.lon*utils.degree, site.lat*utils.degree, site.alt,
		info.pmx*utils.arcsec, info.pmy*utils.arcsec, site.T, site.P, site.hum,
		299792.458/site.freq, site.lapse)
	am = pyfsla.sla_mappa(2000.0, trepr)
	# This involves a transpose operation, which is not optimal
	pyfsla.aomulti(time, coord.T, ao, am)
	return coord

def to_healpix(imap, omap=None, nside=0, order=3, chunk=100000, destroy_input=False):
	"""Project the enmap "imap" onto the healpix pixelization. If omap is given,
	the output will be written to it. Otherwise, a new healpix map will be constructed.
	The healpix map must be in RING order. nside controls the resolution of the output map.
	If 0, nside is chosen such that the output map is higher resolution than the input.
	This is needed to avoid losing information. To go to a lower-resolution output map,
	you should first degrade the input map. The chunk argument affects the speed/memory
	tradeoff of the function. Higher values use more memory, and might (and might not)
	give higher speed. If destroy_input is True, then the input map will be prefiltered
	in-place, which saves memory but modifies its values."""
	import healpy
	if not destroy_input and order > 1: imap = imap.copy()
	if order > 1:
		imap = utils.interpol_prefilter(imap, order=order, inplace=True)
	if omap is None:
		# Generate an output map
		if not nside:
			npix_full_cyl = 4*np.pi/imap.pixsize()
			nside = 2**int(np.floor(np.log2((npix_full_cyl/12)**0.5)))
		npix = 12*nside**2
		omap = np.zeros(imap.shape[:-2]+(npix,),imap.dtype)
	else:
		nside = healpy.npix2nside(omap.shape[-1])
	npix = omap.shape[-1]
	# Interpolate values at output pixel positions
	for i in range(0, npix, chunk):
		pos   = np.array(healpy.pix2ang(nside, np.arange(i, min(npix,i+chunk))))
		# Healpix uses polar angle, not dec
		pos[0] = np.pi/2 - pos[0]
		omap[...,i:i+chunk] = imap.at(pos, order=order, mask_nan=False, prefilter=False)

def multi_pow(mat, exp, axes=[0,1]):
	"""Raise each sub-matrix of mat (ncomp,ncomp,...) to
	the given exponent in eigen-space."""
	return samewcs(utils.eigpow(mat, exp, axes=axes), mat)

"""Render a map field using matplotlib. Less tested and
	maintained than draw_map_field, and supports fewer features.
	Returns an object one can call savefig on to draw."""
	map, color = prepare_map_field(map, args, crange, printer=printer)
	# Set up matplotlib. We do it locally here to
	# avoid having it as a dependency in general
	with printer.time("matplotplib", 3):
		import matplotlib
		matplotlib.use("Agg")
		from matplotlib import pyplot, ticker
		matplotlib.rcParams.update({'font.size': 10})
		dpi, pad = args.mpl_dpi, args.mpl_pad
		winch, hinch = map.shape[2]/dpi, map.shape[1]/dpi
		fig  = pyplot.figure(figsize=(winch+pad,hinch+pad))
		box  = map.box()*180/np.pi
		pyplot.imshow(utils.moveaxis(color,0,2), extent=[box[0,1],box[1,1],box[1,0],box[0,0]])
		# Make conformal in center of image
		pyplot.axes().set_aspect(1/np.cos(np.mean(map.box()[:,0])))
		if args.grid % 2:
			ax = pyplot.axes()
			ticks = np.full(2,1.0)
			ticks[:] = parse_list(args.ticks)
			ax.xaxis.set_major_locator(ticker.MultipleLocator(ticks[1]))
			ax.yaxis.set_major_locator(ticker.MultipleLocator(ticks[0]))
			if args.subticks:
				ax.xaxis.set_minor_locator(ticker.MultipleLocator(args.sub))
				ax.yaxis.set_minor_locator(ticker.MultipleLocator(args.sub))
				pyplot.minorticks_on()
				pyplot.grid(True, which="major", linewidth=2)
				pyplot.grid(True, which="minor", linewidth=1)
			else:
				pyplot.grid(True)

def tele2hor(coord, site, copy=True):
	coord = np.array(coord, copy=copy)
	coord = euler_rot([site.base_az*utils.degree, site.base_tilt*utils.degree, -site.base_az*utils.degree], coord)
	return coord

def read_fits_geometry(fname, hdu=None):
	"""Read an enmap wcs from the specified fits file. By default,
	the map and coordinate system will be read from HDU 0. Use
	the hdu argument to change this. The map must be stored as
	a fits image."""
	if hdu is None: hdu = 0
	with utils.nowarn():
		hdu = astropy.io.fits.open(fname)[hdu]
	if hdu.header["NAXIS"] < 2:
		raise ValueError("%s is not an enmap (only %d axes)" % (fname, hdu.header["NAXIS"]))
	with warnings.catch_warnings():
		wcs = wcsutils.WCS(hdu.header).sub(2)
	shape = tuple([hdu.header["NAXIS%d"%(i+1)] for i in range(hdu.header["NAXIS"])[::-1]])
	return shape, wcs