Query Disc Functions

Query Disc Functions#

This module provides functions for querying pixels within circular regions on the sphere, essential for astronomical applications involving point sources, aperture photometry, and local analysis.

Query disc implementation for HEALPix spherical disc queries.

jax_healpy._query_disc.query_disc(nside, vec, radius, inclusive=False, fact=4, nest=False, max_length=None)[source]#

Find pixels within a disc on the sphere.

This function supports both single and batched queries. It is fully JIT-compatible and differentiable with respect to vec and radius parameters.

Parameters:

  • nside (int) – The resolution parameter of the HEALPix map

  • vec (array-like) – Either a single three-component unit vector (3,) defining the center of the disc, or a batch of vectors (B, 3) defining B disc centers

  • radius (float) – The radius of the disc in radians

  • inclusive (bool, optional) – If False (default), return pixels whose centers lie within the disc. Results are guaranteed to match healpy exactly for inclusive=False. If True, return all pixels that overlap with the disc. Note: inclusive=True may produce slightly different results compared to healpy due to algorithm differences in determining pixel overlap.

  • fact (int, optional) – For inclusive queries, the pixelization factor (default: 4)

  • nest (bool, optional) – If True, assume NESTED pixel ordering, otherwise RING ordering (default: False)

  • max_length (int, optional) – Maximum number of pixels to return per disc. If None, defaults to npix. For batched inputs, this limits memory usage by returning (max_length, B) instead of (npix, B).

Returns:

ipix – For single vector input (3,): returns (max_length,) or (npix,) if max_length is None For batch vector input (B, 3): returns (max_length, B) Pixels outside the disc are marked as npix (sentinel value). This allows direct indexing like: map.at[disc].set(value)

Return type:

array

Raises:

NotImplementedError – If nest=True (nested ordering not yet supported)

Notes

This function currently only supports RING ordering. The function is JIT-compatible and differentiable when compiled with static_argnums=(0,) for the nside parameter. The returned array has fixed size for JIT compatibility - pixels outside the disc have value npix.

When indexing a JAX array, pixels outside the disc should be ignored. If indexing a numpy array, this will raise an out-of-bounds error.

Examples

Single disc query:

>>> import jax_healpy as hp
>>> import jax.numpy as jnp
>>> import jax
>>> nside = 16
>>> vec = jnp.array([1.0, 0.0, 0.0])  # Point on equator
>>> radius = 0.1  # ~5.7 degrees
>>> disc = hp.query_disc(nside, vec, radius)
>>> # Use directly for indexing
>>> map = jnp.zeros(hp.nside2npix(nside))
>>> map = map.at[disc].set(1.0)

Batch disc query:

>>> vecs = jnp.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])  # (2, 3) - two centers
>>> discs = hp.query_disc(nside, vecs, radius, max_length=1000)  # (1000, 2)
>>> # Each column contains pixels for one disc
>>> map1 = map.at[discs[:, 0]].set(1.0)  # First disc
>>> map2 = map.at[discs[:, 1]].set(2.0)  # Second disc

For JIT compilation:

>>> jit_query_disc = jax.jit(
...     lambda n, v, r: hp.query_disc(n, v, r)
... )
>>> disc_jit = jit_query_disc(nside, vec, radius)
jax_healpy._query_disc.estimate_disc_pixel_count(nside, radius)[source]#

Estimate number of pixels in a disc of given radius.

Uses the analytical approximation: n_approx = 6 * nside^2 * (1 - cos(radius))

This approximation is exact at radius = π (full sphere) and loses precision for smaller radii. It provides a good estimate for setting the max_length parameter in query_disc functions.

Parameters:

  • nside (int) – HEALPix nside parameter

  • radius (float) – Disc radius in radians

Returns:

pixel_count – Estimated number of pixels in the disc

Return type:

int

Notes

The approximation assumes uniform pixel density across the sphere, which is accurate for large discs but less precise for small discs where pixel shape variations matter more.

Examples

>>> import jax_healpy as hp
>>> nside = 64
>>> radius = 0.1  # ~5.7 degrees
>>> estimated_count = hp.estimate_disc_pixel_count(nside, radius)
>>> # Use as max_length for query_disc
>>> actual_disc = hp.query_disc(nside, [1, 0, 0], radius, max_length=estimated_count)
jax_healpy._query_disc.estimate_disc_radius(nside, pixel_count)[source]#

Estimate radius needed for a disc containing given pixel count.

Inverse of estimate_disc_pixel_count. Uses the analytical relationship: radius = arccos(1 - pixel_count / (6 * nside^2))

Parameters:

  • nside (int) – HEALPix nside parameter

  • pixel_count (int) – Desired number of pixels in the disc

Returns:

radius – Estimated disc radius in radians

Return type:

float

Notes

This is the inverse of estimate_disc_pixel_count and has the same accuracy characteristics: exact at full sphere, less precise for small discs.

Examples

>>> import jax_healpy as hp
>>> nside = 64
>>> target_pixels = 1000
>>> estimated_radius = hp.estimate_disc_radius(nside, target_pixels)
>>> # Verify the relationship
>>> back_pixels = hp.estimate_disc_pixel_count(nside, estimated_radius)
>>> abs(back_pixels - target_pixels) < 10  # Should be close
True

Disc Query Functions#

jax_healpy.query_disc(nside, vec, radius, inclusive=False, fact=4, nest=False, max_length=None)[source]#

Find pixels within a disc on the sphere.

This function supports both single and batched queries. It is fully JIT-compatible and differentiable with respect to vec and radius parameters.

Parameters:

  • nside (int) – The resolution parameter of the HEALPix map

  • vec (array-like) – Either a single three-component unit vector (3,) defining the center of the disc, or a batch of vectors (B, 3) defining B disc centers

  • radius (float) – The radius of the disc in radians

  • inclusive (bool, optional) – If False (default), return pixels whose centers lie within the disc. Results are guaranteed to match healpy exactly for inclusive=False. If True, return all pixels that overlap with the disc. Note: inclusive=True may produce slightly different results compared to healpy due to algorithm differences in determining pixel overlap.

  • fact (int, optional) – For inclusive queries, the pixelization factor (default: 4)

  • nest (bool, optional) – If True, assume NESTED pixel ordering, otherwise RING ordering (default: False)

  • max_length (int, optional) – Maximum number of pixels to return per disc. If None, defaults to npix. For batched inputs, this limits memory usage by returning (max_length, B) instead of (npix, B).

Returns:

ipix – For single vector input (3,): returns (max_length,) or (npix,) if max_length is None For batch vector input (B, 3): returns (max_length, B) Pixels outside the disc are marked as npix (sentinel value). This allows direct indexing like: map.at[disc].set(value)

Return type:

array

Raises:

NotImplementedError – If nest=True (nested ordering not yet supported)

Notes

This function currently only supports RING ordering. The function is JIT-compatible and differentiable when compiled with static_argnums=(0,) for the nside parameter. The returned array has fixed size for JIT compatibility - pixels outside the disc have value npix.

When indexing a JAX array, pixels outside the disc should be ignored. If indexing a numpy array, this will raise an out-of-bounds error.

Examples

Single disc query:

>>> import jax_healpy as hp
>>> import jax.numpy as jnp
>>> import jax
>>> nside = 16
>>> vec = jnp.array([1.0, 0.0, 0.0])  # Point on equator
>>> radius = 0.1  # ~5.7 degrees
>>> disc = hp.query_disc(nside, vec, radius)
>>> # Use directly for indexing
>>> map = jnp.zeros(hp.nside2npix(nside))
>>> map = map.at[disc].set(1.0)

Batch disc query:

>>> vecs = jnp.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])  # (2, 3) - two centers
>>> discs = hp.query_disc(nside, vecs, radius, max_length=1000)  # (1000, 2)
>>> # Each column contains pixels for one disc
>>> map1 = map.at[discs[:, 0]].set(1.0)  # First disc
>>> map2 = map.at[discs[:, 1]].set(2.0)  # Second disc

For JIT compilation:

>>> jit_query_disc = jax.jit(
...     lambda n, v, r: hp.query_disc(n, v, r)
... )
>>> disc_jit = jit_query_disc(nside, vec, radius)

Disc Estimation Functions#

jax_healpy.estimate_disc_pixel_count(nside, radius)[source]#

Estimate number of pixels in a disc of given radius.

Uses the analytical approximation: n_approx = 6 * nside^2 * (1 - cos(radius))

This approximation is exact at radius = π (full sphere) and loses precision for smaller radii. It provides a good estimate for setting the max_length parameter in query_disc functions.

Parameters:

  • nside (int) – HEALPix nside parameter

  • radius (float) – Disc radius in radians

Returns:

pixel_count – Estimated number of pixels in the disc

Return type:

int

Notes

The approximation assumes uniform pixel density across the sphere, which is accurate for large discs but less precise for small discs where pixel shape variations matter more.

Examples

>>> import jax_healpy as hp
>>> nside = 64
>>> radius = 0.1  # ~5.7 degrees
>>> estimated_count = hp.estimate_disc_pixel_count(nside, radius)
>>> # Use as max_length for query_disc
>>> actual_disc = hp.query_disc(nside, [1, 0, 0], radius, max_length=estimated_count)

Estimate the number of pixels within a circular disc on the sphere.

jax_healpy.estimate_disc_radius(nside, pixel_count)[source]#

Estimate radius needed for a disc containing given pixel count.

Inverse of estimate_disc_pixel_count. Uses the analytical relationship: radius = arccos(1 - pixel_count / (6 * nside^2))

Parameters:

  • nside (int) – HEALPix nside parameter

  • pixel_count (int) – Desired number of pixels in the disc

Returns:

radius – Estimated disc radius in radians

Return type:

float

Notes

This is the inverse of estimate_disc_pixel_count and has the same accuracy characteristics: exact at full sphere, less precise for small discs.

Examples

>>> import jax_healpy as hp
>>> nside = 64
>>> target_pixels = 1000
>>> estimated_radius = hp.estimate_disc_radius(nside, target_pixels)
>>> # Verify the relationship
>>> back_pixels = hp.estimate_disc_pixel_count(nside, estimated_radius)
>>> abs(back_pixels - target_pixels) < 10  # Should be close
True

Estimate the angular radius needed to contain a given number of pixels.

Find all pixels within a circular region on the sphere.

This function identifies all HEALPix pixels whose centers lie within a circular region defined by a center direction and angular radius.

Parameters#

nsideint

HEALPix resolution parameter

vecarray_like

Unit vector pointing to the disc center, shape (3,)

radiusfloat

Angular radius of the disc in radians

inclusivebool, optional

If True, include pixels that are partially within the disc. If False, only include pixels whose centers are within the disc. Default: False

factint, optional

Factor for subdividing pixels when inclusive=True. Higher values give more accurate results but are slower. Default: 4

nestbool, optional

If True, return pixel indices in NESTED ordering. If False, return in RING ordering. Default: False

Returns#

pixelsarray_like

Array of pixel indices within the specified disc

Examples#

Basic disc query:

import jax.numpy as jnp
import jax_healpy as hp

# Define disc center (North pole)
center_vec = jnp.array([0.0, 0.0, 1.0])

# Query pixels within 10 degree radius
radius = jnp.radians(10.0)  # Convert to radians
nside = 64

pixels = hp.query_disc(nside, center_vec, radius)
print(f"Found {len(pixels)} pixels in disc")

Using angular coordinates:

# Convert from theta, phi to unit vector
theta, phi = jnp.radians(45.0), jnp.radians(30.0)  # 45° from pole, 30° azimuth
center_vec = hp.ang2vec(theta, phi)

# Query with inclusive boundary
pixels = hp.query_disc(nside, center_vec, radius, inclusive=True)

Batch processing multiple discs:

# Define multiple disc centers
centers = jnp.array([
    [0.0, 0.0, 1.0],    # North pole
    [1.0, 0.0, 0.0],    # X-axis
    [0.0, 1.0, 0.0]     # Y-axis
])

# Query each disc
all_pixels = []
for center in centers:
    pixels = hp.query_disc(nside, center, radius)
    all_pixels.append(pixels)

Mathematical Background#

The query_disc function finds pixels whose centers satisfy:

\[\cos(\text{radius}) \leq \vec{v}_{\text{center}} \cdot \vec{v}_{\text{pixel}}\]

where \(\vec{v}_{\text{center}}\) is the disc center unit vector and \(\vec{v}_{\text{pixel}}\) is the pixel center unit vector.

For inclusive queries, the function also considers pixels that are partially within the disc boundary by subdividing pixel boundaries and checking if any subdivision points fall within the disc.

Applications#

Common use cases for disc queries include:

Point Source Analysis

# Find pixels around a known source position
source_ra, source_dec = 83.6331, 22.0145  # Crab Nebula (degrees)

# Convert to HEALPix coordinates
theta = jnp.radians(90.0 - source_dec)  # Convert Dec to colatitude
phi = jnp.radians(source_ra)
source_vec = hp.ang2vec(theta, phi)

# Query 1-degree radius around source
aperture_radius = jnp.radians(1.0)
source_pixels = hp.query_disc(nside, source_vec, aperture_radius)

Aperture Photometry

# Extract flux within aperture
sky_map = jnp.random.normal(0, 1, hp.nside2npix(nside))

aperture_pixels = hp.query_disc(nside, source_vec, aperture_radius)
aperture_flux = jnp.sum(sky_map[aperture_pixels])

print(f"Total flux in aperture: {aperture_flux}")

Local Statistics

# Compute local statistics around multiple positions
positions = jnp.array([
    hp.ang2vec(jnp.radians(30), jnp.radians(0)),
    hp.ang2vec(jnp.radians(60), jnp.radians(90)),
    hp.ang2vec(jnp.radians(90), jnp.radians(180))
])

local_stats = []
for pos in positions:
    pixels = hp.query_disc(nside, pos, aperture_radius)
    local_mean = jnp.mean(sky_map[pixels])
    local_std = jnp.std(sky_map[pixels])
    local_stats.append((local_mean, local_std))

Masking and Exclusion

# Create exclusion mask around bright sources
bright_sources = [
    hp.ang2vec(jnp.radians(45), jnp.radians(0)),
    hp.ang2vec(jnp.radians(135), jnp.radians(90))
]

exclusion_mask = jnp.ones(hp.nside2npix(nside))
exclusion_radius = jnp.radians(5.0)  # 5-degree exclusion

for source in bright_sources:
    excluded_pixels = hp.query_disc(nside, source, exclusion_radius)
    exclusion_mask = exclusion_mask.at[excluded_pixels].set(0)

Performance Considerations#

  • Query disc is optimized for GPU execution when JAX is configured with CUDA/ROCm

  • Performance scales with disc size and resolution (nside)

  • For many small discs, consider batching operations:

@jax.jit
def batch_query_disc(nside, centers, radius):
    # Vectorized implementation for multiple discs
    return jax.vmap(lambda c: hp.query_disc(nside, c, radius))(centers)

  • The inclusive parameter significantly affects computation time

  • For very large discs (> 60°), consider alternative approaches

Accuracy Notes#

  • Standard queries (inclusive=False) are exact for pixel centers

  • Inclusive queries provide better boundary accuracy at computational cost

  • Higher fact values in inclusive mode improve accuracy but slow computation

  • For critical applications, validate results against known geometric cases

See Also#

Contents