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:

names = ['a', 'b', 'c']
x = np.array([-0.5, 0., 0.5]) * DEG
y = np.zeros(3)
gamma = np.array([0,30,60]) * DEG
fp = so3g.proj.FocalPlane.from_xieta(names, x, y, gamma)

This particular function, from_xieta(), will apply the SO standard coordinate definitions and store a rotation quaternion for each detector. FocalPlane is just a thinly wrapped OrderedDict, where the detector name is the key and the value is the rotation quaternion:

>>> fp['c']
spt3g.core.quat(0.866017,0.00377878,-0.00218168,0.499995)

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

# Get vector of quaternion pointings for detector 'a'
q_total = csl.Q * fp['a']
# 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)
[13.06731917] [-53.12633433]

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)
OrderedDict([('a', array([[ 0.22806774, -0.92722945, -0.9999468 ,
0.01031487]])), ('b', array([[ 0.24261138, -0.92726871, -0.86536489,
-0.5011423 ]])), ('c', array([[ 0.25715457, -0.92720643, -0.48874018,
-0.87242939]]))])

To be clear, coords() now returns a dictionary whose keys are the detector names. Each value is an array 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([13.04], dtype=float32), array([13.88], dtype=float32), array([14.72], dtype=float32)]
>>> [x/DEG for x in pix_dec]
[array([-53.100002], dtype=float32), array([-53.100002], dtype=float32), array([-53.100002], 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([[ 45, 148]], dtype=int32), array([[ 45, 106]], dtype=int32),
array([[45, 64]], 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([45, 45, 45]), array([ 64, 106, 148]))
>>> map_out[map_out!=0]
array([100.,  10.,   1.])

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([200.,  20.,   2.])

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[:,45,106]
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[...,45,106]
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[:,45,106]
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)

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(keyed=False, 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, det_offsets)

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.

• det_offsets (FocalPlane) – The “offsets” of each detector relative to the boresight.

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)

Construct a SightLine from horizon coordinates. This uses high-precision pointing.

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(det_offsets=None, output=None)

Get the celestial coordinates of each detector at each time.

Parameters

• det_offset – A dict or list or array of detector offset tuples. If each tuple has 2 elements or 3 elements, the typles are interpreted as X,Y[,phi] coordinates in the conventional way. If 4 elements, it’s interpreted as a quaternion. If this argument is None, then the boresight pointing is returned.

• output – An optional structure for holding the results. For that to work, each element of output must support the buffer protocol.

Returns

If the det_offset was passed in as a dict, then a dict with the same keys is returned. Otherwise a list is returned. For each detector, the corresponding returned object is an array with shape (n_samp, 4) where the four components correspond to longitude, latitude, cos(psi), sin(psi).

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

This class collects the focal plane positions and orientations for multiple named detectors. The classmethods can be used to construct the object from some common input formats.

classmethod from_xieta(names, xi, eta, gamma=0)

Creates a FocalPlane object for a set of detector positions provided in xieta projection plane coordinates.

Parameters

• names – vector of detector names.

• xi – vector of xi positions, in radians.

• eta – vector of eta positions, in radians.

• gamma – vector or scalar detector orientation, in radians.

The (xi,eta) coordinates are Cartesian projection plane coordinates. When looking at the sky along the telescope boresight, xi 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.

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, offs, shape, wcs, tile_shape=None, active_tiles=None, n_threads=None, offs_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.

• offs (array of quaternions) – 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).

• offs_rep (array of quaternions) – A representative set of detectors, for determining scan direction (if not present then offs 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.

Projectionist

class so3g.proj.Projectionist

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

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.

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.

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 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.

• 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.

• 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.

• threads - the thread assignment, which is a RangesMatrix with shape (n_threads,n_dets,n_samps), used to specify which samples should be treated by each thread in TOD-to-map operations. Such objects should satisfy the condition that threads[x,j]*threads[y,j] is the empty Range for x != y; i.e. each detector-sample is assigned to at most one thread.

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)

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)

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)

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.

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

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

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.

Returns

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 so3g._libso3g_docstring_shells.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.

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.