二维数组重采样求助：傅里叶变换法数据损耗异常

嘿，我来帮你搞定这个地理数据重采样的问题！首先得明确：你用FFT的方法完全走错方向了——fftpack.fft2的shape参数根本不是用来做空间重采样的，它是在频域对信号做零填充（放大时）或截断（缩小时），这和地理数据（比如人口）需要的重采样逻辑完全不搭，所以数据总和大幅下降是必然的结果。

针对地理科学里的这类数据（很多时候需要保持总和不变，比如人口、通量这类守恒型数据），我给你推荐两种靠谱的解决办法：

方法一：分块加权求和（精准保持数据总和） #

这种方法适用于任何比例的重采样，核心思路是计算新网格每个单元对应原始网格的覆盖范围，然后按重叠面积加权求和——如果新尺寸刚好是原始尺寸的整数倍，就会退化成你之前用的等比例方法，非常适配地理数据的需求。

示例代码（100×100 → 60×80） #

import numpy as np

# 复用你创建原始数据的代码
xc1, xc2, yc1, yc2 = 100, 110, 35, 45
XSIZE, YSIZE = 100, 100
lon, lat = np.linspace(xc1, xc2, XSIZE), np.linspace(yc1, yc2, YSIZE)
pop = np.random.uniform(low=1000, high=50000, size=(XSIZE*YSIZE,)).reshape(YSIZE, XSIZE)
original_sum = pop.sum()
print(f"原始数据总和: {original_sum}")

# 目标重采样尺寸
target_shape = (60, 80)  # (纬度方向, 经度方向)

# 计算原始网格和目标网格的像素边界（假设原始坐标是像素中心，所以边界要向外扩展半个步长）
lon_step = (xc2 - xc1) / (XSIZE - 1) if XSIZE > 1 else 0
lat_step = (yc2 - yc1) / (YSIZE - 1) if YSIZE > 1 else 0
lon_edges = np.linspace(xc1 - lon_step/2, xc2 + lon_step/2, XSIZE + 1)
lat_edges = np.linspace(yc1 - lat_step/2, yc2 + lat_step/2, YSIZE + 1)

# 计算目标网格的像素边界
target_lon_step = (xc2 - xc1) / (target_shape[1] - 1) if target_shape[1] > 1 else 0
target_lat_step = (yc2 - yc1) / (target_shape[0] - 1) if target_shape[0] > 1 else 0
target_lon_edges = np.linspace(xc1 - target_lon_step/2, xc2 + target_lon_step/2, target_shape[1] + 1)
target_lat_edges = np.linspace(yc1 - target_lat_step/2, yc2 + target_lat_step/2, target_shape[0] + 1)

# 计算每个目标像素与原始像素的重叠面积权重
# 利用numpy广播计算所有组合的交集
lon_left = np.maximum(lon_edges[:-1, np.newaxis], target_lon_edges[:-1])
lon_right = np.minimum(lon_edges[1:, np.newaxis], target_lon_edges[1:])
lon_overlap = np.maximum(0, lon_right - lon_left)

lat_bottom = np.maximum(lat_edges[:-1, np.newaxis], target_lat_edges[:-1])
lat_top = np.minimum(lat_edges[1:, np.newaxis], target_lat_edges[1:])
lat_overlap = np.maximum(0, lat_top - lat_bottom)

# 权重为重叠面积（原始像素面积为lon_step*lat_step，所以归一化到原始像素的占比）
weights = lon_overlap[np.newaxis, :, :] * lat_overlap[:, np.newaxis, :]
weights /= (lon_step * lat_step)

# 加权求和得到重采样后的数组
resampled_pop = np.tensordot(pop, weights, axes=((0, 1), (0, 1))).reshape(target_shape)

print(f"重采样后数据总和: {resampled_pop.sum()}")
print(f"总和误差（浮点运算导致）: {abs(resampled_pop.sum() - original_sum):.4f}")

这个方法能保证数据总和几乎完全不变（误差仅来自浮点运算），完美适配人口这类需要守恒的地理数据。

方法二：插值重采样（适合非守恒型数据） #

如果你的数据不需要严格保持总和（比如温度、高程这类连续型变量），可以直接用scipy.ndimage.zoom，它支持任意比例的插值重采样，用法很简单：

from scipy.ndimage import zoom

# 计算缩放比例：目标尺寸 / 原始尺寸
zoom_factor = (target_shape[0]/YSIZE, target_shape[1]/XSIZE)
# order=1是双线性插值，order=0是最近邻插值（适合离散型数据）
resampled_pop_zoom = zoom(pop, zoom_factor, order=1)

print(f"插值后数据总和: {resampled_pop_zoom.sum()}")

注意：插值方法会改变数据总和，因为它是基于邻域插值而非面积加权求和，所以只适合不需要守恒的场景。

再说说你的FFT方法为什么不行 #

fftpack.fft2(pop, shape=(60,80))的逻辑是：先对原始100×100数组做FFT转换到频域，然后直接把频域数组截断到60×80（丢弃大部分高频分量），再逆FFT转回空间域得到60×80的数组。这个操作本质是频域低通滤波+强制缩尺寸，完全没有考虑地理数据的空间分布逻辑，所以总和会大幅丢失，完全不适合这类重采样场景。

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