Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def make_projectable_map_cyl(map, verbose=False):
"""Given an enmap in a cylindrical projection, return a map with
the same pixelization, but extended to cover a whole band in phi
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
def get_rotated_pixels(shape_source, wcs_source, shape_target, wcs_target,
inverse=False, pos_target=None,
center_target=None, center_source=None):
""" Given a source geometry (shape_source,wcs_source)
return the pixel positions in the target geometry (shape_target,wcs_target)
if the source geometry were rotated such that its center lies on the center
of the target geometry.
WARNING: Only currently tested for a rotation along declination
from one CAR geometry to another CAR geometry.
"""
# what are the center coordinates of each geometries
if center_source is None:
center_source = enmap.pix2sky(
shape_source, wcs_source,
(shape_source[0] / 2., shape_source[1] / 2.))
if center_target is None:
center_target = enmap.pix2sky(
shape_target, wcs_target,
(shape_target[0] / 2., shape_target[1] / 2.))
decs, ras = center_source
dect, rat = center_target
# what are the angle coordinates of each pixel in the target geometry
if pos_target is None:
pos_target = enmap.posmap(shape_target, wcs_target)
#del pos_target
# recenter the angle coordinates of the target from the target center
# to the source center
if inverse:
transfun = lambda x: coordinates.decenter(x, (rat, dect, ras, decs))
def make_projectable_map_cyl(map, verbose=False):
"""Given an enmap in a cylindrical projection, return a map with
the same pixelization, but extended to cover a whole band in phi
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 = [], []
center_target=None, center_source=None):
""" Given a source geometry (shape_source,wcs_source)
return the pixel positions in the target geometry (shape_target,wcs_target)
if the source geometry were rotated such that its center lies on the center
of the target geometry.
WARNING: Only currently tested for a rotation along declination
from one CAR geometry to another CAR geometry.
"""
# what are the center coordinates of each geometries
if center_source is None:
center_source = enmap.pix2sky(
shape_source, wcs_source,
(shape_source[0] / 2., shape_source[1] / 2.))
if center_target is None:
center_target = enmap.pix2sky(
shape_target, wcs_target,
(shape_target[0] / 2., shape_target[1] / 2.))
decs, ras = center_source
dect, rat = center_target
# what are the angle coordinates of each pixel in the target geometry
if pos_target is None:
pos_target = enmap.posmap(shape_target, wcs_target)
#del pos_target
# recenter the angle coordinates of the target from the target center
# to the source center
if inverse:
transfun = lambda x: coordinates.decenter(x, (rat, dect, ras, decs))
else:
transfun = lambda x: coordinates.recenter(x, (rat, dect, ras, decs))
res = coordinates.transform_meta(transfun, pos_target[1::-1], fields=["ang"])
pix_new = enmap.sky2pix(shape_source, wcs_source, res.ocoord[1::-1])
def get_pixsize_rect(shape, wcs):
"""Return the exact pixel size in steradians for the rectangular cylindrical
projection given by shape, wcs. Returns area[ny], where ny = shape[-2] is the
number of rows in the image. All pixels on the same row have the same area."""
ymin = enmap.sky2pix(shape, wcs, [-np.pi / 2, 0])[0]
ymax = enmap.sky2pix(shape, wcs, [np.pi / 2, 0])[0]
y = np.arange(shape[-2])
x = y * 0
dec1 = enmap.pix2sky(shape, wcs, [np.maximum(ymin, y - 0.5), x])[0]
dec2 = enmap.pix2sky(shape, wcs, [np.minimum(ymax, y + 0.5), x])[0]
area = np.abs((np.sin(dec2) - np.sin(dec1)) * wcs.wcs.cdelt[0] * np.pi / 180)
return area
def get_pixsize_rect(shape, wcs):
"""Return the exact pixel size in steradians for the rectangular cylindrical
projection given by shape, wcs. Returns area[ny], where ny = shape[-2] is the
number of rows in the image. All pixels on the same row have the same area."""
ymin = enmap.sky2pix(shape, wcs, [-np.pi / 2, 0])[0]
ymax = enmap.sky2pix(shape, wcs, [np.pi / 2, 0])[0]
y = np.arange(shape[-2])
x = y * 0
dec1 = enmap.pix2sky(shape, wcs, [np.maximum(ymin, y - 0.5), x])[0]
dec2 = enmap.pix2sky(shape, wcs, [np.minimum(ymax, y + 0.5), x])[0]
area = np.abs((np.sin(dec2) - np.sin(dec1)) * wcs.wcs.cdelt[0] * np.pi / 180)
return area