"""
get_interp_val — interpolation bilinéaire sur une grille HEALPix (NESTED)
via healpix-geo pour les conversions de coordonnées.
Équivalent à healpy.get_interp_val / get_interp_weights.
Deux régimes, selon `ellipsoid` :
- ellipsoid="sphere" (défaut) : reproduit EXACTEMENT l'algorithme RING de
référence de healpy (`T_Healpix_Base::get_interpol`, healpix_base.cc),
transcrit ici en pur NumPy (aucune dépendance à `healpy` au runtime).
Les 4 pixels renvoyés sont rigoureusement identiques à
`healpy.get_interp_weights(nside, lon, lat, lonlat=True, nest=True)`
(schéma NESTED, obtenu par conversion RING->NESTED via
`healpix_geo.ring.to_nested`) ; les poids concordent à ~1e-12 près
(bruit de calcul flottant inévitable tant qu'on ne relie pas le même
binaire C++ — largement sous la précision utile de n'importe quelle
application scientifique).
- ellipsoid != "sphere" (ex. "WGS84") : healpy n'a aucune notion
d'ellipsoïde, donc "identique à healpy" n'a pas de sens dans ce cas.
On délègue alors directement à `healpix_geo.nested.bilinear_interpolation`,
qui implémente une interpolation bilinéaire géométriquement correcte
(vérifiée indépendamment, voir la discussion qui a mené à ce fichier) et
supporte nativement les ellipsoïdes de référence.
Dépendances :
pip install healpix-geo # inclut cdshealpix comme dépendance
pip install numpy
Auteurs : Claude (Anthropic)
"""
import numpy as np
from healpix_geo import ring as _hg_ring
from healpix_geo.nested import bilinear_interpolation as _nested_bilinear_interpolation
# ---------------------------------------------------------------------------
# Portage NumPy fidèle de T_Healpix_Base<I>::get_interpol (schéma RING)
#
# Transcrit ligne à ligne depuis la source de référence :
# https://github.com/healpy/healpixmirror/blob/main/src/cxx/Healpix_cxx/healpix_base.cc
# (fonctions ring_above, get_ring_info2, get_interpol). Ne PAS modifier ces
# formules sans revalider contre healpy (voir tests associés) : plusieurs
# constantes (fact1_, fact2_, le "shift" par demi-pixel) sont non triviales
# et une erreur y est silencieuse (résultat plausible mais faux).
# ---------------------------------------------------------------------------
def _ring_above(z: np.ndarray, nside: int) -> np.ndarray:
"""Numéro de l'anneau immédiatement au nord (ou au pôle) de z=cos(theta)."""
az = np.abs(z)
twothird = 2.0 / 3.0
out = np.empty(z.shape, dtype=np.int64)
eq = az <= twothird
# I(nside*(2-1.5*z)) ; toujours positif => trunc == floor
out[eq] = np.trunc(nside * (2.0 - 1.5 * z[eq])).astype(np.int64)
pol = ~eq
iring = np.trunc(nside * np.sqrt(3.0 * (1.0 - az[pol]))).astype(np.int64)
out[pol] = np.where(z[pol] > 0, iring, 4 * nside - iring - 1)
return out
def _get_ring_info2(ring: np.ndarray, nside: int):
"""
Renvoie, pour chaque numéro d'anneau RING (1 <= ring <= 4*nside-1) :
startpix, ringpix (nb de pixels sur l'anneau), theta (colatitude du
centre de l'anneau), shifted (les pixels sont-ils décalés d'un
demi-pas de phi par rapport à phi=0 ?).
"""
ncap = 2 * nside * (nside - 1)
npix = 12 * nside * nside
fact1 = 2.0 / (3.0 * nside)
fact2 = 1.0 / (3.0 * nside * nside)
northring = np.where(ring > 2 * nside, 4 * nside - ring, ring)
theta = np.empty(ring.shape, dtype=np.float64)
ringpix = np.empty(ring.shape, dtype=np.int64)
startpix = np.empty(ring.shape, dtype=np.int64)
shifted = np.empty(ring.shape, dtype=bool)
cap = northring < nside
nr = northring[cap].astype(np.float64)
tmp = nr * nr * fact2
costheta = 1.0 - tmp
sintheta = np.sqrt(tmp * (2.0 - tmp))
theta[cap] = np.arctan2(sintheta, costheta)
ringpix[cap] = (4 * northring[cap]).astype(np.int64)
shifted[cap] = True
startpix[cap] = (2 * northring[cap] * (northring[cap] - 1)).astype(np.int64)
eqz = ~cap
theta[eqz] = np.arccos((2 * nside - northring[eqz]) * fact1)
ringpix[eqz] = 4 * nside
shifted[eqz] = ((northring[eqz] - nside) & 1) == 0
startpix[eqz] = ncap + (northring[eqz] - nside) * ringpix[eqz]
south = northring != ring
theta[south] = np.pi - theta[south]
startpix[south] = npix - startpix[south] - ringpix[south]
return startpix, ringpix, theta, shifted
def _trunc_cast(tmp: np.ndarray) -> np.ndarray:
"""Equivalent de I(tmp) en C++ pour un double : troncature vers zéro."""
return np.trunc(tmp).astype(np.int64)
def _get_interpol_ring(theta: np.ndarray, phi: np.ndarray, nside: int):
"""
Portage vectorisé de T_Healpix_Base<I>::get_interpol (schéma RING).
Parameters
----------
theta : colatitude en radians, shape (N,), dans [0, pi].
phi : longitude en radians, shape (N,), déjà repliée dans [0, 2*pi).
nside : int
Returns
-------
pix : np.ndarray uint64, shape (N, 4) — pixels en schéma RING.
wgt : np.ndarray float64, shape (N, 4).
"""
N = theta.shape[0]
npix = 12 * nside * nside
twopi = 2.0 * np.pi
z = np.cos(theta)
ir1 = _ring_above(z, nside)
ir2 = ir1 + 1
pix = np.zeros((N, 4), dtype=np.int64)
wgt = np.zeros((N, 4), dtype=np.float64)
theta1 = np.zeros(N, dtype=np.float64)
theta2 = np.zeros(N, dtype=np.float64)
m1 = ir1 > 0
if m1.any():
sp, nr, th1, shift = _get_ring_info2(ir1[m1], nside)
theta1[m1] = th1
dphi = twopi / nr
half_shift = np.where(shift, 0.5, 0.0)
tmp = phi[m1] / dphi - half_shift
i1 = np.where(tmp < 0, _trunc_cast(tmp) - 1, _trunc_cast(tmp))
w1 = (phi[m1] - (i1 + half_shift) * dphi) / dphi
i2 = i1 + 1
i1 = np.where(i1 < 0, i1 + nr, i1)
i2 = np.where(i2 >= nr, i2 - nr, i2)
pix[m1, 0] = sp + i1
pix[m1, 1] = sp + i2
wgt[m1, 0] = 1 - w1
wgt[m1, 1] = w1
m2 = ir2 < 4 * nside
if m2.any():
sp, nr, th2, shift = _get_ring_info2(ir2[m2], nside)
theta2[m2] = th2
dphi = twopi / nr
half_shift = np.where(shift, 0.5, 0.0)
tmp = phi[m2] / dphi - half_shift
i1 = np.where(tmp < 0, _trunc_cast(tmp) - 1, _trunc_cast(tmp))
w1 = (phi[m2] - (i1 + half_shift) * dphi) / dphi
i2 = i1 + 1
i1 = np.where(i1 < 0, i1 + nr, i1)
i2 = np.where(i2 >= nr, i2 - nr, i2)
pix[m2, 2] = sp + i1
pix[m2, 3] = sp + i2
wgt[m2, 2] = 1 - w1
wgt[m2, 3] = w1
north_pole = ir1 == 0
if north_pole.any():
wtheta = theta[north_pole] / theta2[north_pole]
wgt[north_pole, 2] *= wtheta
wgt[north_pole, 3] *= wtheta
fac = (1 - wtheta) * 0.25
wgt[north_pole, 0] = fac
wgt[north_pole, 1] = fac
wgt[north_pole, 2] += fac
wgt[north_pole, 3] += fac
pix[north_pole, 0] = (pix[north_pole, 2] + 2) & 3
pix[north_pole, 1] = (pix[north_pole, 3] + 2) & 3
south_pole = (ir2 == 4 * nside) & ~north_pole
if south_pole.any():
wtheta = (theta[south_pole] - theta1[south_pole]) / (np.pi - theta1[south_pole])
wgt[south_pole, 0] *= (1 - wtheta)
wgt[south_pole, 1] *= (1 - wtheta)
fac = wtheta * 0.25
wgt[south_pole, 0] += fac
wgt[south_pole, 1] += fac
wgt[south_pole, 2] = fac
wgt[south_pole, 3] = fac
pix[south_pole, 2] = ((pix[south_pole, 0] + 2) & 3) + npix - 4
pix[south_pole, 3] = ((pix[south_pole, 1] + 2) & 3) + npix - 4
mid = ~north_pole & ~south_pole
if mid.any():
wtheta = (theta[mid] - theta1[mid]) / (theta2[mid] - theta1[mid])
wgt[mid, 0] *= (1 - wtheta)
wgt[mid, 1] *= (1 - wtheta)
wgt[mid, 2] *= wtheta
wgt[mid, 3] *= wtheta
return pix.astype(np.uint64), wgt
# ---------------------------------------------------------------------------
# API principale
# ---------------------------------------------------------------------------
[docs]
def get_interp_weights(
lon: np.ndarray,
lat: np.ndarray,
depth: int,
ellipsoid: str = "sphere",
) -> tuple[np.ndarray, np.ndarray]:
"""
Renvoie les 4 cellules HEALPix et leurs poids d'interpolation bilinéaire
pour chaque position (lon, lat).
Pour ellipsoid="sphere" (défaut), les résultats (pixels ET poids) sont
rigoureusement identiques à
``healpy.get_interp_weights(nside, lon, lat, lonlat=True, nest=True)``
(aux erreurs d'arrondi flottant près, ~1e-12) : c'est un portage direct
de l'algorithme RING de référence de healpy, pas une réinvention
géométrique approximative.
Pour un ellipsoïde non-sphérique, il n'existe pas de "référence healpy"
(healpy ne connaît que la sphère) : on utilise alors
``healpix_geo.nested.bilinear_interpolation``, qui implémente une
interpolation bilinéaire correcte sur la grille NESTED en tenant compte
de l'ellipsoïde choisi.
Parameters
----------
lon : np.ndarray, shape (N,)
Longitudes en degrés.
lat : np.ndarray, shape (N,)
Latitudes en degrés.
depth : int
Profondeur HEALPix (nside = 2**depth).
ellipsoid : str, optional
Ellipsoïde de référence : "sphere" (défaut, identique à healpy) ou
"WGS84", "GRS80", etc. Voir la doc de healpix-geo pour la liste.
Returns
-------
pixels : np.ndarray of uint64, shape (N, 4)
Indices des 4 cellules HEALPix (schéma NESTED).
weights : np.ndarray of float64, shape (N, 4)
Poids bilinéaires correspondants. Chaque ligne somme à 1.
Raises
------
ValueError
Si lon et lat n'ont pas la même forme.
"""
lon = np.asarray(lon, dtype=np.float64)
lat = np.asarray(lat, dtype=np.float64)
if lon.shape != lat.shape:
raise ValueError(
f"lon et lat doivent avoir la même forme "
f"(got {lon.shape} vs {lat.shape})"
)
lon_flat = lon.ravel()
lat_flat = lat.ravel()
if ellipsoid == "sphere":
nside = 2 ** depth
theta = np.radians(90.0 - lat_flat)
phi = np.radians(lon_flat) % (2.0 * np.pi)
pix_ring, weights = _get_interpol_ring(theta, phi, nside)
pixels = _hg_ring.to_nested(pix_ring.ravel(), depth).reshape(pix_ring.shape)
return pixels.astype(np.uint64), weights
# Ellipsoïdes non-sphériques : healpy n'a pas d'équivalent, on délègue
# à l'implémentation native de healpix-geo.
cell_ids, weights = _nested_bilinear_interpolation(
lon_flat, lat_flat, depth, ellipsoid=ellipsoid
)
pixels = np.asarray(cell_ids.data if hasattr(cell_ids, "data") else cell_ids, dtype=np.uint64)
weights = np.asarray(weights.data if hasattr(weights, "data") else weights, dtype=np.float64)
return pixels, weights
[docs]
def get_interp_val(
hpx_map: np.ndarray,
lon,
lat,
depth: int,
ellipsoid: str = "sphere",
) -> np.ndarray:
"""
Interpolation bilinéaire d'une carte HEALPix aux coordonnées géographiques.
Équivalent à healpy.get_interp_val(m, theta, phi, nest=True, lonlat=True)
mais utilise healpix-geo pour les projections, ce qui permet de travailler
sur des ellipsoïdes de référence (ex. WGS84).
Parameters
----------
hpx_map : np.ndarray, shape (12 * 4**depth,)
Carte HEALPix en ordre NESTED.
lon : float ou np.ndarray
Longitude(s) en degrés.
lat : float ou np.ndarray
Latitude(s) en degrés.
depth : int
Profondeur HEALPix (nside = 2**depth).
ellipsoid : str, optional
Ellipsoïde de référence : "sphere" (défaut) ou "WGS84", etc.
Returns
-------
np.ndarray ou float
Valeurs interpolées. Scalaire si lon/lat sont scalaires, sinon array de
même forme que lon/lat.
Examples
--------
>>> import numpy as np
>>> from healpix_interp import get_interp_val
>>>
>>> depth = 3
>>> nside = 2**depth
>>> npix = 12 * nside**2
>>> hpx_map = np.arange(npix, dtype=float)
>>>
>>> # Point unique
>>> val = get_interp_val(hpx_map, lon=45.0, lat=30.0, depth=depth)
>>> print(val)
>>>
>>> # Grille de points
>>> lons = np.linspace(0, 360, 50)
>>> lats = np.linspace(-80, 80, 40)
>>> lon_grid, lat_grid = np.meshgrid(lons, lats)
>>> vals = get_interp_val(hpx_map, lon_grid, lat_grid, depth=depth)
>>> print(vals.shape) # (40, 50)
>>>
>>> # Avec ellipsoïde WGS84
>>> vals_wgs84 = get_interp_val(
... hpx_map, lon=2.3522, lat=48.8566, depth=depth, ellipsoid="WGS84"
... )
Notes
-----
La carte hpx_map doit contenir exactement 12 * 4**depth éléments et être
en ordre NESTED (le schéma utilisé par healpix-geo).
Si votre carte est en ordre RING (healpy), convertissez-la d'abord :
import healpy as hp
hpx_map_nested = hp.reorder(hpx_map_ring, r2n=True)
"""
lon = np.asarray(lon, dtype=np.float64)
lat = np.asarray(lat, dtype=np.float64)
scalar_input = lon.ndim == 0
original_shape = lon.shape
lon = np.atleast_1d(lon).ravel()
lat = np.atleast_1d(lat).ravel()
pixels, weights = get_interp_weights(lon, lat, depth, ellipsoid=ellipsoid)
# hpx_map[pixels] : (N, 4) — valeurs aux 4 cellules
vals = np.sum(weights * hpx_map[pixels], axis=1)
if scalar_input:
return float(vals[0])
return vals.reshape(original_shape)
# ---------------------------------------------------------------------------
# Test rapide (optionnel)
# ---------------------------------------------------------------------------
if __name__ == "__main__":
import time
print("=== Test get_interp_val avec healpix-geo ===\n")
depth = 5
nside = 2**depth
npix = 12 * nside**2
rng = np.random.default_rng(42)
hpx_map = rng.standard_normal(npix)
# Grille de test
lons = np.linspace(0.1, 359.9, 200)
lats = np.linspace(-89.0, 89.0, 200)
lon_grid, lat_grid = np.meshgrid(lons, lats)
# --- healpix-geo (sphere) ---
t0 = time.perf_counter()
vals_sphere = get_interp_val(hpx_map, lon_grid, lat_grid, depth=depth, ellipsoid="sphere")
t1 = time.perf_counter()
print(f"healpix-geo (sphere) : shape={vals_sphere.shape} "
f"min={vals_sphere.min():.4f} max={vals_sphere.max():.4f} "
f"[{(t1-t0)*1000:.1f} ms]")
# --- healpix-geo (WGS84) ---
t0 = time.perf_counter()
vals_wgs84 = get_interp_val(hpx_map, lon_grid, lat_grid, depth=depth, ellipsoid="WGS84")
t1 = time.perf_counter()
print(f"healpix-geo (WGS84) : shape={vals_wgs84.shape} "
f"min={vals_wgs84.min():.4f} max={vals_wgs84.max():.4f} "
f"[{(t1-t0)*1000:.1f} ms]")
# Différence sphere vs WGS84
diff = np.abs(vals_sphere - vals_wgs84)
print(f"\nDiff sphere vs WGS84 : mean={diff.mean():.6f} max={diff.max():.6f}")
# --- Comparaison avec healpy (si disponible) ---
try:
import healpy as hp
# healpy utilise le schéma RING par défaut ; on convertit la carte
hpx_map_ring = hp.reorder(hpx_map, n2r=True)
lons_flat = lon_grid.ravel()
lats_flat = lat_grid.ravel()
# healpy attend colatitude en radians et longitude en radians
theta = np.radians(90.0 - lats_flat)
phi = np.radians(lons_flat)
vals_healpy = hp.get_interp_val(hpx_map_ring, theta, phi).reshape(lon_grid.shape)
diff_hp = np.abs(vals_sphere - vals_healpy)
print(f"\nDiff vs healpy (sphere): mean={diff_hp.mean():.6f} max={diff_hp.max():.6f}")
except ImportError:
print("\n(healpy non disponible pour comparaison)")
print("\n=== OK ===")