Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# 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.
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
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
# 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
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)