Pointing and Projections (so3g.proj)

Overview

This module contains fast and/or flexible routines for pointing computations that are required in the context of time-ordered data and map-making. The main focus is on rapid conversion from horizon coordinates to celestial coordinates. On-the-fly conversion to projection plane coordinates and accumulation of detector data into maps (without storing long position or pixel vectors) is also supported.

The basic acceleration approach used here is to express the rotation that takes a detector to the celestial sky as a product of a (fixed) rotation that take each detector to some reference position in the focal plane and a rotation that takes the reference position in the focal plane to the celestial sky, as a function of time. This is an approximation insofar as non-linear effects such as atmospheric refraction and aberration add an additional dependence on the horizon and celestial pointing of a detector. But having accounted for the mean effects, at the reference position, the errors incurred on nearby detectors in the focal plane may be small enough to ignore.

Dependencies

To take full advantage of this module, you will want to install these packages:

• https://github.com/arahlin/qpoint : for high-precision astrometry.

• https://github.com/simonsobs/pixell : for rectilinear sky map operations.

Tutorial

Import

Because the Projection module has special requirements, it must be explicitly imported for use. It is likely you’ll also need numpy:

import so3g.proj
import numpy as np

Horizon to equatorial coordinates

The CelestialSightLine class is used to store a time-dependent rotation from focal plane coordinates to celestial coordinates. To create such an object based on the horizon coordinates of a telescope boresight, you must pass in the az, el, and time (as a unix timestamp), as well as a description of the telescope site and the weather conditions:

# Define DEG, carrying conversion from degrees to radians.
DEG = np.pi / 180

# Vectors of time, azimuth, and elevation.
t = [1900000000.]
az = [180. * DEG]
el = [60. * DEG]

# Construct a SightLine from this information.
csl = so3g.proj.CelestialSightLine.az_el(t, az, el,
    site='so', weather='typical')

After instantiation, csl will have a member Q that holds the rotation quaternion from focal plane to celestial coordinates. The coords() method can decompose such quaternions into rotation angles:

>>> csl.Q
spt3g.core.G3VectorQuat([(-0.0384748,0.941783,0.114177,0.313891)])

>>> csl.coords()
array([[ 0.24261138, -0.92726871, -0.99999913, -0.00131952]])

The coords() returns an array with shape (n_time, 4); each 4-tuple contains values (lon, lat, cos(gamma), sin(gamma)). The lon and lat are the celestial longitude and latitude coordinates (typically Right Ascension and Declination), in radians, and gamma is the parallactic angle (the angle between the directions of increasing declination and increasing elevation at that point, with 0 corresponding to North, increasing towards West).

You can get the vectors of RA and dec, in degrees like this:

>>> ra, dec = csl.coords().transpose()[:2] / DEG
>>> print(ra, dec)
[13.90060809] [-53.12858329]

Pointing for many detectors

Create a FocalPlane object, with some detector positions and orientations:

xi  = np.array([-0.5, 0.0, 0.5]) * DEG
eta = np.zeros(3)
gamma = np.array([0,30,60]) * DEG
fp  = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma)

This particular function, from_xieta(), will apply the SO standard coordinate definitions and stores quaternions (.quats) and responsivities (.resps) for each detector. These are a G3VectorQuat wth length ndet and a numpy array with shape [ndet,2] respectively. fp.quats[0] gives the quaternion for the first detector, while fp.resps[0] gives its total intensity and polarization responsivity.:

>>> fp.quats[2]
spt3g.core.quat(0.866017,0.00377878,-0.00218168,0.499995)
>>> fp.resps[2]
array([1., 1.], dtype=float32)

As you can see, the default responsivity is 1 for both total intensty and polarization. To represent detectors with responsivity different from 1, use the T and P arguments to from_xieta() to set the total intensity and polarization responsivity respectively. These can be either single numbers or array-likes with lengths ndet.:

xi  = np.array([-0.5, 0.0, 0.5]) * DEG
eta = np.zeros(3)
gamma = np.array([0,30,60]) * DEG
T   = np.array([1.0, 0.9, 1.1])
P   = np.array([0.5, 0.6, 0.05])
fp2 = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma, T=T, P=P)

Together, gamma, T and P specify the full responsivity of a detector to the T, Q and U Stokes parameters in focal plane coordinates. But as an alternative, it’s also possible to specify these directly. The example above is equivalent to:

xi  = np.array([-0.5, 0.0, 0.5]) * DEG
eta = np.zeros(3)
gamma = np.array([0,30,60]) * DEG
T   = np.array([1.0, 0.9, 1.1])
Q   = np.array([0.5, 0.3, -0.025])
U   = np.array([0.0, 0.51961524, 0.04330127])
fp2 = so3g.proj.FocalPlane.from_xieta(xi, eta, T=T, Q=Q, U=U)

At this point you could get the celestial coordinates for any one of those detectors:

# Get vector of quaternion pointings for detector 0
q_total = csl.Q * fp.quats[0]
# Decompose q_total into lon, lat, roll angles
ra, dec, gamma = so3g.proj.quat.decompose_lonlat(q_total)

As expected, these coordinates are close to the ones computed before, for the boresight:

>>> print(ra / DEG, dec / DEG)
[14.73387143] [-53.1250149]

But the more expedient way to get pointing for multiple detectors is to call coords() with the FocalPlane object as first argument:

>>> csl.coords(fp)
[array([[ 0.25715457, -0.92720643,  0.9999161 ,  0.01295328]]),
 array([[ 0.24261138, -0.92726871,  0.86536489,  0.5011423 ]]),
 array([[ 0.22806774, -0.92722945,  0.50890634,  0.86082189]])]

So coords() returns a list of numpy arrays with shape (n_time,4), and at each time step the 4 elements of the array are: (lon, lat, cos(gamma), sin(gamma)).

Projecting to Maps

Accelerated projection routines for various pixelizations have been written using C++ templates. Abstraction at the Python level is provided by the Projectionist class, which decodes the map WCS description in order to choose the right accelerated pointing projection routines (encapsulated in a C++ class called a ProjEng).

This code relies on pixell library for working with sky maps that are represented as numpy arrays tied to astropy WCS objects describing the rectangular pixelization.

Let’s start by creating a suitable map; we choose a supported projection, and make sure it contains the particular points touched by our example from above:

from pixell import enmap
shape, wcs = enmap.geometry([[-54*DEG, 16.*DEG], [-52*DEG, 12.*DEG]],
  res=.02*DEG, proj='car')

Note that arguments to most enmap functions are specified in (dec, RA) order, because of the typical correspondence between the RA direction and matrix column (the fastest index of the map array). The shape is a simple tuple, and will be used like a numpy array shape. The wcs is an astropy WCS object, but decorated with a short descriptive __repr__ provided by pixell:

>>> shape
(100, 200)
>>> wcs
car:{cdelt:[-0.02,0.02],crval:[14,0],crpix:[101,2701]}

In pixell nomenclature, the combination (shape, wcs) is called a geometry and is used in many different enmap functions.

Let us create a Projectionist that can be used to project between the time-domain of our detectors and the map specified by this geometry:

p = so3g.proj.Projectionist.for_geom(shape, wcs)
asm = so3g.proj.Assembly.attach(csl, fp)

As a first example, let’s project a map into the time domain. For a map to project from, let’s ask enmap to give us the maps that specifiy the coordinate values of every pixell (this will be an array with shape (2,100,200), where the first index (0 or 1) selects the coordinate (dec or RA)):

pmap = enmap.posmap(shape, wcs)

Now projecting into the time domain:

pix_dec = p.from_map(pmap[0], asm)
pix_ra = p.from_map(pmap[1], asm)

Inspecting the values, we see they are close to the coordinates we expect. That makes sense, as these values correspond to the values from the pixels in pmap, which are the coordinates of the centers of the pixels:

>>> [x/DEG for x in pix_ra]
[array([14.740001], dtype=float32), array([13.900001], dtype=float32), array([13.059999], dtype=float32)]
>>> [x/DEG for x in pix_dec]
[array([-53.12], dtype=float32), array([-53.12], dtype=float32), array([-53.12], dtype=float32)]

If you are not getting what you expect, you can grab the pixel indices inferred by the projector – perhaps your pointing is taking you off the map (in which case the pixel indices would return value -1):

>>> p.get_pixels(asm)
[array([[44, 63]], dtype=int32), array([[ 44, 105]], dtype=int32),
array([[ 44, 147]], dtype=int32)]

Let’s project signal into an intensity map using Projectionist.to_map():

# Create dummy signal for our 3 detectors at 1 time points:
signal = np.array([[1.], [10.], [100.]], dtype='float32')

# Project into T-only map.
map_out = p.to_map(signal, asm, comps='T')

Inspecting the map, we see our signal values occupy the three non-zero pixels:

>>> map_out.nonzero()
(array([0, 0, 0]), array([44, 44, 44]), array([ 63, 105, 147]))
>>> map_out[map_out!=0]
array([  1.,  10., 100.])

If we run this projection again, but pass in this map as a starting point, the signal will be added to the map:

>>> p.to_map(signal, asm, output=map_out, comps='T')
array([[[0., 0., 0., ..., 0.]]])
>>> map_out[map_out!=0]
array([  2.,  20., 200.])

If we instead want to treat the signal as coming from polarization-sensitive detectors, we can request components 'TQU':

map_pol = p.to_map(signal, asm, comps='TQU')

Now a 3-dimensional output map is created, with values binned according to the projected detector angle on the sky:

>>> map_pol.shape
(3, 100, 200)
>>> map_pol[:,44,105]
array([10.,  4.97712803, 8.673419])

For the most basic map-making, the other useful operation is the Projectionist.to_weights() method. This is used to compute the weights matrix in each pixel, returned as a (n_comp,n_comp,n_dec,n_ra) map. Mathematically this corresponds to:

\[W = s^T\,P\,P^T\,s\]

where s is the time-ordered signal matrix with value 1 at every position. (If det_weights is specified in the call, then s is filled with the detector weight.)

The weights map can be obtained like this:

weight_out = p.to_weights(asm, comps='TQU')

The shape corresponds to a 3x3 matrix at each pixel; but notice only the upper diagonal has been filled in, for efficiency reasons…:

>>> weight_out.shape
(3, 3, 100, 200)
>>> weight_out[...,44,105]
array([[1.        , 0.49771279, 0.86734194],
       [0.        , 0.24771802, 0.43168718],
       [0.        , 0.        , 0.75228202]])

OpenMP

The three routines of the Projectionist that move data between map and time-domain are: from_map, to_map, and to_weights. The from_map function will automatically use OMP threading. The to_map and to_weights functions need to be instructed on how to use OpenMP safely, or else they will default to a single thread implementation.

The reason that thread-safety is non-trivial for to_map and to_weights is that one must be careful that different threads to not try to update the same pixel value at the ~same time. Some common approaches to dealing with this race condition are:

• Maintain separate copies of the output map for each thread, then add them together at the end. This is memory-intensive, requiring a separate copy of the map for each thread, and requires post-processing to combine the maps.

• Use atomic update operations. This requires locking, which can be very slow to somewhat slow depending on the architecture.

The approach that is implemented in so3g.proj is to assign different regions of the map (i.e. different groups of pixels) to specific threads. In a precomputation step, the code determines and caches, for each detector, which samples correspond to which pixel group. Then, in TOD-to-map operations, each OpenMP thread loops over all detectors but only processes the samples belonging to that thread. The result is that each thread touches a disjoint set of pixels.

The assignment of ranges of samples to threads is called a thread assignment and is passed to Projectionist methods through the threads= argument. The thread assignment is carried in a RangesMatrix object, with shape (n_threads, n_dets, n_samps).

The thread assignment approach requires precomputation, which can itself be somewhat expensive, so this solution is targeting use cases where the projection operations need to be applied multiple times. The best way to assign pixels to threads depends on the particulars of the scan pattern. Below we describe a few ways to generate the thread assignment.

Although it is not perfectly general, for constant elevation scans the domdir method is has proven quite successful so you might want to start there.

The function Projectionist.assign_threads() can be used to compute thread assignments using a few different algorithms. For example:

threads = p.assign_threads(asm, method='domdir')
map_pol2 = p.to_map(signal, asm, comps='TQU', threads=threads)

Inspecting:

>>> threads
RangesMatrix(4,3,1)
>>> map_pol2[:,44,105]
array([10.        ,  4.97712898,  8.67341805])

The same threads result can be passed to p.to_weights.

In the sections below, a bit more detail is provided on the thread assignment options.

Simple thread assignment (simple)

The simple algorithm works by dividing the map into n_threads equal-sized horizontal bands, and assigning each band to a thread. This works well if the scan is horizontal and the data fill the entire map area. It tends to perform poorly otherwise, usually because the weight in the map is heavier in the center than the edges. Access this by passing method='simple' to Projectionist.assign_threads().

Dominant Direction (domdir)

The domdir algorithm is a generalization of the simple algorithm that accounts for the dominant scan direction and for imbalance of weight in the map. Determining the scan direction involves a TOD-to-map projection operation, but can be accomplished with a decimated focal plane. The balancing of weights in threads requires a projection from map-to-TOD, which requires creating a new TOD-sized data structure (temporarily).

You can access this method by passing method='domdir' to Projectionist.assign_threads(). For finer-grained control, see so3g.proj.mapthreads.get_threads_domdir().

Tile-based thread assignment (tiles)

In Tiled mapping, the tiles provide convenient grouping of pixels, and can be used as the units for thread assignments. The tiles can be assigned to a thread, accounting for the number of samples ending up in each tile to balance the work. However, a good balance is probably not possible unless you have managed to end up with many more active tiles than threads.

You can access this method by passing method='tiles' to Projectionist.assign_threads(). The thread assignment in this case is computed by Projectionist.get_active_tiles().

Thread assignment from map

If you want to “manually” assign pixels to threads, consider using Projectionist.assign_threads_from_map(). You must first create a single component map where the value of the pixel is the index of the thread to which you want to assign the pixel. Using this map, the code determines the thread index of every sample in the TOD, and returns the RangesMatrix with the result. (This function is used by the domdir method.)

For example, we might know that slicing the map into wide columns is a good way to assign threads. So we make a map with 0s, 1s, 2s and 3s in broad columns, and than use that to produce a thread assignment RangesMatrix:

# Get a blank map to start from...
thread_map = map_out * 0

# Stripe it.
full_width = thread_map.shape[-1]
n_threads = 4
for i in range(n_threads):
  thread_map[...,i*full_width//n_threads:] = i

# Convert to RangesMatrix
threads = p.assign_threads_from_map(asm, thread_map)

Tiled Maps

This projection code supports tiling of maps. Tiles cover some parent geometry. Suppose the parent geometry is obtained like this:

from pixell import enmap
shape, wcs = enmap.geometry([[-54*DEG, 16.*DEG], [-52.2*DEG, 12.2*DEG]],
  res=.02*DEG, proj='car')

As above, this produces a footprint with a particular shape:

>>> shape
(90, 190)

And we use classmethod Projectionist.for_tiled() to specify the parent geometry as well as the tile size (which is (20, 50) here):

p = so3g.proj.Projectionist.for_tiled(shape, wcs, (20, 50))

The 3rd argument gives the shape of the tile. In cases where the tiles don’t precisely cover the parent footprint, the tiles at the end of each row or column will be smaller than the specified tile size. Tiles are indexed from 0 to n_tile - 1, and in standard display (where the parent shape corner index (0,0) is in the lower left), the tile index increases to the right and then upwards. See the diagram below, which shows how parent shape (90, 190) is tiled with tile shape (20, 50).

For the projection routines between TOD and map space, arguments or returned values representing a map now return a list of maps, with one entry per size of the tile grid. In the example above, we can get a blank map:

>>> tms = p.zeros((1,))
>>> len(tms)
20
>>> tms[0].shape
(1, 20, 50)
>>> tms[-1].shape
(1, 10, 40)

The advantage of tiling is that one does not need to track all the tiles; to specify which tiles you care about, pass them in as a list to for_tiled. For example, if I know I can leave off the upper right and lower left corners:

p = so3g.proj.Projectionist.for_tiled(shape, wcs, (20, 50),
    [0, 1, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14, 15, 18, 19])

Such a list can be generated automatically using get_active_tiles(), which can then be used to construct a new Projectionist:

# Construct without tile list...
p = so3g.proj.Projectionist.for_tiled(shape, wcs, (20, 50))

# Get active tiles (does not require creating a map)
tile_list = p.get_active_tiles(asm)

# Re-construct with tiles selected.
p = so3g.proj.Projectionist.for_tiled(shape, wcs, (20, 50), tile_list)

HEALPix Maps

In addition to rectangular pixelizations, the HEALPix pixelization is supported through the ProjectionistHealpix class. This shares most of the routines of the Projectionist class for rectangular pixelization.

We make a simple ProjectionistHealpix:

nside = 128
p = so3g.proj.ProjectionistHealpix.for_healpix(nside, ordering='NEST')
asm = so3g.proj.Assembly.attach(csl, fp)

The first argument here is HEALPix NSIDE defining the number of pixels. The second is an optional argument defining the pixel ordering; it may be ‘NEST’ or ‘RING’, default ‘NEST’. See HEALPix docs for more info.

As for rectangular pixelization we can use from_map to go from map to time domain:

npix = 12*nside**2
imap = np.arange(npix, dtype=np.float64)  # Make a map
tod = p.from_map(imap, asm)
# We can check the answer with healpy
import healpy as hp
ftod = np.array(tod, dtype=np.int64).flatten()
pix_dec, pix_ra = hp.pix2ang(nside, ftod, nest=True, lonlat=True)

This basic projectionist supports only a very simple threading scheme, 'simple'; this is likely to be highly inefficient for small sky-area maps. It is therefore primarily useful for small cases and from_map. For efficient OpenMP parallelization for to_map and to_weights use a Tiled map:

nside = 128
nside_tile = 4 # power of 2 or 'auto'
p = so3g.proj.ProjectionistHealpix.for_healpix(nside, nside_tile)
threads = p.assign_threads(asm, 'tiles')

Tiles are defined by HEALPix super-pixels at nside_tile, which may be given explicitly or computed automatically for efficient threading. Tiled maps currently only support ‘NEST’ ordering.

To project a map:

# Same dummy signal as above; 3 detectors at 1 time point
signal = np.array([[1.], [10.], [100.]], dtype='float32')
imap = p.to_map(signal, asm, comps='TQU', threads=threads)
weights = p.to_weights(asm, comps='TQU', threads=threads)

A tiled map will be a list of tiles, with None for empty tiles. Converting to a full-sky HEALPix map requires a bit more work, for example:

tile_shape = (3, (nside//nside_tile)**2) # (ncomp, npix_per_tile)
full_map = np.hstack([tile if tile is not None else \
           np.zeros(tile_shape) for tile in imap])

Coordinate Systems

We endeavor to limit the expression of rotations in terms of tuples of angles to three representations: iso, lonlat, and xieta. These are defined and applied in tod2maps_docs/coord_sys. Routines for converting between angle tuples and quaternions are provided in the so3g.proj.quat submodule.

Class and module reference

Several core classes and submodules from so3g.proj should be auto-documented below.

Assembly

class so3g.proj.Assembly(collapse=False)

This class groups together boresight and detector offset information for use by the Projectionist.

The basic implementation simply attachs some number of detectors rigidly to a boresight. But this abstraction layer could facilitate more complex arrangements, eventually.

classmethod attach(sight_line, fplane)

Create an Assembly based on a CelestialSightLine and a FocalPlane.

Parameters:

• sight_line (CelestialSightLine) – The position and orientation of the boresight. This can just be a G3VectorQuat if you don’t have a whole CelestialSightLine handy.

• fplane (FocalPlane) – The “offsets” of each detector relative to the boresight, and their response to intensity and polarization

classmethod for_boresight(sight_line)

Returns an Assembly where a single detector serves as a dummy for the boresight.

CelestialSightLine

class so3g.proj.CelestialSightLine

Carries a vector of celestial pointing data.

Instantiation should result in the pointing quaternion being stored in self.Q.

static decode_site(site=None)

Convert site argument to a Site object. The argument must be one of:

• an EarthlySite object

• a string, corresponding to the name of an internally known site

• None, in which casethe default site info will be loaded.

classmethod naive_az_el(t, az, el, roll=0.0, site=None, weather=None)

Construct a SightLine from horizon coordinates, az and el (in radians) and time t (unix timestamp).

This will be off by several arcminutes… but less than a degree. The weather is ignored.

classmethod az_el(t, az, el, roll=None, site=None, weather=None, **kwargs)

Construct a SightLine from horizon coordinates. This uses high-precision pointing. kwargs are passed to qpoint.Qpoint.

classmethod for_lonlat(lon, lat, psi=0)

Construct a SightLine directly from lonlat angles representing celestial coordinates.

This is appropriate if you want a SightLine that takes focal plane coordinates directly to some celestial pointing. In that case, lon and lat are RA and dec, and psi is the parallactic rotation of the focal plane on the sky at the target position. psi=0 will correspond to focal plane “up” (+eta axis) mapping parallel to lines of longitude; psi > 0 is a clockwise rotation of the focal plane on the sky (i.e. opposite IAU Position Angle direction).

classmethod for_horizon(t, az, el, roll=None, site=None, weather=None)

Construct the trivial SightLine, where “celestial” coordinates are taken to simply be local horizon coordinates (without any effects of atmosphere). Although this accepts site and weather arguments, they are ignored.

coords(fplane=None, output=None)

Get the celestial coordinates of each detector at each time.

Parameters:

fplane – A FocalPlane object representing the detector offsets and responses, or None

output: An optional structure for holding the results. For

that to work, each element of output must support the buffer protocol.

Returns:

If fplane is None, then the result will be [n_samp,{lon,lat,cos2psi,sin2psi}]. Otherwise it will be [n_det,n_samp,{lon,lat,cos2psi,sin2psi}]

EarthlySite

class so3g.proj.EarthlySite(lon, lat, elev, typical_weather=None)

Arguments in degrees E, degrees N, meters above sea level. The typical_weather argument should be an so3g Weather object, or None.

ephem_observer()

Return an ephem.Observer corresponding to this location.

skyfield_site(spice_kernel)

Return a skyfield VectorSum for this location on ‘earth’ (starting from wgs84). The spice_kernel is the thing one might get from jpllib.SpiceKernel(emphemeris_filename).

FocalPlane

class so3g.proj.FocalPlane(quats=None, resps=None, dets=None)

This class represents the detector positions and intensity and polarization responses in the focal plane.

quats

G3VectorQuat representing the pointing quaternions for each detector. Can be turned into a numpy array of coefficients with np.array(). Or call .coeffs() to get them directly.

resps

Array of float32 with shape [ndet,2] representing the total intensity and polarization responsivity of each detector

ndet

The number of detectors (read-only)

(This used to be a subclass of OrderedDict, but this was hard to generalize to per-detector polarization efficiency.)

Construct a FocalPlane from detector quaternions and responsivities.

Parameters:

• quats –

  Detector quaternions. Either:

  • An array-like of floats with shape [ndet,4]

  • An array-like of so3g.proj.quat.quat with shape [ndet]

  • An so3g.proj.quat.G3VectorQuat

  • None, which results in an empty focalplane with no detectors

• resps –

  Detector responsivities. Either:

  • An array-like of floats with shape [ndet,2], where the first and second number in the last axis are the total intensity and polarization response respectively

  • None, which results in a T and P response of 1 for all detectors.

• dets – Deprecated argument temporarily present for backwards compatibility.

coeffs()

Return an [ndet,4] numpy array representing the quaternion coefficients of .quats. Useful for passing the detector pointing to C++ code

classmethod boresight()

Construct a FocalPlane representing a single detector with a responsivity of one pointing along the telescope boresight

classmethod from_xieta(*args, **kwargs)

from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False) Construct a FocalPlane from focal plane coordinates (xi,eta).

For backwards compatibility, the signature from_xieta(names, xi, eta, gamma=0) is also supported; but this will be removed in the future.

These are Cartesian projection plane coordinates. When looking at the sky along the telescope boresight, xi is parallel to increasing azimuth and eta is parallel to increasing elevation. The angle gamma, which specifies the angle of polarization sensitivity, is measured from the eta axis, increasing towards the xi axis.

Parameters:

• xi – An array-like of floats with shape [ndet]

• eta – An array-like of floats with shape [ndet]

• gamma – The detector polarization angles. [ndet], or a scalar

• T – The total intensity and polarization responsivity. Array-like with shape [ndet], or scalar.

• P – The total intensity and polarization responsivity. Array-like with shape [ndet], or scalar.

• Q – The Q- and U-polarization responsivity. Array-like with shape [ndet], or scalar. Q,U are alternatives to P,gamma for specifying the polarization responsivity and angle.

• U – The Q- and U-polarization responsivity. Array-like with shape [ndet], or scalar. Q,U are alternatives to P,gamma for specifying the polarization responsivity and angle.

So there are two ways to specify the polarization angle and responsivity:

1. gamma and P

2. Q and U

Examples, assuming ndet = 2:

from_xieta(xi, eta, gamma=[0,pi/4])

Constructs a FocalPlane with T and P responsivity of 1 and polarization angles of 0 and 45 degrees, representing a Q-sensitive and U-sensitive detector.

from_xieta(xi, eta, gamma=[0,pi/4], P=0.5)

Like the above, but with a polarization responsivity of just 0.5.

from_xieta(xi, eta, gamma=[0,pi/4], T=[1,0.9], P=[0.5,0.6])

Like above, but with a detector-dependent intensity and polarization responsivity. There is no restriction that T > P. For the pseudo-detector timestreams one gets after HWP demodulation, one would have T=0 for the cos-modulated and sin-modulated timestreams, for example.

from_xieta(xi, eta, Q=[1,0], U=[0,1])

Construct the FocalPlane with explicit Q and U responsivity. This example is equivalent to example 1.

Usually one would either use gamma,P or Q,U. If they are combined, then gamma_total = gamma + arctan2(U,Q)/2 and P_tot = P * (Q**2+U**2)**0.5.

items()

Iterate over detector quaternions and responsivities. Yields (spt3g.core.quat, array([Tresp,Presp])) pairs. Unlke the raw entries in the .quats member, which are just numpy arrays with length 4, spt3g.core.quat are proper quaternon objects that support quaternion math.

mapthread

This submodule is for functions that produce “thread assignments”, i.e. disjoint RangesMatrix objects with shape (n_threads, n_dets, n_samps) that are suitable for the threads= argument in Projection code that projects from time to map space using OpenMP parallelization.

so3g.proj.mapthreads.get_num_threads(n_threads=None)

Utility function for computing n_threads. If n_threads is not None, it is returned directly. But if it is None, then the OpenMP thread count is returned. Uses so3g.useful_info().

so3g.proj.mapthreads.get_threads_domdir(sight, fplane, shape, wcs, tile_shape=None, active_tiles=None, n_threads=None, fplane_rep=None, plot_prefix=None)

Assign samples to threads according to the dominant scan direction.

The dominant scan direction is first determined, which requires creating a 3-component map. This projection operation can’t be parallelized safely so it is wise to pass in a decimated detector set through offs_rep. The second step is to assign pixels to threads; to do this a signal-sized array is projected from map space. In the present implementation this uses all the detectors (but takes advantage of threads). In principle this step could be done with a decimated detector set, though with care that the decimated subset covers the array well, spatially.

Parameters:

• sight (CelestialSightLine) – The boresight pointing.

• fplane (FocalPlane) – The detector pointing

• shape (tuple) – The map shape.

• wcs (wcs) – The map WCS.

• tile_shape (tuple) – The tile shape, if this should be done using tiles.

• active_tiles (list) – The list of active tiles. If None, this will be computed along the way.

• n_threads (int) – The number of threads to target (defaults to OMP_NUM_THREADS).

• fplane_rep (FocalPlane) – A representative set of detectors, for determining scan direction (if not present then fplane is used for this).

• plot_prefix (str) – If present, pylab will be imported and some plots saved with this name as prefix.

Returns:

Ranges matrix with shape (n_thread, n_det, n_samp), that can be passes as the “threads” argument to Projectionist routines.

_ProjectionistBase

class so3g.proj.wcs._ProjectionistBase

This is a base class to assist with populating data structures needed for accelerated pointing routines. This class should not be used directly; call instead Projectionist for rectangular pixelizations or ProjectionistHealpix for HEALPix.

On instantiation, it carries information about the relation between celestial spherical coordinates and the Native Spherical coordinates of the projection, and also the pixelization scheme for the projection.

The following symbols are used in docstrings to represent shapes:

• n_dets - Number of detectors in the focal plane.

• n_samps - Number of samples in the signal space of each detector.

• n_threads - Number of threads to use in parallelized computation.

• n_mapidx - Number of indices required to specify a pixel in a map. For simple RectPix maps this is probably 2, with the first index specifying row and the second index specifying column. But for tiled maps n_mapidx is 3, with the first axis specifying the tile number. For HEALPix, n_mapidx is 1 for untiled and 2 for tiled maps.

• n_comp - Number of map components, such as Stokes components or other spin components. This takes the value 3, for example, if projection into TQU space has been requested.

When coordinate computation or projection routines are called, a ProjEng (a wrapped C++ object) is instantiated to perform the accelerated computations. The spin-composition (i.e. is this a ‘T’ map or a ‘TQU’ map) must be specified at this time, to instantiate the correct accelerated object.

Some method parameters are common to many functions and are documented here for consistency:

• assembly - an Assembly, providing a focal plane (quaternion offsets for each detector) as well as a boresight vector (quaternion for each time sample).

• comps - a string specifying the Stokes components of the map, for example ‘T’ or ‘TQU’. When absent, this will be guessed from the map shape; with 1|2|3 mapping to ‘T’|’QU’|’TQU’ respectively.

• threads - the thread assignment, consisting of a list of RangesMatrix objects. Each RangesMatrix object must have shape (n_threads, n_dets, n_samps). The n_threads does not need to be the same for every entry in the list. In TOD-to-map operations, each entry of this list is processed fully before proceeding to the next one. Each entry “ranges” is processed using (up to) the specified number of threads, such that thread i performs operations only on the samples included in ranges[i,:,:]. Most thread assignment routines in this module will return a list of two RangesMatrix objects, [ranges_parallel, ranges_serial]. The first item represents the part of the computation that can be done in parallel, and has shape (n_threads, n_dets, n_samps). The ranges_serial object has shape (1, n_dets, n_samps) and represents any samples that need to be treated in a single thread. The ranges_serial is only non-trivial when interpolation is active.

get_ProjEng(comps='TQU', proj_name=None, get=True, instance=True, interpol=None)

Returns an so3g.ProjEng object appropriate for use with the configured geometry.

get_pixels(assembly)

Get the pixel indices for the provided pointing Assembly. For each detector, an int32 array of shape (n_samps, n_mapidx) is returned. The value of the first slot of the second axis will be -1 if and only if the sample’s projected position is outside the map footprint.

See class documentation for description of standard arguments.

get_pointing_matrix(assembly)

Get the pointing matrix information, which is to say both the pixel indices and the “spin projection factors” for the provided pointing Assembly. Returns (pixel_indices, spin_proj) where pixel_indices is as returned by get_pixels() and spin_proj is a list of float arrays of shape [n_samps, n_comp]. As an alternative to on-the-fly computation, this information can be cached and used for very fast projection operations.

See class documentation for description of standard arguments.

get_coords(assembly, use_native=False, output=None)

Get the spherical coordinates for the provided pointing Assembly. For each detector, a float64 array of shape [n_samps,4] is returned. In the right-most index, the first two components are the longitude and latitude in radians. The next two components are the cosine and sine of the parallactic angle.

This routine uses a trivial CAR projection to return results from the celestial coordinate system (i.e. (lon,lat) = (RA,dec)), and the parallactic angle will be relative to that projection. If you are interested in the parallactic angle of the native spherical coordinate system (e.g. if you’re doing a tangent plane projection), make a second call specifying use_native=True. In this case you might also take a look at the get_planar() routine.

See class documentation for description of standard arguments.

get_planar(assembly, output=None)

Get projection plane coordinates for all detectors at all times. For each detector, a float64 array of shape [n_samps,4] is returned. The first two elements are the x and y projection plane coordiantes, similar to the “intermediate world coordinates”, in FITS language. Insofar as FITS ICW has units of degrees, these coordinates have units of radians. Indices 2 and 3 carry the cosine and sine of the detector parallactic angle.

zeros(super_shape)

Return a map, filled with zeros, with shape (super_shape,) + self.shape. For tiled maps, this will be a list of map tiles (with shape (super_shape, ) + tile_shape. If any tiles are not active, None is returned in the corresponding slots.

to_map(signal, assembly, output=None, det_weights=None, threads=None, comps=None)

Project signal into a map.

Parameters:

• signal (Signal-like) – The signal to project.

• output (Map-like) – The map into which to accumulate the projected signal. If None, a map will be initialized internally.

• det_weights (shape n_det) – Weights to use for each detector.

See class documentation for description of standard arguments.

to_weights(assembly, output=None, det_weights=None, threads=None, comps=None)

Project pointing into a weights map.

Parameters:

• output (Map-like) – The weights map into which to accumulate the projected signal. If None, a map will be initialized internally.

• det_weights (shape n_det) – Weights to use for each detector.

See class documentation for description of standard arguments.

from_map(src_map, assembly, signal=None, comps=None)

De-project from a map, returning a Signal-like object.

Parameters:

• src_map (Map-like) – The map from which to sample.

• signal (Signal-like) – The object into which to accumulate the signal. If not provided, a suitable object will be created and initialized to zero.

See class documentation for description of standard arguments.

assign_threads(assembly, method='domdir', n_threads=None)

Get a thread assignment RangesMatrix. Different algorithms can be chosen, though this method does not provide fine-grained control of algorithm parameters and instead seeks to apply reasonable defaults.

The methods exposed here are:

• 'simple': Divides the map into n_threads horizontal bands, and assigns the pixels in each band to a single thread. This tends to be bad for scans that aren’t horizontally correlated, or that don’t have equal weight throughout the map. The ‘domdir’ algorithm addresses both of those problems.

• 'domdir': For constant elevation scans, determines the dominant scan direction in the map and divides it into bands parallel to that scan direction. The thickness of bands is adjusted so roughly equal numbers of samples fall into each band.

• 'tiles': In Tiled projections, assign each tile to a thread. Some effort is made to balance the total number of samples over the threads.

The returned object can be passed to the threads argument of the projection methods in this class.

See class documentation for description of standard arguments.

Parameters:

method (str) – Algorithm to use.

assign_threads_from_map(assembly, tmap, n_threads=None)

Assign threads based on a map.

The returned object can be passed to the threads argument of the projection methods in this object.

See class documentation for description of standard arguments.

Parameters:

tmap – Map structure with shape (1,m,n) where each pixel, when converted to integer, species the index of the thread to which that pixel should be assigned.

get_active_tiles(assembly, assign=False)

For a tiled Projection, figure out what tiles are hit by an assembly.

See class documentation for description of standard arguments.

Parameters:

assign (bool or int) – If True or an int, then the function will also compute a thread assignment, based on assigning each tile to a particular thread. If this is an int, it sets the thread count; if it is simply True then OMP_NUM_THREADS is used.

Returns:

The basic analysis results are:

• 'active_tiles' (list of int): sorted, non-repeating list of tiles that are hit by the assembly.

• 'hit_counts' (list of int): the number of hits in each tile returned in ‘active_tiles’, respectively.

If the thread assignment took place then the dict will also contain:

• 'group_tiles' (list of lists of ints): There is one entry per thread, and the entry lists the tiles that have been assigned to that thread.

• 'group_ranges' (RangesMatrix): The thread assignments; shape (n_thread, n_det, n_samp).

Return type:

dict

Projectionist

class so3g.proj.Projectionist

This class assists with analyzing WCS information to populate data structures needed for accelerated pointing routines for rectangular pixelization.

See also the methods and parameters defined in the base class _ProjectionistBase.

As in pixell, the code and discussion here uses the term “geometry” to refer to the combination of an astropy.WCS object and a 2-d array shape.

Some common method parameters specific to rectangular pixelization are documented here for consistency.

• proj_name - a string specifying a projection. The nomenclature is mostly the same as the FITS CTYPE identifiers. Accepted values: ARC, CAR, CEA, TAN, ZEA, Flat, Quat.

• shape - the shape describing the celestial axes of the map. (For tiled representation this specifies the parent footprint.)

• wcs - the WCS describing the celestial axes of the map. Together with shape this is a geometry; see pixell.enmap documentation.

• interpol: How positions that fall between pixel centers will be handled. Options are “nearest” (default): Use Nearest Neighbor interpolation, so a sample takes the value of whatever pixel is closest; or “bilinear”: linearly interpolate between the four closest pixels. bilinear is slower (around 60%) but avoids problems caused by a discontinuous model.

naxis

2-element integer array specifying the map shape (for the 2 celestial axes).

cdelt

2-element float array specifying the pixel pitch.

crpix

2-element float array specifying the pixel coordinates of the reference point.

proj_name

string, name of the projection.

q_celestial_to_native

quaternion rotation taking celestial coordinates to the native spherical coordinates of the projection.

tile_shape

2-element integer array specifying the tile shape (this is the shape of one full-size sub-tile, not the decimation factor on each axis).

tiling

a _Tiling object if the map is tiled, None otherwise.

static get_q(wcs)

Analyze a wcs object to compute the quaternion rotation from celestial to native spherical coordinates.

classmethod for_geom(shape, wcs, interpol=None)

Construct a Projectionist for use with the specified “geometry”.

The shape and wcs are the core information required to prepare a Projectionist, so this method is likely to be called by other constructors.

classmethod for_map(emap, wcs=None, interpol=None)

Construct a Projectionist for use with maps having the same geometry as the provided enmap.

Parameters:

• emap – enmap from which to extract shape and wcs information. It is acceptable to pass a bare ndarray here (or anything with shape attribute), provided that wcs is provided separately.

• wcs – optional WCS object to use instead of emap.wcs.

classmethod for_source_at(alpha0, delta0, gamma0=0.0, proj_name='TAN')

Return a pointing-only Projectionist where some particular equatorial position will be put at the North Pole of the Native spherical coordinates.

classmethod for_tiled(shape, wcs, tile_shape, active_tiles=True, interpol=None)

Construct a Projectionist for use with the specified geometry (shape, wcs), cut into tiles of shape tile_shape.

See class documentation for description of standard arguments.

Parameters:

• tile_shape – tuple of ints, giving the shape of each tile.

• active_tiles – bool or list of ints. Specifies which tiles should be considered active. If True, all tiles are populated. If None or False, this will remain uninitialized and attempts to generate blank maps (such as calls to zeros or to projection functions without a target map set) will fail. Otherwise this must be a list of integers specifying which tiles to populate on such calls.

from_map(src_map, assembly, signal=None, comps=None)

De-project from a map, returning a Signal-like object. See parent class documentation for full description.

ProjectionistHealpix

class so3g.proj.ProjectionistHealpix

Projectionist for Healpix maps. See base class _ProjectionistBase for more methods and explanation of common method parameters.

nside

int, nside of the map, power of 2; 0 < nside <= 8192.

nside_tile

None, int, or ‘auto’, nside of tiling. Ntile will be 12*nside_tile**2. None for untiled. If ‘auto’, an appropriate nside_tile will be computed and set on calling compute_nside_tile().

ordering

str, ‘NEST’ or ‘RING’. Only NEST supported for tiled maps.

classmethod for_healpix(nside, nside_tile=None, active_tiles=None, ordering='NEST', interpol=None)

Construct a Projectionist for Healpix maps.

See class documentation for description of standard arguments.

compute_nside_tile(assembly, nActivePerThread=5, nThreads=None)

Automatically compute and set self.nside_tile for good balancing over threads.

Parameters:

• nActivePerThread – int, how many active pixels do you want per thread (minimum)

• nThreads – int, number of threads to optimize for (takes from OMP_NUM_THREADS if None)

get_active_tiles(assembly, assign=False)

For a tiled Projection, figure out what tiles are hit by an assembly. See parent class documentation for full description.

get_coords(assembly, use_native=False, output=None)

Get the spherical coordinates for the provided pointing Assembly. See parent class documentation for full description.

from_map(src_map, assembly, signal=None, comps=None)

De-project from a map, returning a Signal-like object. See parent class documentation for full description.

assign_threads(assembly, method=None, n_threads=None)

Get a thread assignment RangesMatrix. Available methods are 'simple' and 'tiles'. See parent class documentation for full description.

quat

so3g.proj.quat.euler(axis, angle)

The quaternion representing of an Euler rotation.

For example, if axis=2 the computed quaternion(s) will have components:

q = (cos(angle/2), 0, 0, sin(angle/2))

Parameters:

• axis ({0, 1, 2}) – The index of the cartesian axis of the rotation (x, y, z).

• angle (float or 1-d float array) – Angle of rotation, in radians.

Return type:

quat or G3VectorQuat, depending on ndim(angle).

so3g.proj.quat.rotation_iso(theta, phi, psi=None)

Returns the quaternion that composes the Euler rotations:

Qz(phi) Qy(theta) Qz(psi)

Note arguments are in radians.

so3g.proj.quat.rotation_lonlat(lon, lat, psi=0.0)

Returns the quaternion that composes the Euler rotations:

Qz(lon) Qy(pi/2 - lat) Qz(psi)

Note arguments are in radians.

so3g.proj.quat.rotation_xieta(xi, eta, gamma=0)

Returns the quaternion that rotates the center of focal plane to (xi, eta, gamma). This is equivalent to composed Euler rotations:

Qz(phi) Qy(theta) Qz (psi)

where

xi = - sin(theta) * sin(phi) eta = - sin(theta) * cos(phi) gamma = psi + phi

Note arguments are in radians.

so3g.proj.quat.decompose_iso(q)

Decomposes the rotation encoded by q into the product of Euler rotations:

q = Qz(phi) Qy(theta) Qz(psi)

Parameters:

q (quat or G3VectorQuat) – The quaternion(s) to be decomposed.

Returns:

(theta, phi, psi) – The rotation angles, in radians.

Return type:

tuple of floats or of 1-d arrays

so3g.proj.quat.decompose_lonlat(q)

Like decompose_iso, but returns (lon, lat, psi) assuming that the input quaternion(s) were constructed as rotation_lonlat(lon, lat, psi).

so3g.proj.quat.decompose_xieta(q)

Like decompose_iso, but returns (xi, eta, gamma) assuming that the input quaternion(s) were constructed as rotation_xieta(xi, eta, gamma).

Ranges

so3g.proj.Ranges

alias of RangesInt32

class so3g.RangesInt32

A finite series of non-overlapping semi-open intervals on a domain of type: int32_t.

To create an empty object, instantiate with just a sample count: RangesInt32(count).

Alternately, consider convenience methods such as from_mask, from_array, and from_bitmask; see below.

In addition to the methods explained below, note the that following operators have been defined and perform as follows (where r1 and r2 are objects of this class:

• ~r1 is equivalent to r1.complement()

• r1 *= r2 is equivalent to r1.intersect(r2)

• r1 += r2 is equivalent to r1.merge(r2)

• r1 * r2 and r1 + r2 behave as you might expect, returning a new object and leaving r1 and r2 unmodified.

The object also supports slicing. For example, if r1 has count = 100 then r1[10:-5] returns a new object (not a view) that has count = 85. A data member reference keeps track of the history of shifts in the first sample; for example if r1.reference = 0 then r1[10:-5].reference will be -10. This variable can be interpreted as giving the logical index, in the new index system, of where index=0 of the original object would be found. This is useful for bookkeeping in some cases.

add_interval(start, end) → RangesInt32

Merge an interval into the set.

append_interval_no_check(start, end) → RangesInt32

Append an interval to the set without checking for overlap or sequence.

merge(src) → RangesInt32

Merge ranges from another RangesInt32 into this one.

buffer(buff) → RangesInt32

Buffer each interval by an amount specified by buff

buffered(buff) → RangesInt32

Return an interval buffered by buff

close_gaps(gap) → RangesInt32

Remove gaps between ranges less than gap

intersect(src) → RangesInt32

Intersect another RangesInt32 with this one.

complement() → RangesInt32

Return the complement (over domain).

zeros_like() → RangesInt32

Return range of same length but no intervals

ones_like() → RangesInt32

Return range of same length and interval spanning count

ranges() → object

Return the intervals as a 2-d numpy array of ranges.

static from_array(data, count) → RangesInt32

The input data must be an (n,2) shape ndarray of int32. The integer count sets the domain of the object.

static from_bitmask(bitmask_array) → object

Return a list of RangesInt32 extracted from an ndarray encoding a bitmask.

static from_mask(bool_array) → object

Return a list of RangesInt32 extracted from an ndarray of bool.

static bitmask(ranges_list, n_bits) → object

Return an ndarray bitmask from a list ofRangesInt32. n_bits determines the output integer type. Bits are assigned from LSB onwards; use None in the list to skip a bit.

mask() → object

Return a boolean mask from this Ranges object.

copy() → RangesInt32

Get a new object with a copy of the data.

RangesMatrix

class so3g.proj.RangesMatrix(items=[], child_shape=None)

This is a wrapper for multi-dimensional grids of Ranges objects. This can be used to store Ranges objects per-detector (the 2-d case) or per-thread and per-detector (a 3-d case, required if using OpenMP). The right-most axis always corresponds to the sample axis carried by the Ranges objects.

In addition to .shape and multi-dimensional slicing, it supports inversion (the ~ operator) and multiplication against other RangesMatrix or Ranges objects, with broadcasting rules similar to standard numpy array rules.

static concatenate(items, axis=0)

Static method to concatenate multiple RangesMatrix (or Ranges) objects along the specified axis.

close_gaps(gap_size=0)

Call close_gaps(gap_size) on all children. Any ranges that abutt each other within the gap_size are merged into a single entry. Usually a gap_size of 0 is not possible, but if a caller is carefully using append_interval_no_check, then it can happen.

static full(shape, fill_value)

Construct a RangesMatrix with the specified shape, initialized to fill_value.

Parameters:

• shape (tuple of int) – The shape. Must have at least one element. If there is only one element, a Ranges object is returned, not a RangesMatrix.

• fill_value (bool) – True or False. If False, the wrapped Ranges objects will have no intervals. If True, the wrapped Ranges objects will all have a single interval spanning their entire domain.

Returns:

Ranges object (if shape has a single element) or a RangesMatrix object (len(shape) > 1).

See also: zeros, ones.

classmethod zeros(shape)

Equivalent to full(shape, False).

classmethod ones(shape)

Equivalent to full(shape, True).

classmethod from_mask(mask)

Create a RangesMatrix from a boolean mask. The input mask can have any dimensionality greater than 0 but be aware that if ndim==1 then a Ranges object, and not a RangesMatrix, is returned.

Parameters:

mask (ndarray) – Must be boolean array with at least 1 dimension.

Returns:

RangesMatrix with the same shape as mask, with ranges corresponding to the intervals where mask was True.

mask(dest=None)

Return the boolean mask equivalent of this object.

Weather

class so3g.proj.Weather(*args, **kwargs)

This is a thin wrapper around a dict or some other convenient container, to simplify passing telescope atmospheric weather information to various consumers.

Important keys:

• ‘temperature’: ambient temperature, in degrees Celsius.

• ‘pressure’: air pressure, in mBar.

• ‘humidity’: relative humidity, as fraction of 1.

If you find that some consumer prefers some other units (K and atm, for example…), then define a function that makes the necessary conversions.

apply(something)

Helper method to load configure an object with this weather info. Behavior depends on the type of observer. The types below will be checked, in order.

• ephem.Observer – sets pressure and temperature.

to_qpoint()

Extract a dict appropriate for passing to qpoint functions.