Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
width = 2
if atype in ["c","circle"]:
x,y = topix(annot[1:5])
if skippable(x,y): continue
rad = 8
if len(annot) > 5: rad = int(annot[5])
if len(annot) > 6: width = int(annot[6])
if len(annot) > 7: color = annot[7]
antialias = 1 if width < 1 else 4
draw_ellipse(img,
(x-rad,y-rad,x+rad,y+rad),
outline=color,width=width, antialias=antialias)
elif atype in ["l","line"] or atype in ["r","rect"]:
x1,y1 = topix(annot[1:5])
x2,y2 = topix(annot[5:9])
nphi = utils.nint(abs(360/map.wcs.wcs.cdelt[0]))
x1, x2 = utils.unwind([x1,x2], nphi, ref=nphi//2)
if skippable(x1,y1) and skippable(x2,y2): continue
if len(annot) > 9: width = int(annot[9])
if len(annot) > 10: color = annot[10]
if atype[0] == "l":
draw.line((x1,y1,x2,y2), fill=color, width=width)
else:
if x2 < x1: x1,x2 = x2,x1
if y2 < y1: y1,y2 = y2,y1
for i in range(width):
draw.rectangle((x1+i,y1+i,x2-i,y2-i), outline=color)
elif atype in ["t", "text"]:
x,y = topix(annot[1:5])
if skippable(x,y): continue
text = annot[5]
size = 16
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":
geo = Geometry(shape, wcs)
for y in range(0, shape[-2], step):
sub_geo = geo[y:y+step]
pos = posmap(*sub_geo, safe=False)
if domains:
distances.distance_from_points_simple(pos, points, omap=omap[y:y+step], odomains=odomains[y:y+step], domains=True)
else:
distances.distance_from_points_simple(pos, points, omap=omap[y:y+step])
if domains: return omap, odomains
else: return omap
# 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
return np.moveaxis([centers-rpix,center+rpix+1],0,1)
poss = np.asarray(poss)
res = np.zeros([len(poss),2,2])
for i, pos in enumerate(poss):
# Find the coordinate box we need
dec, ra = pos[:2]
dec1, dec2 = max(dec-r,-np.pi/2), min(dec+r,np.pi/2)
with utils.nowarn():
scale = 1/min(np.cos(dec1), np.cos(dec2))
dra = min(r*scale, np.pi)
ra1, ra2 = ra-dra, ra+dra
box = np.array([[dec1,ra1],[dec2,ra2]])
# And get the corresponding pixbox
res[i] = skybox2pixbox(shape, wcs, box)
# Turn ranges into from-inclusive, to-exclusive integers.
res = utils.nint(res)
res = np.sort(res, 1)
res[:,1] += 1
return res
def pixbox_of(iwcs,oshape,owcs):
"""Obtain the pixbox which when extracted from a map with WCS=iwcs
returns a map that has geometry oshape,owcs.
"""
# First check that our wcs is compatible
assert wcsutils.is_compatible(iwcs, owcs), "Incompatible wcs in enmap.extract: %s vs. %s" % (str(iwcs), str(owcs))
# Find the bounding box of the output in terms of input pixels.
# This is simple because our wcses are compatible, so they
# can only differ by a simple pixel offset. Here pixoff is
# pos_input - pos_output. This is equivalent to the coordinates of
pixoff = utils.nint((iwcs.wcs.crpix-owcs.wcs.crpix) - (iwcs.wcs.crval-owcs.wcs.crval)/iwcs.wcs.cdelt)[::-1]
pixbox = np.array([pixoff,pixoff+np.array(oshape[-2:])])
return pixbox
def extract_pixbox(map, pixbox, omap=None, wrap="auto", op=lambda a,b:b, cval=0, iwcs=None, reverse=False):
"""This function extracts a rectangular area from an enmap based on the
given pixbox[{from,to,[stride]},{y,x}]. The difference between this function
and plain slicing of the enmap is that this one supports wrapping around the
sky. This is necessary to make things like fast thumbnail or tile extraction
at the edge of a (horizontally) fullsky map work."""
if iwcs is None: iwcs = map.wcs
pixbox = np.asarray(pixbox)
if omap is None:
oshape, owcs = slice_geometry(map.shape, iwcs, (slice(*pixbox[:,-2]),slice(*pixbox[:,-1])), nowrap=True)
omap = full(map.shape[:-2]+tuple(oshape[-2:]), owcs, cval, map.dtype)
nphi = utils.nint(360/np.abs(iwcs.wcs.cdelt[0]))
# If our map is wider than the wrapping length, assume we're a lower-spin field
nphi *= (nphi+map.shape[-1]-1)//nphi
if utils.streq(wrap, "auto"):
wrap = [0,0] if wcsutils.is_plain(iwcs) else [0,nphi]
else: wrap = np.zeros(2,int)+wrap
for ibox, obox in utils.sbox_wrap(pixbox.T, wrap=wrap, cap=map.shape[-2:]):
islice = utils.sbox2slice(ibox)
oslice = utils.sbox2slice(obox)
if reverse: map [islice] = op(map[islice], omap[oslice])
else: omap[oslice] = op(omap[oslice], map[islice])
return omap
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
omaps = enmap.zeros((nsrc,)+imap.shape[:-2]+oshape, owcs, imap.dtype)
for si, pixbox in enumerate(pixboxes):
if oversample > 1:
# Make the pixbox fft-friendly
for i in range(2):
pixbox[1,i] = pixbox[0,i] + fft.fft_len(pixbox[1,i]-pixbox[0,i], direction="above", factors=[2,3,5])
ithumb = imap.extract_pixbox(pixbox)
if extensive: ithumb /= ithumb.pixsizemap()
ithumb = ithumb.apod(apod_pix, fill="median")
if filter is not None: ithumb = filter(ithumb)
if verbose:
print("%4d/%d %6.2f %6.2f %8.2f %dx%d" % (si+1, nsrc, coords[si,0]/utils.degree, coords[si,1]/utils.degree, np.max(ithumb), ithumb.shape[-2], ithumb.shape[-1]))
# Oversample using fourier if requested. We do this because fourier
# interpolation is better than spline interpolation overall
if oversample > 1:
fshape = utils.nint(np.array(oshape[-2:])*oversample)
ithumb = ithumb.resample(fshape, method="fft")
# I apologize for the syntax. There should be a better way of doing this
ipos = coordinates.transform("cel", ["cel",[[0,0,coords[si,1],coords[si,0]],False]], opos[::-1], pol=pol)
ipos, rest = ipos[1::-1], ipos[2:]
omaps[si] = ithumb.at(ipos, order=order)
# Apply the polarization rotation. The sign is flipped because we computed the
# rotation from the output to the input
if pol: omaps[si] = enmap.rotate_pol(omaps[si], -rest[0])
if extensive: omaps *= omaps.pixsizemap()
# Restore original dimension
omaps = omaps.reshape(ishape + omaps.shape[1:])
return omaps