How to use the pixell.utils.nint 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 / enplot.py View on Github external
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
github simonsobs / pixell / pixell / enmap.py View on Github external
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
github simonsobs / pixell / pixell / curvedsky.py View on Github external
# 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
github simonsobs / pixell / pixell / enmap.py View on Github external
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
github simonsobs / pixell / pixell / enmap.py View on Github external
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
github simonsobs / pixell / pixell / enmap.py View on Github external
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
github simonsobs / pixell / pixell / curvedsky.py View on Github external
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
github simonsobs / pixell / pixell / reproject.py View on Github external
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