How to use the pixell.curvedsky.ShapeError function in pixell

To help you get started, we’ve selected a few pixell examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github simonsobs / pixell / pixell / curvedsky.py View on Github external
# First check if we need flipping. Sharp wants theta,phi increasing,
	# which means dec decreasing and ra increasing.
	flipx = map.wcs.wcs.cdelt[0] < 0
	flipy = map.wcs.wcs.cdelt[1] > 0
	if flipx: map = map[...,:,::-1]
	if flipy: map = map[...,::-1,:]
	# Then check if the map satisfies the lat-ring requirements
	ny, nx = map.shape[-2:]
	vy,vx = enmap.pix2sky(map.shape, map.wcs, [np.arange(ny),np.zeros(ny)])
	hy,hx = enmap.pix2sky(map.shape, map.wcs, [np.zeros(nx),np.arange(nx)])
	dx    = hx[1:]-hx[:-1]
	dx    = dx[np.isfinite(dx)] # Handle overextended coordinates

	if not np.allclose(dx,dx[0]): raise ShapeError("Map must have constant phi spacing")
	nphi = utils.nint(2*np.pi/dx[0])
	if not np.allclose(2*np.pi/nphi,dx[0]): raise ShapeError("Pixels must evenly circumference")
	if not np.allclose(vx,vx[0]): raise ShapeError("Different phi0 per row indicates non-cylindrical enmap")
	phi0 = vx[0]
	# Make a map with the same geometry covering a whole band around the sky
	# We can do this simply by extending it in the positive pixel dimension.
	oshape = map.shape[:-1]+(nphi,)
	owcs   = map.wcs
	# Our input map could in theory cover multiple copies of the sky, which
	# 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.
github simonsobs / pixell / pixell / curvedsky.py View on Github external
def rand_map(shape, wcs, ps, lmax=None, dtype=np.float64, seed=None, oversample=2.0, spin=[0,2], method="auto", direct=False, verbose=False):
	"""Generates a CMB realization with the given power spectrum for an enmap
	with the specified shape and WCS. This is identical to enlib.rand_map, except
	that it takes into account the curvature of the full sky. This makes it much
	slower and more memory-intensive. The map should not cross the poles."""
	# Ensure everything has the right dimensions and restrict to relevant dimensions
	ps = utils.atleast_3d(ps)
	if not ps.shape[0] == ps.shape[1]: raise ShapeError("ps must be [ncomp,ncomp,nl] or [nl]")
	if not (len(shape) == 2 or len(shape) == 3): raise ShapeError("shape must be (ncomp,ny,nx) or (ny,nx)")
	ncomp = 1 if len(shape) == 2 else shape[-3]
	ps = ps[:ncomp,:ncomp]

	ctype = np.result_type(dtype,0j)
	if verbose: print("Generating alms with seed %s up to lmax=%d dtype %s" % (str(seed), lmax, np.dtype(ctype).char))
	alm   = rand_alm_healpy(ps, lmax=lmax, seed=seed, dtype=ctype)
	if verbose: print("Allocating output map shape %s dtype %s" % (str((ncomp,)+shape[-2:]), np.dtype(dtype).char))
	map   = enmap.empty((ncomp,)+shape[-2:], wcs, dtype=dtype)
	alm2map(alm, map, spin=spin, oversample=oversample, method=method, direct=direct, verbose=verbose)
	if len(shape) == 2: map = map[0]
	return map
github simonsobs / pixell / pixell / curvedsky.py View on Github external
def rand_map(shape, wcs, ps, lmax=None, dtype=np.float64, seed=None, oversample=2.0, spin=[0,2], method="auto", direct=False, verbose=False):
	"""Generates a CMB realization with the given power spectrum for an enmap
	with the specified shape and WCS. This is identical to enlib.rand_map, except
	that it takes into account the curvature of the full sky. This makes it much
	slower and more memory-intensive. The map should not cross the poles."""
	# Ensure everything has the right dimensions and restrict to relevant dimensions
	ps = utils.atleast_3d(ps)
	if not ps.shape[0] == ps.shape[1]: raise ShapeError("ps must be [ncomp,ncomp,nl] or [nl]")
	if not (len(shape) == 2 or len(shape) == 3): raise ShapeError("shape must be (ncomp,ny,nx) or (ny,nx)")
	ncomp = 1 if len(shape) == 2 else shape[-3]
	ps = ps[:ncomp,:ncomp]

	ctype = np.result_type(dtype,0j)
	if verbose: print("Generating alms with seed %s up to lmax=%d dtype %s" % (str(seed), lmax, np.dtype(ctype).char))
	alm   = rand_alm_healpy(ps, lmax=lmax, seed=seed, dtype=ctype)
	if verbose: print("Allocating output map shape %s dtype %s" % (str((ncomp,)+shape[-2:]), np.dtype(dtype).char))
	map   = enmap.empty((ncomp,)+shape[-2:], wcs, dtype=dtype)
	alm2map(alm, map, spin=spin, oversample=oversample, method=method, direct=direct, verbose=verbose)
	if len(shape) == 2: map = map[0]
	return map
github simonsobs / pixell / pixell / curvedsky.py View on Github external
# which means dec decreasing and ra increasing.
	flipx = map.wcs.wcs.cdelt[0] < 0
	flipy = map.wcs.wcs.cdelt[1] > 0
	if flipx: map = map[...,:,::-1]
	if flipy: map = map[...,::-1,:]
	# Then check if the map satisfies the lat-ring requirements
	ny, nx = map.shape[-2:]
	vy,vx = enmap.pix2sky(map.shape, map.wcs, [np.arange(ny),np.zeros(ny)])
	hy,hx = enmap.pix2sky(map.shape, map.wcs, [np.zeros(nx),np.arange(nx)])
	dx    = hx[1:]-hx[:-1]
	dx    = dx[np.isfinite(dx)] # Handle overextended coordinates

	if not np.allclose(dx,dx[0]): raise ShapeError("Map must have constant phi spacing")
	nphi = utils.nint(2*np.pi/dx[0])
	if not np.allclose(2*np.pi/nphi,dx[0]): raise ShapeError("Pixels must evenly circumference")
	if not np.allclose(vx,vx[0]): raise ShapeError("Different phi0 per row indicates non-cylindrical enmap")
	phi0 = vx[0]
	# Make a map with the same geometry covering a whole band around the sky
	# We can do this simply by extending it in the positive pixel dimension.
	oshape = map.shape[:-1]+(nphi,)
	owcs   = map.wcs
	# Our input map could in theory cover multiple copies of the sky, which
	# 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
github simonsobs / pixell / pixell / curvedsky.py View on Github external
around the sky. Also returns the slice required to recover the
	input map from the output map."""
	# First check if we need flipping. Sharp wants theta,phi increasing,
	# which means dec decreasing and ra increasing.
	flipx = map.wcs.wcs.cdelt[0] < 0
	flipy = map.wcs.wcs.cdelt[1] > 0
	if flipx: map = map[...,:,::-1]
	if flipy: map = map[...,::-1,:]
	# Then check if the map satisfies the lat-ring requirements
	ny, nx = map.shape[-2:]
	vy,vx = enmap.pix2sky(map.shape, map.wcs, [np.arange(ny),np.zeros(ny)])
	hy,hx = enmap.pix2sky(map.shape, map.wcs, [np.zeros(nx),np.arange(nx)])
	dx    = hx[1:]-hx[:-1]
	dx    = dx[np.isfinite(dx)] # Handle overextended coordinates

	if not np.allclose(dx,dx[0]): raise ShapeError("Map must have constant phi spacing")
	nphi = utils.nint(2*np.pi/dx[0])
	if not np.allclose(2*np.pi/nphi,dx[0]): raise ShapeError("Pixels must evenly circumference")
	if not np.allclose(vx,vx[0]): raise ShapeError("Different phi0 per row indicates non-cylindrical enmap")
	phi0 = vx[0]
	# Make a map with the same geometry covering a whole band around the sky
	# We can do this simply by extending it in the positive pixel dimension.
	oshape = map.shape[:-1]+(nphi,)
	owcs   = map.wcs
	# Our input map could in theory cover multiple copies of the sky, which
	# 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)