Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def extent(shape, wcs, nsub=None, signed=False, method="auto"):
"""Returns the area of a patch with the given shape
and wcs, in steradians."""
if method == "auto":
if wcsutils.is_plain(wcs): method = "intermediate"
elif wcsutils.is_cyl(wcs): method = "cylindrical"
else: method = "subgrid"
if method in ["inter","intermediate"]: return extent_intermediate(shape, wcs, signed=signed)
elif method in ["cyl", "cylindrical"]: return extent_cyl(shape, wcs, signed=signed)
elif method in ["sub", "subgrid"]: return extent_subgrid(shape, wcs, nsub=nsub, signed=signed)
else: raise ValueError("Unrecognized method '%s' in extent()" % method)
"""Return an enmap where each entry is the coordinate of that entry,
such that posmap(shape,wcs)[{0,1},j,k] is the {y,x}-coordinate of
pixel (j,k) in the map. Results are returned in radians, and
if safe is true (default), then sharp coordinate edges will be
avoided. separable controls whether a fast calculation that assumes that
ra is only a function of x and dec is only a function of y is used.
The default is "auto", which determines this based on the wcs, but
True or False can also be passed to control this manually.
For even greater speed, and to save memory, consider using posaxes directly
for cases where you know that the wcs will be separable. For separable cases,
separable=True is typically 15-20x faster than separable=False, while posaxes
is 1000x faster.
"""
res = zeros((2,)+tuple(shape[-2:]), wcs, dtype)
if separable == "auto": separable = wcsutils.is_cyl(wcs)
if separable:
# If posmap could return a (dec,ra) tuple instead of an ndmap,
# we could have returned np.broadcast_arrays(dec, ra) instead.
# That would have been as fast and memory-saving as broadcast-arrays.
dec, ra = posaxes(shape, wcs, safe=safe, corner=corner)
res[0] = dec[:,None]
res[1] = ra[None,:]
else:
rowstep = int((bsize+shape[-1]-1)//shape[-1])
for i1 in range(0, shape[-2], rowstep):
i2 = min(i1+rowstep, shape[-2])
pix = np.mgrid[i1:i2,:shape[-1]]
res[:,i1:i2,:] = pix2sky(shape, wcs, pix, safe, corner)
return res
def area(shape, wcs, nsamp=1000, method="auto"):
"""Returns the area of a patch with the given shape
and wcs, in steradians."""
if method == "auto":
if wcsutils.is_plain(wcs): method = "intermediate"
elif wcsutils.is_cyl(wcs): method = "cylindrical"
else: method = "contour"
if method in ["inter","intermediate"]: return area_intermediate(shape, wcs)
elif method in ["cyl", "cylindrical"]: return area_cyl(shape, wcs)
elif method in ["cont", "contour"]: return area_contour(shape, wcs, nsamp=nsamp)
else: raise ValueError("Unrecognized method '%s' in area()" % method)
If separable is True, then the map will be assumed to be in a cylindircal
projection, where the pixel size is only a function of declination.
This makes the calculation dramatically faster, and the resulting array
will also use much less memory due to numpy striding tricks. The default,
separable=auto", determines whether to use this shortcut based on the
properties of the wcs.
Normally the function returns a ndmap of shape [ny,nx]. If broadcastable
is True, then it is allowed to return a smaller array that would still
broadcast correctly to the right shape. This can be useful to save work
if one is going to be doing additional manipulation of the pixel size
before using it. For a cylindrical map, the result would have shape [ny,1]
if broadcastable is True.
"""
if separable == True or (separable == "auto" and wcsutils.is_cyl(wcs)):
psize = np.product(pixshapes_cyl(shape, wcs),0)[:,None]
# Expand to full shape unless we are willing to accept an array
# with smaller size that is still broadcastable to the right result
if not broadcastable:
psize = np.broadcast_to(psize, shape[-2:])
return ndmap(psize, wcs)
else:
return np.product(pixshapemap(shape, wcs, bsize=bsize, separable=separable),0)
def distance_from(shape, wcs, points, omap=None, odomains=None, domains=False, method="bubble", rmax=None, step=1024):
"""Find the distance from each pixel in the geometry (shape, wcs) to the
nearest of the points[{dec,ra},npoint], returning a [ny,nx] map of distances.
If domains==True, then it will also return a [ny,nx] map of the index of the point
that was closest to each pixel. If rmax is specified and the method is "bubble", then
distances will only be computed up to rmax. Beyond that distance will be set to rmax
and domains to -1. This can be used to speed up the calculation when one only cares
about nearby areas."""
from pixell import distances
if wcsutils.is_plain(wcs): warnings.warn("Distance functions are not tested on plain coordinate systems.")
if omap is None: omap = empty(shape[-2:], wcs)
if domains and odomains is None: odomains = empty(shape[-2:], wcs, np.int32)
if wcsutils.is_cyl(wcs):
dec, ra = posaxes(shape, wcs)
if method == "bubble":
point_pix = utils.nint(sky2pix(shape, wcs, points))
return distances.distance_from_points_bubble_separable(dec, ra, points, point_pix, rmax=rmax, omap=omap, odomains=odomains, domains=domains)
elif method == "simple":
return distances.distance_from_points_simple_separable(dec, ra, points, omap=omap, odomains=odomains, domains=domains)
else: raise ValueError("Unknown method '%s'" % str(method))
else:
# We have a general geometry, so we need the full posmap. But to avoid wasting memory we
# can loop over chunks of the posmap.
if method == "bubble":
# Not sure how to slice bubble. Just do it in one go for now
pos = posmap(shape, wcs, safe=False)
point_pix = utils.nint(sky2pix(shape, wcs, points))
return distances.distance_from_points_bubble(pos, points, point_pix, rmax=rmax, omap=omap, odomains=odomains, domains=domains)
elif method == "simple":