How to use the pixell.wcsutils.WCS function in pixell

def fullsky_geometry(res=None, shape=None, dims=(), proj="car"):
	"""Build an enmap covering the full sky, with the outermost pixel centers
	at the poles and wrap-around points. Assumes a CAR (clenshaw curtis variant)
	projection for now."""
	assert proj == "car", "Only CAR fullsky geometry implemented"
	if shape is None:
	 res   = np.zeros(2)+res
	 shape = utils.nint(([1*np.pi,2*np.pi]/res) + (1,0))
	ny, nx = shape
	assert abs(res[0] * (ny-1) - np.pi) < 1e-8, "Vertical resolution does not evenly divide the sky; this is required for SHTs."
	assert abs(res[1] * nx - 2*np.pi)   < 1e-8, "Horizontal resolution does not evenly divide the sky; this is required for SHTs."
	wcs   = wcsutils.WCS(naxis=2)
	# Note the reference point is shifted by half a pixel to keep
	# the grid in bounds, from ra=180+cdelt/2 to ra=-180+cdelt/2.
	wcs.wcs.crval = [res[0]/2/utils.degree,0]
	wcs.wcs.cdelt = [-360./nx,180./(ny-1)]
	wcs.wcs.crpix = [nx//2+0.5,ny//2+1]
	wcs.wcs.ctype = ["RA---CAR","DEC--CAR"]
	return dims+(ny,nx), wcs

"""Read an enmap from the specified hdf file. Two formats
	are supported. The old enmap format, which simply used
	a bounding box to specify the coordinates, and the new
	format, which uses WCS properties. The latter is used if
	available. With the old format, plate carree projection
	is assumed. Note: some of the old files have a slightly
	buggy wcs, which can result in 1-pixel errors."""
	import h5py
	with h5py.File(fname,"r") as hfile:
		data = hfile["data"]
		hwcs = hfile["wcs"]
		header = astropy.io.fits.Header()
		for key in hwcs:
			header[key] = fix_python3(hwcs[key].value)
		if wcs is None:
			wcs = wcsutils.WCS(header).sub(2)
		proxy = ndmap_proxy_hdf(data, wcs, fname=fname, threshold=sel_threshold)
		return read_helper(proxy, sel=sel, box=box, pixbox=pixbox, geometry=geometry, wrap=wrap, mode=mode, delayed=delayed)

decrange = (decrange-np.mean(decrange))*1.1+np.mean(decrange)
	decrange = np.array([max(-np.pi/2,decrange[0]),min(np.pi/2,decrange[1])])
	decrange /= utils.degree
	wdec = np.abs(decrange[1]-decrange[0])
	# The shortest wavelength in the alm is about 2pi/lmax. We need at least
	# two samples per mode.
	res = 180./lmax/oversample
	# Set up an intermediate coordinate system for the SHT. We will use
	# CAR coordinates conformal on the equator, with a pixel on each pole.
	# This will give it clenshaw curtis pixelization.
	nx    = utils.nint(360/res)
	nytot = utils.nint(180/res)
	# First set up the pixelization for the whole sky. Negative cdelt to
	# make sharp extra happy. Not really necessary, but makes some things
	# simpler later.
	wcs   = wcsutils.WCS(naxis=2)
	wcs.wcs.ctype = ["RA---CAR","DEC--CAR"]
	wcs.wcs.crval = [ra_ref,0]
	wcs.wcs.cdelt = [360./nx,-180./nytot]
	wcs.wcs.crpix = [nx/2.0+1,nytot/2.0+1]
	# Then find the subset that includes the dec range we want
	y1= utils.nint(wcs.wcs_world2pix(0,decrange[0],0)[1])
	y2= utils.nint(wcs.wcs_world2pix(0,decrange[1],0)[1])
	y1, y2 = min(y1,y2), max(y1,y2)
	# Offset wcs to start at our target range
	ny = y2-y1
	wcs.wcs.crpix[1] -= y1
	# Construct the map. +1 to put extra pixel at pole when we are fullsky
	if verbose: print("Allocating shape %s dtype %s intermediate map" % (dims+(ny+1,nx),np.dtype(dtype).char))
	tmap = enmap.zeros(dims+(ny+1,nx),wcs,dtype=dtype)
	return tmap

def read_fits(fname, hdu=None, sel=None, box=None, pixbox=None, geometry=None, wrap="auto", mode=None, sel_threshold=10e6, wcs=None, delayed=False, verbose=False):
	"""Read an enmap 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 sel is specified, it should be a slice
	that will be applied to the image before reading. This avoids
	reading more of the image than necessary. Instead of sel,
	a coordinate box [[yfrom,xfrom],[yto,xto]] can be specified."""
	if hdu is None: hdu = 0
	hdu = astropy.io.fits.open(fname)[hdu]
	ndim = len(hdu.shape)
	if hdu.header["NAXIS"] < 2:
		raise ValueError("%s is not an enmap (only %d axes)" % (fname, hdu.header["NAXIS"]))
	if wcs is None:
		with warnings.catch_warnings():
			wcs = wcsutils.WCS(hdu.header).sub(2)
	proxy = ndmap_proxy_fits(hdu, wcs, fname=fname, threshold=sel_threshold, verbose=verbose)
	return read_helper(proxy, sel=sel, box=box, pixbox=pixbox, geometry=geometry, wrap=wrap, mode=mode, delayed=delayed)

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