Merge branch 'master' into gp/feat/meson-support

.github/workflows/build.yml

@@ -32,7 +32,7 @@ jobs:
 
   test-macos:
     name: "Run tests on MacOS"
-    runs-on: macos-latest
+    runs-on: macos-12
     env:
       # LDFLAGS: "-ld64" # For MacOS 13 and above (XCode CLT 15 and above.)
       CC: gcc-12
@@ -125,7 +125,7 @@ jobs:
 
   build_wheels_macos:
     name: Build wheels on MacOS
-    runs-on: macos-latest
+    runs-on: macos-12
     env:
       CC: gcc-12
       CXX: gcc-12
@@ -135,7 +135,7 @@ jobs:
       # Ensure that a wheel builder finishes even if another fails
       fail-fast: false
       matrix:
-        os: [macos-latest]
+        os: [macos-12]
         # We don't build 3.7 wheels for MacOS because that's x86 only.
         cp: [cp38, cp39, cp310, cp311]
         include:

README.rst

@@ -25,7 +25,7 @@ pixell
 Dependencies
 ------------
 
-* Python>=3.7
+* Python>=3.7 but Python 3.12 is not currently supported
 * gcc/gfortran or Intel compilers (clang might not work out of the box), if compiling from source
 * ducc0_, healpy, Cython, astropy, numpy, scipy, matplotlib, pyyaml, h5py, Pillow (Python Image Library)
 

pixell/aberration.py

@@ -100,7 +100,7 @@ def demodulate_map(map, dir=dir_equ, beta=beta,
 
 class Aberrator:
 	def __init__(self, shape, wcs, dir=dir_equ, beta=beta, spin=[0,2],
-			nthread=None, coord_dtype=np.float64, scale_pix=True, boundary="auto"):
+			nthread=None, coord_dtype=np.float64, boundary="auto"):
 		"""Construct an Aberrator object, that can be used to more efficiently
 		aberrate a map. E.g.
 
@@ -117,21 +117,17 @@ def __init__(self, shape, wcs, dir=dir_equ, beta=beta, spin=[0,2],
 		# 1. Calculate the aberration field. These are tiny
 		alm_dpos = calc_boost_field(-beta, dir, nthread=nthread)
 		# 2. Evaluate these on our target geometry. Hardcoded float64 because of get_deflected_angles
-		deflect = enmap.zeros(alm_dpos.shape[:-1]+shape[-2:], wcs, np.float64)
+		deflect = enmap.zeros(alm_dpos.shape[:-1]+shape[-2:], wcs, coord_dtype)
 		curvedsky.alm2map(alm_dpos.astype(coord_ctype, copy=False), deflect, spin=1, nthread=nthread)
 		# 3. Calculate the offset angles.
-		# get_deflected_angles only supports float64 :(
 		rinfo = curvedsky.get_ring_info(shape, wcs)
 		dphi  = np.full(shape[-2], wcs.wcs.cdelt[0]*utils.degree)
-		tmp   = ducc0.misc.get_deflected_angles(theta=rinfo.theta, phi0=rinfo.phi0,
+		odec, ora, gamma = ducc0.misc.get_deflected_angles(theta=rinfo.theta, phi0=rinfo.phi0,
 			nphi=rinfo.nphi, dphi=dphi, ringstart=rinfo.offsets, nthreads=nthread, calc_rotation=True,
 			deflect=np.asarray(deflect).reshape(2,-1).T).T
 		del deflect
-		# We drop down to the target precision as early as possible
-		odec, ora, gamma = tmp.astype(coord_dtype, copy=False)
 		odec = np.pi/2-odec
 		gamma= enmap.ndmap(gamma.reshape(shape[-2:]), wcs)
-		del tmp
 		# 4. Calculate pixel coordinates of offset angles. In general this would use
 		# enmap.sky2pix, but that's slow. Much faster for our typical projections.
 		# Probably worth it to make overrides.
@@ -140,16 +136,12 @@ def __init__(self, shape, wcs, dir=dir_equ, beta=beta, spin=[0,2],
 		fast_rewind(pix[1].reshape(-1), rinfo.nphi[0])
 		del odec, ora
 		# 5. Evaluate the map at these locations using nufft
-		if scale_pix:
-			pix[0] /= shape[-2]
-			pix[1] /= shape[-1]
 		# Boundary condition
 		if boundary == "auto": boundary = "fullsky" if fully(shape, wcs) else "periodic"
 		if boundary not in ["periodic", "fullsky"]:
 			raise ValueError("Unrecognized boundary '%s'" % str(boundary))
 		# 6. Store for when the map is passed in later
 		self.pix       = pix
-		self.scale_pix = scale_pix
 		self.gamma     = gamma
 		self.nthread   = nthread
 		self.spin      = spin
@@ -159,7 +151,7 @@ def __call__(self, map, spin=None, nthread=None):
 		if spin    is None: spin    = self.spin
 		shape, wcs = map.shape, map.wcs
 		ydouble    = self.boundary == "fullsky"
-		map = interpol_map(map, self.pix, nthread=nthread, scaled=self.scale_pix, ydouble=ydouble)
+		map = interpol_map(map, self.pix, nthread=nthread, ydouble=ydouble)
 		map = enmap.ndmap (map.reshape(shape), wcs)
 		# 6. Apply polarization rotation. ducc0.misc.lensing_rotate can do this,
 		# but for some reason it operates on complex numbers instead of a QU field.
@@ -186,10 +178,11 @@ def __init__(self, shape, wcs, dir=dir_equ, beta=beta,
 		compute it from scratch each time."""
 		nthread = int(utils.fallback(utils.getenv("OMP_NUM_THREADS",nthread),0))
 		# 1. Calculate the aberration field. These are tiny
-		alm_dpos, alm_mod = calc_boost_field(-beta, dir, nthread=nthread, modulation=True, mod_exp=-1)
+		alm_mod = calc_boost_field(-beta, dir, nthread=nthread, modulation=True, mod_exp=-1)[1]
+		alm_mod = alm_mod.astype(utils.complex_dtype(dtype), copy=False)
 		# 2: Apply modulation
 		A = enmap.zeros(alm_mod.shape[:-1]+shape[-2:], wcs, dtype)
-		curvedsky.alm2map(alm_mod.astype(utils.complex_dtype(dtype)), A, spin=0, nthread=nthread)
+		curvedsky.alm2map(alm_mod, A, spin=0, nthread=nthread)
 		# Store for __call__
 		self.nthread = nthread;  self.A     = A;     self.modulation = modulation
 		self.T0      = T0;       self.freq  = freq;  self.dipole     = dipole
@@ -250,7 +243,7 @@ def calc_boost_field(beta, dir, lmax=None, nthread=None, modulation=False, mod_e
 	else:
 		return alm
 
-def interpol_map(imap, pixs, epsilon=None, nthread=None, scaled=False, ydouble=False):
+def interpol_map(imap, pixs, epsilon=None, nthread=None, ydouble=False):
 	ny, nx = imap.shape[-2:]
 	if ydouble:
 		# Make y doubled map. This is necessary for the correct boundary condition for
@@ -260,11 +253,7 @@ def interpol_map(imap, pixs, epsilon=None, nthread=None, scaled=False, ydouble=F
 		dmap[...,ny:,:] = np.roll(imap[...,::-1], nx//2, -1)
 	else:
 		dmap = imap
-	if not scaled:
-		pixs[0] /= imap.shape[-2]
-		pixs[1] /= imap.shape[-1]
-	elif ydouble:
-		pixs[0] /= 2
+	periodicity = np.array(dmap.shape[-2:])
 	nthread = int(utils.fallback(utils.getenv("OMP_NUM_THREADS",nthread),0))
 	pflat   = pixs.reshape(2,-1).T
 	if epsilon is None:
@@ -273,18 +262,11 @@ def interpol_map(imap, pixs, epsilon=None, nthread=None, scaled=False, ydouble=F
 	for I in utils.nditer(imap.shape[:-2]):
 		fmap = enmap.fft(dmap[I], normalize=False)
 		oarr[I] = ducc0.nufft.u2nu(grid=np.asarray(fmap), coord=pflat, forward=False,
-			epsilon=epsilon, nthreads=nthread, periodicity=1.0, fft_order=True).real
+			epsilon=epsilon, nthreads=nthread, periodicity=periodicity, fft_order=True).real
 		del fmap
 	# Restore predims
 	oarr = oarr.reshape(imap.shape[:-2]+oarr.shape[-1:])
 	oarr/= dmap.npix
-	# Scale back. Could have worked on a copy instead, but that would
-	# use more memory
-	if not scaled:
-		pixs[0] *= imap.shape[-2]
-		pixs[1] *= imap.shape[-1]
-	elif ydouble:
-		pixs[0] *= 2
 	return oarr
 
 @numba.njit(nogil=True)

pixell/coordinates.py

@@ -286,50 +286,30 @@ def euler_rot(euler_angles, coords, kind="zyz"):
 	co     = utils.rect2ang(rect, False)
 	return co.reshape(coords.shape)
 
-
-def euler_mat(euler_angles, kind="zyz"):
-    """Defines the rotation matrix M for a ABC euler rotation,
-    such that M = A(alpha)B(beta)C(gamma), where euler_angles =
-    [alpha,beta,gamma]. The default kind is ABC=ZYZ."""
-    alpha, beta, gamma = euler_angles
-    R1 = utils.rotmatrix(gamma, kind[2])
-    R2 = utils.rotmatrix(beta,  kind[1])
-    R3 = utils.rotmatrix(alpha, kind[0])
-    return np.einsum("...ij,...jk->...ik",np.einsum("...ij,...jk->...ik",R3,R2),R1)
-
-def euler_rot(euler_angles, coords, kind="zyz"):
-    coords = np.asarray(coords)
-    co     = coords.reshape(2,-1)
-    M      = euler_mat(euler_angles, kind)
-    rect   = utils.ang2rect(co, False)
-    rect   = np.einsum("...ij,j...->i...",M,rect)
-    co     = utils.rect2ang(rect, False)
-    return co.reshape(coords.shape)
-
 def recenter(angs, center, restore=False):
-    """Recenter coordinates "angs" (as ra,dec) on the location given by "center",
-    such that center moves to the north pole."""
-    # Performs the rotation E(0,-theta,-phi). Originally did
-    # E(phi,-theta,-phi), but that is wrong (at least for our
-    # purposes), as it does not preserve the relative orientation
-    # between the boresight and the sun. For example, if the boresight
-    # is at the same elevation as the sun but 10 degrees higher in az,
-    # then it shouldn't matter what az actually is, but with the previous
-    # method it would.
-    #
-    # Now supports specifying where to recenter by specifying center as
-    # lon_from,lat_from,lon_to,lat_to
-    if len(center) == 4: ra0, dec0, ra1, dec1 = center
-    elif len(center) == 2: ra0, dec0, ra1, dec1 = center[0], center[1], center[0]*0, center[1]*0+np.pi/2
-    if restore: ra1 += ra0
-    return euler_rot([ra1,dec0-dec1,-ra0], angs, kind="zyz")
+	"""Recenter coordinates "angs" (as ra,dec) on the location given by "center",
+	such that center moves to the north pole."""
+	# Performs the rotation E(0,-theta,-phi). Originally did
+	# E(phi,-theta,-phi), but that is wrong (at least for our
+	# purposes), as it does not preserve the relative orientation
+	# between the boresight and the sun. For example, if the boresight
+	# is at the same elevation as the sun but 10 degrees higher in az,
+	# then it shouldn't matter what az actually is, but with the previous
+	# method it would.
+	#
+	# Now supports specifying where to recenter by specifying center as
+	# lon_from,lat_from,lon_to,lat_to
+	if len(center) == 4: ra0, dec0, ra1, dec1 = center
+	elif len(center) == 2: ra0, dec0, ra1, dec1 = center[0], center[1], center[0]*0, center[1]*0+np.pi/2
+	if restore: ra1 += ra0
+	return euler_rot([ra1,dec0-dec1,-ra0], angs, kind="zyz")
 
 def decenter(angs, center, restore=False):
-    """Inverse operation of recenter."""
-    if len(center) == 4: ra0, dec0, ra1, dec1 = center
-    elif len(center) == 2: ra0, dec0, ra1, dec1 = center[0], center[1], center[0]*0, center[1]*0+np.pi/2
-    if restore: ra1 += ra0
-    return euler_rot([ra0,dec1-dec0,-ra1],  angs, kind="zyz")
+	"""Inverse operation of recenter."""
+	if len(center) == 4: ra0, dec0, ra1, dec1 = center
+	elif len(center) == 2: ra0, dec0, ra1, dec1 = center[0], center[1], center[0]*0, center[1]*0+np.pi/2
+	if restore: ra1 += ra0
+	return euler_rot([ra0,dec1-dec0,-ra1],  angs, kind="zyz")
 
 def nohor(sys): return sys if sys not in ["altaz","tele","bore"] else "icrs"
 def getsys(sys): return str2sys[sys.lower()] if isinstance(sys,basestring) else sys

pixell/curvedsky.py

@@ -676,8 +676,12 @@ def alm2cl(alm, alm2=None, ainfo=None):
 	ainfo = alm_info(nalm=alm.shape[-1]) if ainfo is None else ainfo
 	return ainfo.alm2cl(alm, alm2=alm2)
 
+euler_angs={}
+euler_angs[("gal","equ")] = np.array([57.06793215,  62.87115487, -167.14056929])*utils.degree
+euler_angs[("equ","gal")] = -euler_angs[("gal","equ")][::-1]
 def rotate_alm(alm, psi, theta, phi, lmax=None, method="auto", nthread=None, inplace=False):
-	"""Rotate the given alm[...,:] via the zyz rotations given by psi, theta and phi.
+	"""Rotate the given alm[...,:] via the zyz rotations given by euler angles
+	psi, theta and phi. See curvedsky.euler_angs for some predefined angles.
 	The underlying implementation is provided by ducc0 or healpy. This is controlled
 	with the "method" argument, which can be "ducc0", "healpy" or "auto". For "auto"
 	it uses ducc0 if available, otherwise healpy. The resulting alm is returned.
@@ -818,6 +822,7 @@ def map2alm_cyl(map, alm=None, ainfo=None, minfo=None, lmax=None, spin=[0,2], we
 			weights/= minfo.ducc_geo.nx
 		else:
 			weights = map.pixsizemap(separable=True, broadcastable=True)[:,0]
+			if minfo.flip: weights = weights[::-1]
 		weights = weights.astype(map.dtype, copy=False)
 	# Loop over pre-pre-dimensions. ducc usually doesn't do anything clever with
 	# these, so looping in python is cheap
@@ -1111,7 +1116,7 @@ def minres_inverse(forward, approx_backward, y, epsilon=1e-6, maxiter=100, zip=N
 	not the fastest choice when only moderate accuracy is needed.
 	"""
 	rhs   = approx_backward(y)
-	rtype = utils.real_dtype(rhs)
+	rtype = utils.real_dtype(rhs.dtype)
 	if zip is None:
 		def zip(a): return a.view(rtype).reshape(-1)
 	if unzip is None:

pixell/enmap.py

@@ -640,6 +640,16 @@ def insert_at(omap, pix, imap, wrap="auto", op=lambda a,b:b, cval=0, iwcs=None):
 	extract_pixbox(omap, pixbox, imap, wrap=wrap, op=op, cval=cval, iwcs=iwcs, reverse=True)
 	return omap
 
+def map_union(map1, map2):
+	"""Given two maps with compatible wcs but possibly covering different
+	parts of the sky, return a new map that contains all pixels of both maps.
+	If the input maps overlap, then those pixels will have the sum of the two maps"""
+	oshape, owcs = union_geometry([map1.geometry, map2.geometry])
+	omap = enmap.zeros(map1.shape[:-2]+oshape[-2:], owcs, map1.dtype)
+	omap.insert(map1)
+	omap.insert(map2, op=lambda a,b:a+b)
+	return omap
+
 def overlap(shape, wcs, shape2_or_pixbox, wcs2=None, wrap="auto"):
 	"""Compute the overlap between the given geometry (shape, wcs) and another *compatible*
 	geometry. This can be either another shape, wcs pair or a pixbox[{from,to},{y,x}].
@@ -1843,6 +1853,7 @@ def distance_from(shape, wcs, points, omap=None, odomains=None, domains=False, m
 	if omap is None: omap = empty(shape[-2:], wcs)
 	if domains and odomains is None: odomains = empty(shape[-2:], wcs, np.int32)
 	points = np.asarray(points)
+	if points.ndim == 1: points = points[:,None]
 	assert points.ndim == 2 and len(points) == 2, "points must be [{dec,ra},npoint]"
 	# Handle case where no points are specified
 	if points.size == 0:
@@ -2397,6 +2408,19 @@ def read_map_geometry(fname, fmt=None, hdu=None, address=None):
 		shape, wcs = slice_geometry(shape, wcs, sel)
 	return shape, wcs
 
+def read_map_dtype(fname, fmt=None, hdu=None, address=None):
+	toks = fname.split(":")
+	fname = toks[0]
+	if fmt == None:
+		if   fname.endswith(".hdf"):     fmt = "hdf"
+		elif fname.endswith(".fits"):    fmt = "fits"
+		elif fname.endswith(".fits.gz"): fmt = "fits.gz"
+		else: fmt = "fits"
+	if   fmt == "fits":    return read_fits_dtype(fname, hdu=hdu)
+	elif fmt == "fits.gz": return read_fits_dtype(fname, hdu=hdu, quick=False)
+	elif fmt == "hdf":     return read_hdf_dtype (fname, address=address)
+	else: raise ValueError
+
 def write_map_geometry(fname, shape, wcs, fmt=None):
 	"""Write an enmap geometry to file. The file type is inferred
 	from the file extension, unless fmt is passed.
@@ -2466,11 +2490,7 @@ def read_fits(fname, hdu=None, sel=None, box=None, pixbox=None, geometry=None, w
 	proxy = ndmap_proxy_fits(hdu, wcs, fname=fname, threshold=sel_threshold, verbose=verbose)
 	return read_helper(proxy, sel=sel, box=box, pixbox=pixbox, geometry=geometry, wrap=wrap, mode=mode, delayed=delayed)
 
-def read_fits_geometry(fname, hdu=None, quick=True):
-	"""Read an enmap wcs from the specified fits file. By default,
-	the map and coordinate system will be read from HDU 0. Use
-	the hdu argument to change this. The map must be stored as
-	a fits image."""
+def read_fits_header(fname, hdu=None, quick=True):
 	if hdu is None: hdu = 0
 	if hdu == 0 and quick:
 		# Read header only, without body
@@ -2483,13 +2503,29 @@ def read_fits_geometry(fname, hdu=None, quick=True):
 	else:
 		with utils.nowarn():
 			header = astropy.io.fits.open(fname)[hdu].header
+	return header
+
+def read_fits_geometry(fname, hdu=None, quick=True):
+	"""Read an enmap wcs from the specified fits file. By default,
+	the map and coordinate system will be read from HDU 0. Use
+	the hdu argument to change this. The map must be stored as
+	a fits image."""
+	header = read_fits_header(fname, hdu=hdu, quick=quick)
 	if header["NAXIS"] < 2:
 		raise ValueError("%s is not an enmap (only %d axes)" % (str(fname), header["NAXIS"]))
 	with warnings.catch_warnings():
 		wcs = wcsutils.WCS(header).sub(2)
 	shape = tuple([header["NAXIS%d"%(i+1)] for i in range(header["NAXIS"])[::-1]])
 	return shape, wcs
 
+def read_fits_dtype(fname, hdu=None, quick=True):
+	header = read_fits_header(fname, hdu=hdu, quick=quick)
+	if "BITPIX" not in header: raise KeyError("BITPIX not defined in fits file")
+	bitpix = header["BITPIX"]
+	table  = {-32:np.float32, -64:np.float64, 8:np.int8, 16:np.int16, 32:np.int32, 64:np.int64}
+	if bitpix not in table: raise ValueError("Unrecognized BITPIX %d" % bitpix)
+	return table[bitpix]
+
 def write_hdf(fname, emap, address=None, extra={}):
 	"""Write an enmap as an hdf file, preserving all the WCS
 	metadata.
@@ -2574,6 +2610,13 @@ def read_hdf_geometry(fname, address=None):
 		shape = hfile["data"].shape
 	return shape, wcs
 
+def read_hdf_dtype(fname, address=None):
+	import h5py
+	with h5py.File(fname,"r") as hfile:
+		if address is not None:
+			hfile = hfile[address]
+		return hfile["data"].dtype
+
 def read_npy(fname, hdu=None, sel=None, box=None, pixbox=None, geometry=None, wrap="auto", mode=None, sel_threshold=10e6, wcs=None, delayed=False, address=None):
 	"""Read an enmap from the specified npy file. Only minimal support.
 	No wcs information."""
@@ -2958,7 +3001,7 @@ def padtiles(*maps, tshape=600, pad=60, margin=60, mode="auto", start=0, step=1)
 		if   len(maps) == 0: mode = ""
 		elif len(maps) == 1: mode = "r"
 		else:                mode = "r"*(len(maps)-1)+"w"
-	tiler = Padtiler(tshape=tshape, pad=pad, margin=margin)
+	tiler = Padtiler(tshape=tshape, pad=pad, margin=margin, start=start, step=step)
 	iters = []
 	for map, io in zip(maps, mode):
 		if   io == "r": iters.append(tiler.read (map))

pixell/enplot.py

@@ -250,6 +250,7 @@ def add_argument(*args, default=None, **kwargs):
 	add_argument("--sub",   type=str, help="Slice a map based on dec1:dec2,ra1:ra2.")
 	add_argument("-H", "--hdu",  type=int, default=0, help="Header unit of the fits file to use")
 	add_argument("--op", type=str, help="Apply this general operation to the map before plotting. For example, 'log(abs(m))' would give you a lograithmic plot.")
+	add_argument("--op2", type=str, help="Like op, but allows multiple statements")
 	add_argument("-d", "--downgrade", type=str, default="1", help="Downsacale the map by this factor before plotting. This is done by averaging nearby pixels. See --upgrade for syntax.")
 	add_argument("--prefix", type=str, default="", help="Specify a prefix for the output file. See --oname.")
 	add_argument("--suffix", type=str, default="", help="Specify a suffix for the output file. See --oname.")
@@ -404,6 +405,11 @@ def get_map(ifile, args, return_info=False, name=None):
 			m1 = m
 			if args.op is not None:
 				m = eval(args.op, {"m":m,"enmap":enmap,"utils":utils,"np":np},np.__dict__)
+			if args.op2 is not None:
+				loc  = {"m":m}
+				glob = {"enmap":enmap,"utils":utils,"np":np}
+				exec(args.op2, glob, loc)
+				m = loc["m"]
 			# Scale if requested
 			scale = parse_list(args.upgrade, int)
 			if np.any(np.array(scale)>1):