#!/usr/bin/env python3
"""
plot_cmb_map.py
Downloads the Planck 2018 SMICA CMB temperature map (nosz version) from IRSA,
reads it with Astropy + Healpy, projects it, and produces a 4096×4096 PNG
using every 24-bit RGB color exactly once—mapping pixel intensities
by sorting the full RGB spectrum by luminance to create a meaningful heat
scale and assigning colors by rank order of CMB temperature.
The projected map (~2048×4096) is resampled (each row repeated) to
exactly fill a square 4096×4096 canvas before color mapping.
"""
import os
import requests
from astropy.io import fits
import healpy as hp
import numpy as np
from PIL import Image
# IRSA URL for the SMICA “nosz” full-mission CMB map
URL = (
"https://irsa.ipac.caltech.edu/data/Planck/release_3/"
"all-sky-maps/maps/component-maps/cmb/"
"COM_CMB_IQU-smica-nosz_2048_R3.00_full.fits"
)
FILENAME = os.path.basename(URL)
OUTFILE = "cmb_rgb_4096.png"
def download_map(url: str, filename: str) -> None:
"""Download the FITS file if not present."""
if not os.path.exists(filename):
r = requests.get(url, stream=True)
r.raise_for_status()
with open(filename, "wb") as f:
for chunk in r.iter_content(4096):
f.write(chunk)
def read_cmb_map(filename: str):
"""Read I_STOKES from the FITS file."""
with fits.open(filename, memmap=True) as hdul:
return hdul[1].data["I_STOKES"]
def project_map(cmb_data, xsize=4096):
"""Project CMB data to a plate carrée array (~2048×xsize)."""
return hp.cartview(
cmb_data,
return_projected_map=True,
xsize=xsize,
cbar=False,
title=None,
notext=True
)
def resample_to_square(proj):
"""Resample projected map by repeating each row to reach square shape."""
h, w = proj.shape
# Compute repeat factor (should be integer, here 2)
factor = w // h
return np.repeat(proj, factor, axis=0)
def map_to_rgb(proj_square):
"""
Map resampled data to every 24-bit RGB color exactly once.
Flatten the data, rank-order by temperature,
generate all RGB triplets sorted by luminance,
and assign each pixel a unique color.
"""
h, w = proj_square.shape
N = h * w
flat = proj_square.ravel()
order = np.argsort(flat)
idx = np.arange(N, dtype=np.uint32)
r = (idx >> 16) & 0xFF
g = (idx >> 8) & 0xFF
b = idx & 0xFF
lum = 0.2126 * r + 0.7152 * g + 0.0722 * b
color_order = np.argsort(lum)
rgb_flat = np.empty((N, 3), dtype=np.uint8)
colors = np.stack((r, g, b), axis=1)[color_order]
rgb_flat[order] = colors
return rgb_flat.reshape((h, w, 3))
def save_image(rgb_img, output: str = OUTFILE):
"""Save the RGB image at exact resolution."""
img = Image.fromarray(rgb_img, mode='RGB')
img.save(output)
def main():
download_map(URL, FILENAME)
cmb = read_cmb_map(FILENAME)
proj = project_map(cmb, xsize=4096)
proj_sq = resample_to_square(proj)
rgb = map_to_rgb(proj_sq)
save_image(rgb)
print(f"Saved RGB image: {OUTFILE} (4096×4096, full 24-bit palette)")
if __name__ == "__main__":
main()
