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 5 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
191 changes: 191 additions & 0 deletions app/routes/docs.client.samples.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,9 @@ handle:
breadcrumb: Sample Code
---

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

Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved
# Sample Code

Here is a collection of functions and example code that may be useful
Expand All @@ -12,6 +15,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 +170,188 @@ for message in consumer.consume(end[0].offset - start[0].offset, timeout=1):
continue
print(message.value())
```

## HEALPix Sky Maps

HEALPix (<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.
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved
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.
jracusin marked this conversation as resolved.
Show resolved Hide resolved
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved
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
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.
For localization of events, the multi-messenger community uses the standard [HEALPix-in-FITS](https://healpix.sourceforge.io/data/examples/healpix_fits_specs.pdf) 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. This is generally the format that is used in GCN Notices.


### 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)
Vidushi-GitHub 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)
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved
```

#### Most Probable Sky Location

The index of the highest probability point

```python
hp_index = np.argmax(skymap['PROBDENSITY'])
uniq = skymap[hq_index]['UNIQ']
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved

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
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved

```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')
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved

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 pixel, calculating the area of each pixel, sorting by probability, and then summing the probability of each pixel until 90% is reached.
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved

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

Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved
```

#### Flat Resolution HEALPix Sky maps
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved

Let's say you have a sky map fits file, read in Python with Healpy:
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
Let's say you have a sky map fits file, read in Python with Healpy:
Let's say you have a sky map fits file. To read in Python with Healpy:

Also, add some transition text to make it clear that you have switched here to flat-resolution files.


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

# Read both the HEALPix image data and the FITS header
hpx, header = hp.read_map('*.fits.gz,0', h=True)
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved

# Plotting a Mollweide-projection all-sky map
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved
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
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved

```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
theta, phi = hp.pix2ang(nside, ipix_max)
ra = np.rad2deg(phi)
dec = np.rad2deg(0.5 * np.pi - theta)
ra, dec
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved
```

#### 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
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved
ra = 180.0
dec = -45.0
radius = 2.5

# We convert to spherical polar coordinates
theta = 0.5 * np.pi - np.deg2rad(dec)
phi = np.deg2rad(ra)
radius = np.deg2rad(radius)

# Calculate Cartesian coordinates of the center of the Circle
xyz = hp.ang2vec(theta, phi)
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved

# Obtain an array of indices for the pixels inside the circle
ipix_disc = hp.query_disc(nside, xyz, radius)
Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved

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

Vidushi-GitHub marked this conversation as resolved.
Show resolved Hide resolved
```

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