如何使用HealPy将FITS格式的掩膜/地图分割为两部分？

Healpix格式的FITS掩膜确实不会直接存储经纬度数据，但每个像素的位置是由nside参数和像素索引唯一确定的，咱们可以用HealPy自带的坐标转换工具，把像素索引转成黄道坐标系的经纬度，再按区域筛选分割。具体步骤如下：

1. 读取原掩膜并获取基础参数 #

先读取你的掩膜文件，拿到Healpix网格的分辨率参数nside和所有像素的索引：

import healpy as hp
import numpy as np

# 读取原掩膜
wmap_map_I = hp.read_map('mask.fits')
nside = hp.get_nside(wmap_map_I)
# 生成所有像素的索引
all_pix = np.arange(hp.nside2npix(nside))

2. 将像素索引转换为黄道坐标系的经纬度 #

原绘图用了coord=["E"]（黄道坐标系），所以咱们需要把默认的赤道坐标转成黄道坐标：

# 获取赤道坐标系下的极角(theta)和方位角(phi)
theta_equ, phi_equ = hp.pix2ang(nside, all_pix)
# 创建坐标转换器：赤道(C)转黄道(E)
rot = hp.Rotator(coord=['C', 'E'])
# 转换为黄道坐标系的极角和方位角
theta_ecl, phi_ecl = rot(theta_equ, phi_equ)
# 转换为更直观的黄道纬度(范围-90~90°)和经度(范围0~360°)
lat_ecl = 90 - np.degrees(theta_ecl)
lon_ecl = np.degrees(phi_ecl)

3. 根据区域范围筛选像素 #

对照你的掩膜图，定义A、B区域的黄道经纬度条件（以下是示例条件，你可以根据图中实际区域调整数值）：

# 示例：A区域为黄道经度0°~180°、纬度-30°~30°
condition_A = (lon_ecl >= 0) & (lon_ecl <= 180) & (lat_ecl >= -30) & (lat_ecl <= 30)
# 示例：B区域为黄道经度180°~360°、纬度-30°~30°
condition_B = (lon_ecl >= 180) & (lon_ecl <= 360) & (lat_ecl >= -30) & (lat_ecl <= 30)

4. 生成并保存新掩膜 #

复制原掩膜，把不符合区域条件的像素设为无效值（原掩膜用-1表示无效，这里保持一致）：

# 生成A区域掩膜
mask_A = wmap_map_I.copy()
mask_A[~condition_A] = -1
# 生成B区域掩膜
mask_B = wmap_map_I.copy()
mask_B[~condition_B] = -1

# 保存为新的FITS文件
hp.write_map('mask_A.fits', mask_A, overwrite=True)
hp.write_map('mask_B.fits', mask_B, overwrite=True)

验证分割结果 #

你可以用原代码的projview函数分别绘制新掩膜，确认是否只保留了目标区域：

from healpy.newvisufunc import projview

projview(mask_A,
         graticule=True, graticule_labels=True, xlabel="longitude", ylabel="latitude",
         projection_type="mollweide", coord=["E"],
         title="Mask A (Ecliptic)", norm="hist", min=-1, max=1)

内容的提问来源于stack exchange，提问作者NeStack