如何基于自定义分箱生成天体物理研究用2D直方图？

针对你提到的自定义等面积分箱无法统计数据生成直方图的问题，核心思路是先将天体坐标匹配到自定义分箱，统计每个分箱的天体数量，再基于分箱形状填充对应颜色，以下是具体实现步骤和优化方案：

1. 分箱归属判断（核心步骤） #

根据你分箱的定义方式（平面投影坐标/球面坐标），选择对应的点归属判断方法：

方法A：平面投影坐标下的批量匹配 #

如果你的自定义分箱是基于Aitoff等投影的平面多边形，可将天体坐标转换为同一投影的xy值，用矢量方法批量判断点是否在分箱内：

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
from shapely.vectorized import contains
from matplotlib.colors import Normalize

# 1. 准备数据（示例）
# 假设你已生成自定义分箱的平面边界：每个分箱是一组(x,y)坐标对，转为Shapely Polygon
bins_polygons = [Polygon(boundary_xy) for boundary_xy in your_custom_bin_boundaries]
# 天体的赤经、赤纬数组（单位：度）
ra, dec = your_object_ra, your_object_dec

# 2. 将天体坐标转换为Aitoff投影的xy值
def ra_dec_to_aitoff(ra_deg, dec_deg):
    ra_rad = np.radians(ra_deg)
    dec_rad = np.radians(dec_deg)
    # 利用matplotlib的Aitoff投影转换逻辑
    ax = plt.gca(projection='aitoff')
    x, y = ax(ra_rad, dec_rad)
    return x, y

x, y = ra_dec_to_aitoff(ra, dec)

# 3. 批量统计每个分箱的天体数量
bin_counts = np.zeros(len(bins_polygons), dtype=int)
for idx, poly in enumerate(bins_polygons):
    # 批量判断所有点是否在当前分箱内
    in_bin_mask = contains(poly, x, y)
    bin_counts[idx] = np.sum(in_bin_mask)

方法B：球面坐标下的精确匹配 #

如果分箱是基于球面定义的（避免投影误差），用Astropy的天区工具直接在球面坐标系判断归属：

import numpy as np
from astropy.coordinates import SkyCoord, PolygonSkyRegion
import astropy.units as u

# 1. 定义球面分箱
# 每个分箱的顶点赤经、赤纬列表（单位：度）
bin_vertices_list = your_spherical_bin_vertices  # 格式：[[ra1, ra2,...], [dec1, dec2,...]]
bin_regions = []
for ra_verts, dec_verts in bin_vertices_list:
    vertices = SkyCoord(ra=ra_verts*u.deg, dec=dec_verts*u.deg)
    bin_regions.append(PolygonSkyRegion(vertices=vertices))

# 2. 天体坐标转为SkyCoord
obj_coords = SkyCoord(ra=your_object_ra*u.deg, dec=your_object_dec*u.deg)

# 3. 统计分箱计数
bin_counts = np.zeros(len(bin_regions), dtype=int)
for idx, region in enumerate(bin_regions):
    in_bin_mask = region.contains(obj_coords)
    bin_counts[idx] = np.sum(in_bin_mask)

2. 绘制带颜色映射的2D直方图 #

基于统计得到的bin_counts，用你熟悉的plt.fill或plt.fill_between绘制分箱，并通过颜色映射展示计数：

# 初始化颜色归一化和色图
norm = Normalize(vmin=bin_counts.min(), vmax=bin_counts.max())
cmap = plt.get_cmap('viridis')

# 创建Aitoff投影画布
fig, ax = plt.subplots(subplot_kw={'projection': 'aitoff'}, figsize=(10, 6))

# 遍历分箱绘制并上色
for idx, poly in enumerate(bins_polygons):
    poly_x, poly_y = poly.exterior.xy
    # 填充分箱，颜色对应计数，添加白色边界区分
    ax.fill(poly_x, poly_y, color=cmap(norm(bin_counts[idx])), edgecolor='white', linewidth=0.3)

# 添加颜色条
sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
sm.set_array([])
fig.colorbar(sm, ax=ax, label='天体数量', pad=0.05)

# 调整样式
ax.grid(True, linestyle='--', alpha=0.6)
ax.set_title('自定义等面积分箱河外天体2D直方图', y=1.05)
plt.show()

3. 优化建议 #

• 若天体数量超过10万级，优先用shapely.vectorized.contains或Astropy的批量匹配，避免单循环遍历，提升效率；
• 分箱边界建议提前保存为数组或Shapely对象，重复使用时无需重新生成；
• 若需要对数刻度的颜色映射，将Normalize替换为LogNorm即可。

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