使用numpy histogram计算SST随TCWV平均时结果异常的问题排查

我尝试用numpy的histogram函数，以可降水量（TCWV）为自变量，计算对应分组的海表温度（SST）平均值，核心代码如下：

# 计算分箱内的SST平均值
num, _   = np.histogram(tcwv, bins=bins)
sstsum, _ = np.histogram(tcwv, bins=bins,weights=sst)
out=np.zeros_like(sstsum)
out[:]=np.nan
sstav  = np.divide(sstsum,num,out=out, where=num>100)

完整可复现代码如下，但绘制原始数据散点图和计算出的平均值点时，发现平均值点远偏离原始数据分布，找不到原因，怀疑是舍入误差导致？

import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset

# 若netcdf库版本较新可直接访问
url = ('http://clima-dods.ictp.it/Users/tompkins/CRM/data/WRF_1min_mem3_grid4.nc#mode=bytes')
ds=Dataset(url)

### 否则需下载后使用：
###ifile="WRF_1min_mem3_grid4.nc"
###ds=Dataset(idir+ifile)


# 轴的分箱
bins=np.linspace(40,80,21)

iran1,iran2=40,60

# 可存入字典循环处理
sst=ds.variables["sst"][iran1:iran2+1,:,:]
tcwv=ds.variables["tcwv"][iran1:iran2+1,:,:]

# 不需要展平，只是尝试看看是否有用（没用）
sst=sst.flatten()
tcwv=tcwv.flatten()

# 计算分箱内的SST平均值
num, _   = np.histogram(tcwv, bins=bins)
sstsum, _ = np.histogram(tcwv, bins=bins,weights=sst)
out=np.zeros_like(sstsum)
out[:]=np.nan
sstav  = np.divide(sstsum,num,out=out,where=num>100)

# 分箱中心值
avbins=(np.array(bins[1:])+np.array(bins[:-1]))/2

#绘图
subsam=2
fig,(ax)=plt.subplots()
plt.scatter(tcwv.flatten()[::subsam],sst.flatten()[::subsam],s=0.05,marker=".")
plt.scatter(avbins,sstav,s=3,color="red")
plt.ylim(299,303)
plt.savefig("scatter.png")

1. 优先排查异常值干扰 #

原始数据中可能存在超出绘图范围的极端SST值，这些值会拉高/拉低分箱均值，但散点图因ylim(299,303)被隐藏，导致均值看起来偏离。先过滤异常值：

# 过滤SST和TCWV的空值与异常值
mask = (sst >= 299) & (sst <= 303) & (~np.isnan(sst)) & (~np.isnan(tcwv))
sst = sst[mask]
tcwv = tcwv[mask]

2. 验证分箱计算逻辑 #

用手动计算对比numpy结果，确认分组逻辑是否正确：

# 手动计算第一个分箱的均值验证
bin_mask = (tcwv >= bins[0]) & (tcwv < bins[1])
print("手动计算均值：", sst[bin_mask].mean())
print("numpy计算均值：", sstav[0])

如果两者结果一致，说明分箱逻辑没问题，问题出在数据本身；如果不一致，检查bins定义是否存在边界错误。

3. 用pandas交叉验证 #

用pandas的分组功能替代numpy histogram，验证结果是否一致：

import pandas as pd
df = pd.DataFrame({"tcwv": tcwv, "sst": sst})
df["bin"] = pd.cut(df["tcwv"], bins=bins, labels=avbins)
sstav_pd = df.groupby("bin")["sst"].mean().values

# 对比两种方法的结果
print("numpy结果：", sstav)
print("pandas结果：", sstav_pd)

若结果差异大，说明numpy代码存在逻辑疏漏；若一致，说明原始数据分布本身存在均值偏离散点集中区的情况（比如分箱内数据分布极不均衡）。

4. 舍入误差排除 #

舍入误差不会导致均值大幅偏离——SST（300量级）和TCWV（60量级）的数值差异不会引发精度问题，无需考虑此因素。

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