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 25 commits into
base: main
Choose a base branch
from
Open
Changes from all commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
b13d14f
Healpix doc
Vidushi-GitHub Jul 29, 2024
7d88cc8
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 18, 2024
e87d9a4
remove repetition
Vidushi-GitHub Dec 18, 2024
b3f7b25
add links
Vidushi-GitHub Dec 18, 2024
b4056be
alternatives link
Vidushi-GitHub Dec 18, 2024
9525aa0
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
eb52937
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
3719156
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
559aa88
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
7f32700
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
e9a1340
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
b9f1b53
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
d95568f
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
4bea18f
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
1dc841a
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
a001d16
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
85a616e
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
3578afa
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
7c252ef
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
acd6b55
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
ed50953
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
439556a
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
fa6f9b8
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
2276168
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
8d0bd2c
Update app/routes/docs.client.samples.md
Vidushi-GitHub Dec 22, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
179 changes: 179 additions & 0 deletions app/routes/docs.client.samples.md
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,9 @@ and some samples from the FAQs section of the [gcn-kafka-python](https://github.

To contribute your own ideas, make a GitHub pull request to add it to [the Markdown source for this document](https://github.com/nasa-gcn/gcn.nasa.gov/blob/CodeSamples/app/routes/docs.client.samples.md), or [contact us](/contact).

- [Working with Kafka messages](#parsing)
- [HEALPix Sky Maps](#healpix-sky-maps)

## Parsing

Within your consumer loop, use the following functions to convert the
Expand Down Expand Up @@ -164,3 +167,179 @@ for message in consumer.consume(end[0].offset - start[0].offset, timeout=1):
continue
print(message.value())
```

## HEALPix Sky Maps

[HEALPix](https://healpix.sourceforge.io) (<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.

### Multi-Order Sky Maps

GCN strongly encourages the use of multi-order sky maps. They 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.

#### 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 can be used for this purpose; a number of [alternatives](#references) are listed at the bottom of this page)

```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.multiorder.fits')
```

#### Most Probable Sky Location

The index of the highest probability point

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

level, ipix = ah.uniq_to_level_ipix(uniq)
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 = ah.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)
```

#### 90% Probability Region

The estimation of a 90% probability region involves sorting the pixels, calculating the area of each pixel, and then summing the probability of each pixel until 90% is reached.

```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)
```

### Flat Resolution HEALPix Sky maps

Let's say you have a sky map fits file. To read in Python with Healpy:

```python
import healpy as hp
import numpy as np

# Read both the HEALPix image data and the FITS header
hpx, header = hp.read_map('skymap.fits.gz', h=True)

# Plot a Mollweide-projection all-sky image
np.mollview(hpx)
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
np.mollview(hpx)
hp.mollview(hpx)

Copy link
Member

Choose a reason for hiding this comment

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

Also, if you actually want the plot to appear, then you have to add:

from matplotlib import pyplot as plt

plt.show()

```

#### Most Probable Sky Location

```python
# Find the highest probability pixel
ipix_max = np.argmax(hpx)

# Probability density per square degree at that position
hpx[ipix_max] / hp.nside2pixarea(nside, degrees=True)
Copy link
Member

Choose a reason for hiding this comment

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

nside has not been defined.


# Highest probability pixel on the sky
ra, dec = hp.pix2ang(nside, ipix_max, lonlat=True)
ra, dec
```

#### Integrated probability in a Circle

We call hp.query_disc to obtain an array of indices for the pixels inside the circle.

```python
# First define the Cartesian coordinates of the center of the circle
ra = 180.0
dec = -45.0
radius = 2.5

# Calculate Cartesian coordinates of the center of the Circle
xyz = hp.ang2vec(ra, dec, lonlat=True)

# Obtain an array of indices for the pixels inside the circle
ipix_disc = hp.query_disc(nside, xyz, np.deg2rad(radius))

# Sum the probability in all of the matching pixels:
hpx[ipix_disc].sum()
```

#### Integrated probability in a Polygon

We can use the `hp.query_polygon` function to find the pixels inside a polygon and compute the probability that the source is inside the polygon by adding up the pixel values.

```python
# Indices of the pixels within a polygon (defined by the Cartesian coordinates of its vertices)

xyz = [[-0.5, -0.4, -0.5],
[-0.4, -0.4, -0.6],
[-0.6, -0.3, -0.6],
[-0.7 , -0.4, -0.5]]
Comment on lines +311 to +314
Copy link
Member

Choose a reason for hiding this comment

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

These are not unit vectors.

ipix_poly = hp.query_polygon(nside, xyz)
hpx[ipix_poly].sum()
```

##### Other Documentation and HEALPix Packages

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

[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/): Integrates HEALPix with Astropy for data manipulation and analysis

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

## References

(1) Calabretta, M. R., & Roukema, B. F. 2007, Mon. Notices Royal Astron. Soc., 381, 865. [doi: 10.1111/j.1365-2966.2007.12297.x](https://doi.org/10.1111/j.1365-2966.2007.12297.x)

(2) Górski, K.M., Hivon, E., Banday, A.J., et al. 2005, Astrophys. J., 622, 759. [doi: 10.1086/427976](https://doi.org/10.1086/427976)

(3) Górski, K. M., Wandelt, B. D., et al. 1999. [doi: 10.48550/arXiv.astro-ph/9905275](https://doi.org/10.48550/arXiv.astro-ph/9905275)

(4) Fernique, P., Allen, et al. 2015, Astron. Astrophys., 578, A114. [doi: 10.1051/0004-6361/201526075](https://doi.org/10.1051/0004-6361/201526075)

(5) Fernique, P., Boch, T., et al. 2014, IVOA Recommendation. [doi: 10.48550/arXiv.1505.02937](https://doi.org/10.48550/arXiv.1505.02937)

(6) Martinez-Castellanos, I., Singer, L. P., et al. 2022, Astrophys. J., 163, 259. [doi: 10.3847/1538-3881/ac6260](https://doi.org/10.3847/1538-3881/ac6260)

(7) Singer, L. P., & Price, L. R. 2016, Phys. Rev. D, 93, 024013. [doi: 10.1103/PhysRevD.93.024013](https://doi.org/10.1103/PhysRevD.93.024013)
Loading