如何在2D numpy数组创建25x25像素网格并统计各单元格数值

这里提供两种主流实现路径，按需选择即可，不需要引入geopandas这类非天文领域的重依赖库。

方案1：纯NumPy实现（最高效，无额外依赖） #

直接对NumPy数组做分块处理，是性能最优的方案，适合只需要统计结果不需要生成天文兼容区域文件的场景：

import numpy as np
import matplotlib.pyplot as plt

# 替换为你自己读取到的FITS图像2D数组
fits_data = np.random.randint(0, 100, size=(500, 500)) 
grid_size = 25

rows, cols = fits_data.shape
# 计算可被25整除的有效尺寸，不足的部分可自行选择截断或补0
valid_rows = (rows // grid_size) * grid_size
valid_cols = (cols // grid_size) * grid_size
cropped_data = fits_data[:valid_rows, :valid_cols]

# 把原数组切分为(网格行数, 网格列数, 25, 25)的分块数组
grid_blocks = cropped_data.reshape(valid_rows//grid_size, grid_size, valid_cols//grid_size, grid_size).swapaxes(1,2)

# 统计每个网格的数值，可按需替换为mean()/max()/std()/np.median()等统计逻辑
grid_stats = grid_blocks.sum(axis=(2,3))

# 直接生成分布映射图
plt.imshow(grid_stats, origin='lower')
plt.colorbar(label='网格内计数统计值')
plt.title('25x25像素网格统计分布')
plt.show()

方案2：pyregion批量生成网格（适合需要保留DS9兼容区域文件的场景） #

如果你需要导出网格区域给DS9等天文软件交互使用，可以批量生成所有网格区域一次性统计，不需要循环单个区域：

import numpy as np
import pyregion
from astropy.io import fits
import matplotlib.pyplot as plt

fits_data = np.random.randint(0, 100, size=(500, 500)) 
grid_size = 25
rows, cols = fits_data.shape

# 批量生成DS9格式的网格区域字符串
region_content = "image\n"
for y_center in range(grid_size//2, rows, grid_size):
    for x_center in range(grid_size//2, cols, grid_size):
        # DS9正方形区域格式：box(中心x, 中心y, 宽度, 高度, 旋转角)，DS9像素坐标从1开始计数
        region_content += f"box({x_center+1}, {y_center+1}, {grid_size}, {grid_size}, 0)\n"

# 读取批量生成的区域
regs = pyregion.parse(region_content).as_imagecoord(fits_data.shape)
# 一次性生成所有网格的掩膜数组，形状为(网格总数, 图像行数, 图像列数)
mask = regs.get_mask(hdu=fits.PrimaryHDU(fits_data))

# 统计每个网格的数值
grid_stats = np.array([fits_data[mask[i]].sum() for i in range(mask.shape[0])])
# 重塑为二维网格形状用于绘图
grid_stats = grid_stats.reshape(rows//grid_size, cols//grid_size)

# 绘图逻辑同方案1
plt.imshow(grid_stats, origin='lower')
plt.colorbar(label='网格内计数统计值')
plt.title('25x25像素网格统计分布')
plt.show()

不推荐使用geopandas的说明 #

geopandas是面向地理空间矢量数据的处理库，用于天文像素网格统计属于过度适配，需要额外做像素坐标转矢量多边形、空间匹配等冗余步骤，性能远低于前两种方案，没有必要引入。

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