The HealPix Module (gdt.core.healpix)

Introduction

HEALPix, or Hierarchical Equal Area isoLatitude Pixelation of a sphere, is an algorithm for dividing up the sky into equal-area pixels. This is incredibly useful for analyzing data that cover the entire sky or large patches of sky. In the GDT, there are two main purposes supported: for storing, computing, and plotting the all-sky effective area (sensitivity) of a gamma-ray instrument; and for storing, calculating, plotting, and functionality for large-area localizations on the sky. Other purposes can be developed using the HealPix base class from which the other two purposes are derived, but we will focus on examples serving effective area and localization.

HEALPix for Effective Area

The HealPixEffectiveArea class provides a way to create, store, and read files that contain effective area information. HEALPix is a natural algorithm for doing this when the effective area of an instrument is over a large region of the sky (or the full sky). While this class can be subclassed, there are some generic capabilities available. For example, we can create an effective area object with a simple cosine angular response:

>>> from gdt.core.healpix import HealPixEffectiveArea
>>> # cosine centered at spacecraft Az=60, Zen=30 with 100 cm^2 peak eff. area
>>> eff_area1 = HealPixEffectiveArea.from_cosine(60.0, 30., 100.0)
>>> eff_area1
<HealPixEffectiveArea: NSIDE=64>

By default, this created a HEALPix map with an NSIDE resolution of 64. While this is a standard cosine response, the cosine coefficient can be adjusted to make the angular response wider or narrower by specifying the coeff keyword.

We can also create a (trivial) uniform response:

>>> # uniform effective area of 50 cm^2.
>>> eff_area2 = HealPixEffectiveArea.from_uniform(50.0)

And the effective areas can be summed:

>>> eff_area_sum = HealPixEffectiveArea.sum(eff_area1, eff_area2, output_nside=64)

You can directly access the HEALPix array via the eff_area attribute:

>>> eff_area_sum.eff_area
array([137.21163256, 136.76058684, 135.97935273, ...,  50.        ,
        50.        ,  50.        ])

Or you can calculate the effective area at any arbitrary point:

>>> # effective area just off the normal
>>> eff_area_sum.effective_area(61.0, 31.0)
149.97909654227354

>>> # effective area way off the normal
>>> eff_area_sum.effective_area(240.0, 120.0)
50.0

Primarily used for plotting purposes, we can have the HEALPix array map to a grid of Azimuth and Zenith:

>>> eff_area_grid, az_grid, zen_grid = eff_area_sum.make_grid()
>>> eff_area_grid.shape
(180, 360)
>>> zen_grid.shape, az_grid.shape
((180,), (360,))

There are a default number of grid points along the Azimuth and Zenith directions that can be changed by specifying the numpts_az and numpts_zen keywords. This function returns a grid of effective area values that has the same dimensions as a grid in the zenith direction and the azimuth direction. The GDT uses these arrays with matplotlib.pyplot.pcolormesh to make the effective area plot.

HEALPix for Localizations

Similar to the HEALPix class for effective area, the HealPixLocalization class provides a way to create, store, analyze, and read localizations files. There are a few different ways to create a HealPixLocalization object. For example, a generic Gaussian localization can be generated by specifying the RA and Dec of the centroid of the localization and the Gaussian standard deviation:

>>> from gdt.core.healpix import HealPixLocalization
>>> # Gaussian centered at RA=50.0, Dec=-30.0, with 5 deg std. dev.
>>> gauss_loc = HealPixLocalization.from_gaussian(50.0, -30.0, 5.0)
>>> gauss_loc
<HealPixLocalization:
 NSIDE=128; trigtime=None;
 centroid=(50.2734375, -30.000000000000018)>

There are a few things to note here. First, there are optional arguments we didn’t set, such as a trigger time (which is why trigtime=None). Secondly, the centroid is not precisely what we defined because the HEALPix resolution is finite and it is telling us the pixel that contains the centroid. Finally, the NSIDE resolution is automatically determined by the the width of the Gaussian (smaller widths -> higher resolutions), but that can be overridden by specifying the nside keyword on creation.

There are a number of attributes available:

>>> # the number of HEALPix pixels
>>> gauss_loc.npix
196608

>>> # the HEALPix resolution
>>> gauss_loc.nside
128

>>> # the area of each pixel in square degrees
>>> gauss_loc.pixel_area
0.20982341130279172

>>> # the localization centroid
>>> gauss_loc.centroid
(50.2734375, -30.000000000000018)

You can also retrieve the localization probability density and confidence level at a given point:

>>> # probability density at RA=49.0, Dec=-29.0
>>> gauss_loc.probability(49.0, -29.0)
0.006086688085425917

>>> # confidence level at RA=49.0, Dec=-29.0
>>> gauss_loc.confidence(49.0, -29.0)
0.04367065341721976

The area of a confidence region can be computed in (square degrees):

>>> # 90% confidence region
>>> gauss_loc.area(0.9)
360.0569737955906

And similar to the HealPixEffectiveArea.make_grid() function, the localization can be mapped to a grid that is useful for plotting:

>>> prob_grid, ra_grid, dec_grid = gauss_loc.prob_array()
>>> prob_grid.shape
(180, 360)
>>> dec_grid.shape, ra_grid.shape
((180,), (360,))

Similar to the generating a Gaussian localization, we can generate a Gaussian-width annulus on the sky, such as one that results from measuring the difference in the time-of-arrival of photons detected by two different instruments separated by some distance.

>>> annulus_loc = HealPixLocalization.from_annulus(300.0, -30, 80.0, 3.0)
>>> annulus_loc
<HealPixLocalization:
 NSIDE=128; trigtime=None;
 centroid=(52.734375, 24.296477994805652)>

In this case, the we generated an annulus that has a center at RA=300, Dec=-30, an angular radius of 80 degrees, and a Gaussian standard deviation width of 3 degrees. Note that since this is a Gaussian-width annulus, the centroid is meaningless because the localization probability density is approximately a constant along the center of the annulus.

Localizations can be multiplied if, for example, the localizations are believed to be of the same source.

>>> mult_loc = HealPixLocalization.multiply(gauss_loc, annulus_loc)
<HealPixLocalization:
 NSIDE=128; trigtime=None;
 centroid=(59.765625, 27.615881983844705)>

We can also calculate the association probability between our localization and a known source on the sky. For example, if we have a known source that appears to be near our localizations:

>>> # known source at RA=55.0, Dec=25.0
>>> mult_loc.source_probability(55.0, 25.0)
0.9952878990151246

which is nearly a 3-sigma association. Alternatively, if we did not know for certain that our two localizations are from the same source, we can calculate the probability that the two localizations are spatially associated:

>>> gauss_loc.region_probability(annulus_loc)
1.6169670432806807e-07

Another way to create a HealPixLocalization object is by providing a list of vertices to create a localization polygon with uniform probability. For example, we can create a localization represented by a triangle:

>>> ra_pts = [270.0, 180.0, 90.0]
>>> dec_pts = [15.0, 60.0, 15.0]
>>> poly_loc = HealPixLocalization.from_vertices(ra_pts, dec_pts)
>>> poly_loc
<HealPixLocalization
 NSIDE=64; trigtime=None;
 centroid=(178.92857142857144, 58.91977535280316)>

Again, the centroid value for this localization may not make sense because the probability is uniformly distributed. Finally, a kernel can be convolved with the localization as long as the kernel can be expressed as a combination of Gaussians. For example, we can define a kernel as follows:

>>> import numpy as np
>>>
>>> def kernel_example(sigma_deg1, sigma_deg2):
>>>     """A 50/50 mixture of two Gaussians with different widths"""
>>>     sigma1 = np.deg2rad(sigma_deg1)
>>>     sigma2 = np.deg2rad(sigma_deg2)
>>>     frac1 = 0.5
>>>     return ([sigma1, sigma2], [frac1])

This example allows a user to specify the width of two different Gaussians that each represent 50% of the model. This can be used in the following way to produce a new HealPixLocalization object that is the result of the original localization convolved with the kernel:

>>> convolved_loc = poly_loc.convolve(kernel_example, 1.0, 3.0)

Reference/API

gdt.core.healpix Module

Classes

HealPix()	Base class for HEALPix files
HealPixEffectiveArea()	Class for effective area HEALPix files
HealPixLocalization()	Class for localization HEALPix files

Class Inheritance Diagram