Healpix doc

app/routes/docs.client.samples.md

@@ -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
@@ -162,3 +165,188 @@ 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
+from matplotlib import pyplot as plt
+
+# 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)
+plt.show()
+```
+
+#### Most Probable Sky Location
+
+```python
+# Reading Sky Maps with Healpy
+healpix_image = hp.read_map('bayestar.fits.gz,0')
+npix = len(hpx)
+
+# Lateral resolution of the HEALPix map
+nside = hp.npix2nside(npix)
+
+# 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)
+
+# 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, 0, 0],
+       [1, 0, 0],
+       [1, 1, 0],
+       [0, 1, 0]]
+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_skymaps.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)