fix bug where map.resample still used the old corner version

Author: amaurea

pixell/enmap.py

@@ -98,7 +98,7 @@ def preflat(self):
 	def npix(self): return np.prod(self.shape[-2:])
 	@property
 	def geometry(self): return self.shape, self.wcs
-	def resample(self, oshape, off=(0,0), method="fft", border="wrap", corner=False, order=3): return resample(self, oshape, off=off, method=method, border=border, corner=corner, order=order)
+	def resample(self, oshape, off=(0,0), method="fft", border="wrap", corner=True, order=3): return resample(self, oshape, off=off, method=method, border=border, corner=corner, order=order)
 	def project(self, shape, wcs, mode="spline", order=3, border="constant", cval=0, safe=True): return project(self, shape, wcs, order, mode=mode, border=border, cval=cval, safe=safe)
 	def extract(self, shape, wcs, omap=None, wrap="auto", op=lambda a,b:b, cval=0, iwcs=None, reverse=False): return extract(self, shape, wcs, omap=omap, wrap=wrap, op=op, cval=cval, iwcs=iwcs, reverse=reverse)
 	def extract_pixbox(self, pixbox, omap=None, wrap="auto", op=lambda a,b:b, cval=0, iwcs=None, reverse=False): return extract_pixbox(self, pixbox, omap=omap, wrap=wrap, op=op, cval=cval, iwcs=iwcs, reverse=reverse)
@@ -478,7 +478,7 @@ def pix2sky(shape, wcs, pix, safe=True, corner=False, bcheck=False):
 	coords = np.asarray(wcs.wcs_pix2world(*(tuple(pflat)[::-1]+(0,)))[::-1])*get_unit(wcs)
 	coords = coords.reshape(pix.shape)
 	if safe and not wcsutils.is_plain(wcs):
-		coords = utils.unwind(coords, refmode="middle")
+		coords[1] = utils.unwind(coords[1], refmode="middle")
 	return coords
 
 def sky2pix(shape, wcs, coords, safe=True, corner=False, bcheck=False):

pixell/reproject.py

@@ -8,8 +8,8 @@
 except NameError: basestring = str
 
 def thumbnails(imap, coords, r=5*utils.arcmin, res=None, proj=None, apod=2*utils.arcmin,
-		order=3, oversample=4, pol=None, oshape=None, owcs=None, extensive=False, verbose=False,
-		filter=None,pixwin=False,pixwin_order=0):
+		method="mixed", order=3, oversample=4, pol=None, oshape=None, owcs=None,
+		extensive=False, verbose=False, filter=None,pixwin=False,pixwin_order=0):
 	"""Given an enmap [...,ny,nx] and a set of coordinates in a numpy array
 	coords with shape (n,2) and ordering [n,{dec,ra}], extract a set
 	of thumbnail images [n,...,thumby,thumbx] centered on each set of
@@ -73,7 +73,7 @@ def thumbnails(imap, coords, r=5*utils.arcmin, res=None, proj=None, apod=2*utils
 	# Define our output maps, which we will fill below
 	omaps = enmap.zeros((nsrc,)+imap.shape[:-2]+oshape, owcs, imap.dtype)
 	for si, pixbox in enumerate(pixboxes):
-		if oversample > 1:
+		if method == "nufft" or method == "mixed" and 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])
@@ -86,13 +86,16 @@ def thumbnails(imap, coords, r=5*utils.arcmin, res=None, proj=None, apod=2*utils
 			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:
+		if method == "mixed" and 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
+		# I apologize for the syntax. There should be a better way to do 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, mode="spline", order=order)
+		if   method in ["spline", "mixed"]: ipol_mode = "spline"
+		elif method in ["nufft"]:           ipol_mode = "nufft"
+		else: raise ValueError("Unrecognized method '%s'" % (str(method)))
+		omaps[si] = ithumb.at(ipos, mode=ipol_mode, 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])

pixell/utils.py

@@ -232,6 +232,7 @@ def unwind(a, period=2*np.pi, axes=[-1], ref=0, refmode="left", mask_nan=False):
 	period-wrapping. I.e. [0.07,0.02,6.25,6.20] would
 	become [0.07,0.02,-0.03,-0.08] with the default period
 	of 2*pi."""
+	if a.ndim == 0: return a
 	res = rewind(a, period=period, ref=ref)
 	for axis in axes:
 		with flatview(res, axes=[axis]) as flat:

tests/test_pixell.py

@@ -1078,8 +1078,8 @@ def test_thumbnails(self):
         # with gnomonic/tangent projection
         proj = "tan"
         r = 10*u.arcmin
-        ret = reproject.thumbnails(omap, srcs[:,:2], r=r, res=res, proj=proj, 
-            apod=2*u.arcmin, order=3, oversample=2,pixwin=False)
+        ret = reproject.thumbnails(omap, srcs[:,:2], r=r, res=res, proj=proj,
+            apod=2*u.arcmin, order=3, oversample=4, pixwin=False)
 
         # Create a reference source at the equator to compare this against
         ishape,iwcs = enmap.geometry(shape=ret.shape,res=res,pos=(0,0),proj=proj)