Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Healpix doc #2498

Open
wants to merge 13 commits into
base: main
Choose a base branch
from
110 changes: 110 additions & 0 deletions app/routes/docs.notices.healpix.mdx
Original file line number Diff line number Diff line change
@@ -0,0 +1,110 @@
---
handle:
breadcrumb: HEALPix Maps
---

import { Highlight } from '~/components/Highlight'
import { Tab, Tabs } from '~/components/tabs'

# HEALPix Sky Maps

HEALPix stands for <b>H</b>ierarchical, <b>E</b>qual <b>A</b>rea, and iso-<b>L</b>atitude <b>Pix</b>elisation is a scheme for indexing positions on the unit sphere.
For localization of events, the multi-messenger community uses the standard [HEALPix](https://healpix.sourceforge.io) format with the file extension `.fits.gz`, as well as multi-resolution HEALPix format with the file extension `.multiorder.fits`. The preferred format is the multi-resolution HEALPix format, that is only format explicitly listed in the GCN.

## Multi-Order Sky Maps

GCN strongly encourages the use of multi-order sky maps. These are a more complex form of healpix sky maps that utilize a variable resolution, with higher probability regions having higher resolution and lower probability regions being encoded with a lower resolution. This is signicantly more efficient than single-resolution healpix sky maps with respect to storage footprint and read speed. However, interpreting these multi-order sky maps can be more complex.

## Working with HEALPix Maps

HEALPix data is distrubuted in standard format with the file extension `.fits.gz`. These files are FITS image files and can be utilized with many common FITS tools. These files usually stored as HEALPix projection,

### Reading Sky Maps

Sky maps can be parsed using python; to start, import a handful of packages (note: while this documentation covers the use of `astropy-healpix`, there are several packages that are ab)
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved
tylerbarna marked this conversation as resolved.
Show resolved Hide resolved

```python
import astropy_healpix as ah
import numpy as np

from astropy import units as u
from astropy.table import QTable
```

A given sky map can then be read in as

```python
skymap = QTable.read('skymap.multiofits)
```

### Most Probable Sky Location

The index of the highest probability point can be found by doing the following

```python
hp_index = np.argmax(skymap['PROBDENSITY'])
uniq = skymap[hq_index]['UNIQ']

level, ipix = ah.uniq_to_level_ipix(uniq)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should introduce a UNIQ scheme a bit.

nside = ah.level_to_nside(level)

ra, dec = ah.healpix_to_lonlat(ipix, nside, order='nested')
```

### Probability Density at a Known Position

Similarly, one can calculate the probability density at a known position

```python
ra, dec = 197.450341598 * u.deg, -23.3814675445 * u.deg

level, ipix = ah.uniq_to_level_ipix(skymap['UNIQ'])
nside = ah.level_to_nside(level)

match_ipix = ad.lonlat_to_healpix(ra, dec, nside, order='nested')

match_index = np.flatnonzero(ipix == match_ipix)[0]

prob_density = skymap[match_index]['PROBDENSITY'].to_value(u.deg**-2)
```

### Find the 90% Probability Region

The estimation of a 90% probability region, involves sorting the pixel, calculating the area of each pixel, and then proceeding by the probability of each pixel.

```python
#Sort the pixels by decending probability density
skymap.sort('PROBDENSITY', reverse=True)

#Area of each pixel
level, ipix = ah.uniq_to_level_ipix(skymap['UNIQ'])
pixel_area = ah.nside_to_pixel_area(ah.level_to_nside(level))

#Pixel area times the probability
prob = pixel_area * skymap['PROBDENSITY']

#Cummulative sum of probability
cumprob = np.cumsum(prob)

#Pixels for which cummulative is 0.9
i = cumprob.searchsorted(0.9)

#Sum of the areas of the pixels up to that one
area_90 = pixel_area[:i].sum()
area_90.to_value(u.deg**2)

```

## Further documentation

Additional information can be found the [LIGO website](https://emfollow.docs.ligo.org/userguide/tutorial/multiorder_sky maps.html)

## Libraries for reference

[healpy](https://healpy.readthedocs.io/en/latest/): Official python library for handling the pixlated data on sphere

[astropy-healpix](https://pypi.org/project/astropy-healpix/)

[mhealpy](https://mhealpy.readthedocs.io/en/latest/): Object-oriented wrapper of healpy for handling the multi-resolution maps

[MOCpy](https://cds-astro.github.io/mocpy/): Python library allowing easy creation, parsing and manipulation of Multi-Order Coverage maps.
3 changes: 3 additions & 0 deletions app/routes/docs.tsx
Original file line number Diff line number Diff line change
Expand Up @@ -177,6 +177,9 @@ export default function () {
<NavLink key="schema" to="notices/schema">
Unified Schema
</NavLink>,
<NavLink key="healpix" to="notices/healpix">
HEALPix
</NavLink>,
<NavLink key="archive" to="notices/archive">
Archive
</NavLink>,
Expand Down
Loading