aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--fg21sim/utils/healpix.py58
1 files changed, 36 insertions, 22 deletions
diff --git a/fg21sim/utils/healpix.py b/fg21sim/utils/healpix.py
index 2c5d9ba..fcdd848 100644
--- a/fg21sim/utils/healpix.py
+++ b/fg21sim/utils/healpix.py
@@ -40,7 +40,8 @@ logger = logging.getLogger(__name__)
def healpix2hpx(data, append_history=None, append_comment=None):
- """Reorganize the HEALPix data (1D array as FITS table) into 2D FITS
+ """
+ Reorganize the HEALPix data (1D array as FITS table) into 2D FITS
image in HPX coordinate system.
Parameters
@@ -83,7 +84,8 @@ def healpix2hpx(data, append_history=None, append_comment=None):
def hpx2healpix(data, append_history=None, append_comment=None):
- """Revert the reorganization and turn the 2D image in HPX format
+ """
+ Revert the reorganization and turn the 2D image in HPX format
back into HEALPix data as 1D array.
Parameters
@@ -118,13 +120,14 @@ def hpx2healpix(data, append_history=None, append_comment=None):
if ((hpx_header["CTYPE1"], hpx_header["CTYPE2"]) !=
("GLON-HPX", "GLAT-HPX")):
raise ValueError("only Galactic 'HPX' projection currently supported")
+
# Calculate Nside
nside = round(hpx_header["NAXIS1"] / 5)
nside2 = round(90 / np.sqrt(2) / hpx_header["CDELT2"])
if nside != nside2:
raise ValueError("Cannot determine the Nside value")
logger.info("Determined HEALPix Nside=%d" % nside)
- #
+
npix = hp.nside2npix(nside)
logger.info("Calculating the HPX indexes ...")
hpx_idx = _calc_hpx_indexes(nside).flatten()
@@ -141,7 +144,8 @@ def hpx2healpix(data, append_history=None, append_comment=None):
@nb.jit(nb.int64[:](nb.int64, nb.int64, nb.int64), nopython=True)
def _calc_hpx_row_idx(nside, facet, jmap):
- """Calculate the HEALPix indexes for one row of a facet.
+ """
+ Calculate the HEALPix indexes for one row of a facet.
NOTE
----
@@ -153,7 +157,7 @@ def _calc_hpx_row_idx(nside, facet, jmap):
"""
I0 = [1, 3, -3, -1, 0, 2, 4, -2, 1, 3, -3, -1]
J0 = [1, 1, 1, 1, 0, 0, 0, 0, -1, -1, -1, -1]
- #
+
n2side = 2 * nside
n8side = 8 * nside
nside1 = nside - 1
@@ -162,7 +166,7 @@ def _calc_hpx_row_idx(nside, facet, jmap):
# double-pixelization pixel coordinates of the center of the facet
i0 = nside * I0[facet]
j0 = nside * J0[facet]
- #
+
row_idx = np.zeros(nside, dtype=np.int64)
for imap in range(nside):
# (i, j) are 0-based, double-pixelization pixel coordinates.
@@ -174,7 +178,7 @@ def _calc_hpx_row_idx(nside, facet, jmap):
if i < 0:
i += n8side
i += 1
- #
+
if j > nside:
# north polar regime
if j == n2side:
@@ -211,7 +215,8 @@ def _calc_hpx_row_idx(nside, facet, jmap):
@nb.jit(nb.int64[:, :](nb.int64), nopython=True)
def _calc_hpx_indexes(nside):
- """Calculate HEALPix element indexes for the HPX projection scheme.
+ """
+ Calculate HEALPix element indexes for the HPX projection scheme.
Parameters
----------
@@ -234,6 +239,7 @@ def _calc_hpx_indexes(nside):
"""
# number of horizontal/vertical facet
nfacet = 5
+
# Facets layout of the HPX projection scheme.
# Note that this appears to be upside-down, and the blank facets
# are marked with "-1".
@@ -248,10 +254,10 @@ def _calc_hpx_indexes(nside):
FACETS_LAYOUT[2, :] = [-1, 0, 4, 11, -1]
FACETS_LAYOUT[3, :] = [-1, -1, 3, 7, 10]
FACETS_LAYOUT[4, :] = [-1, -1, -1, 2, 6]
- #
+
shape = (nfacet*nside, nfacet*nside)
indexes = -np.ones(shape, dtype=np.int64)
- #
+
# Loop vertically facet-by-facet
for jfacet in range(nfacet):
# Loop row-by-row
@@ -267,12 +273,14 @@ def _calc_hpx_indexes(nside):
idx = _calc_hpx_row_idx(nside, facet, j)
col = ifacet * nside
indexes[row, col:(col+nside)] = idx
- #
+
return indexes
def _make_hpx_header(header, append_history=None, append_comment=None):
- """Make the FITS header for the HPX image."""
+ """
+ Make the FITS header for the HPX image.
+ """
header = header.copy(strip=True)
nside = header["NSIDE"]
# set pixel transformation parameters
@@ -312,7 +320,7 @@ def _make_hpx_header(header, append_history=None, append_comment=None):
]
for comment in comments:
header.add_comment(comment)
- #
+
if append_history is not None:
logger.info("HPX FITS header: append history")
for history in append_history:
@@ -326,7 +334,9 @@ def _make_hpx_header(header, append_history=None, append_comment=None):
def _make_healpix_header(header, nside,
append_history=None, append_comment=None):
- """Make the FITS header for the HEALPix data."""
+ """
+ Make the FITS header for the HEALPix data.
+ """
header = header.copy(strip=True)
# set HEALPix parameters
header["PIXTYPE"] = ("HEALPIX", "HEALPix pixelization")
@@ -353,7 +363,8 @@ def _make_healpix_header(header, nside,
@nb.jit(nb.int64(nb.int64), nopython=True)
def nside2npix(nside):
- """Calculate the number of pixels for the given Nside resolution.
+ """
+ Calculate the number of pixels for the given Nside resolution.
NOTE
----
@@ -364,7 +375,8 @@ def nside2npix(nside):
@nb.jit(nb.int64(nb.int64, nb.float64, nb.float64), nopython=True)
def ang2pix_ring_single(nside, theta, phi):
- """Calculate the pixel indexes in RING ordering scheme for one single
+ """
+ Calculate the pixel indexes in RING ordering scheme for one single
pair of angular coordinate on the sphere.
Parameters
@@ -415,18 +427,18 @@ def ang2pix_ring_single(nside, theta, phi):
iring = jp + jm + 1
ip = int(tt * iring)
ip = np.remainder(ip, 4*iring)
- #
if z > 0:
ipix = 2 * iring * (iring-1) + ip
else:
ipix = 12 * nside * nside - 2 * iring * (iring+1) + ip
- #
+
return ipix
@nb.jit(nb.types.UniTuple(nb.float64, 2)(nb.int64, nb.int64), nopython=True)
def pix2ang_ring_single(nside, ipix):
- """Calculate the angular coordinate on the sphere for one pixel index
+ """
+ Calculate the angular coordinate on the sphere for one pixel index
in the RING ordering scheme.
Parameters
@@ -485,7 +497,7 @@ def pix2ang_ring_single(nside, ipix):
iphi = 4 * iring + 1 - (ip - 2 * iring * (iring-1))
z = iring * iring * fact2 - 1.0
phi = (iphi - 0.5) * np.pi / (2 * iring)
- #
+
theta = np.arccos(z)
return (theta, phi)
@@ -494,7 +506,8 @@ def pix2ang_ring_single(nside, ipix):
nb.int64[:, :](nb.int64, nb.float64[:, :], nb.float64[:, :])],
nopython=True)
def ang2pix_ring(nside, theta, phi):
- """Calculate the pixel indexes in RING ordering scheme for each
+ """
+ Calculate the pixel indexes in RING ordering scheme for each
pair of angular coordinates on the sphere.
Parameters
@@ -530,7 +543,8 @@ def ang2pix_ring(nside, theta, phi):
nb.types.UniTuple(nb.float64[:, :], 2)(nb.int64, nb.int64[:, :])],
nopython=True)
def pix2ang_ring(nside, ipix):
- """Calculate the angular coordinates on the sphere for each pixel
+ """
+ Calculate the angular coordinates on the sphere for each pixel
index in the RING ordering scheme.
Parameters