如何在Python中对空间网格内的空间点数据求平均值

我有一组包含坐标及对应数值的数据点，希望将这些数据点划分至空间网格中，对每个网格内的所有点计算平均值。以下是我实现该功能的代码，想请教是否存在更简洁高效的实现方式。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

# ***********************************************
# *** 获取1°x1°经纬度网格内的均值数据 ***
    
# 数据集
lats =  np.array((40.1, 42.1, 40.7, 43.1, 41.0, 42.7, 40.7, 43.1, 41.2, 42.7, 40.5, 43.1, 40.2))
lons =  np.array(( 1.2,  2.6,  3.1,  1.7,  2.1,  1.9,  3.4,  2.7,  3.4,  1.3,  2.1,  2.9,  3.7))
temps = np.array((28.3, 25.1, 24.3, 27.4, 26.3, 21.3, 25.0, 21.9, 22.8, 21.7, 23.8, 21.5, 20.1))

# 初始化网格与统计数组
latbins = np.arange(-90, 90.1, 1)
lonbins = np.arange(0, 360.1, 1)
sum_temp = np.zeros((latbins.size-1, lonbins.size-1))
occ_temp = np.zeros((latbins.size-1, lonbins.size-1))

# 逐点累加网格内的温度总和与观测数
for i, temp in enumerate(temps):
    # 手动查找当前点所属的网格索引
    i_latbin = 0
    while lats[i] > latbins[i_latbin]:
        i_latbin += 1
    i_lonbin = 0
    while lons[i] > lonbins[i_lonbin]:
        i_lonbin += 1  
    
    # 累加数据
    sum_temp[i_latbin-1, i_lonbin-1] += temp
    occ_temp[i_latbin-1, i_lonbin-1] += 1

# 计算均值并处理空网格
FILL_VALUE = -9999.
mean_temp = np.ones((latbins.size-1, lonbins.size-1))*FILL_VALUE
mean_temp[occ_temp != 0] = sum_temp[occ_temp != 0]/occ_temp[occ_temp != 0]
mean_temp = np.ma.masked_where(mean_temp==FILL_VALUE, mean_temp)


# ****************
# *** 绘制可视化图 ***

plt.figure(figsize=(5, 10))

# 绘制温度均值图
ax = plt.subplot(211)
ax.set_facecolor('0.8')
plt.pcolormesh(lonbins, latbins, mean_temp, cmap=cm.RdBu_r)
plt.colorbar(label='温度 (°C)')
plt.xlim(1, 4)
plt.ylim(40, 44)
plt.xlabel('经度')
plt.ylabel('纬度')

# 绘制观测数图
ax = plt.subplot(212)
plt.pcolormesh(lonbins, latbins, occ_temp)
plt.colorbar(label='观测数量')
plt.xlim(1, 4)
plt.ylim(40, 44)
plt.xlabel('经度')
plt.ylabel('纬度')

plt.show()

更简洁高效的实现方式 #

利用numpy的向量化操作替代手动循环，能大幅提升代码效率和可读性，尤其是数据量较大时效果更明显：

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

# 数据集
lats = np.array((40.1, 42.1, 40.7, 43.1, 41.0, 42.7, 40.7, 43.1, 41.2, 42.7, 40.5, 43.1, 40.2))
lons = np.array((1.2, 2.6, 3.1, 1.7, 2.1, 1.9, 3.4, 2.7, 3.4, 1.3, 2.1, 2.9, 3.7))
temps = np.array((28.3, 25.1, 24.3, 27.4, 26.3, 21.3, 25.0, 21.9, 22.8, 21.7, 23.8, 21.5, 20.1))

# 定义网格区间
latbins = np.arange(-90, 90.1, 1)
lonbins = np.arange(0, 360.1, 1)

# 批量获取每个点所属网格的索引
lat_indices = np.digitize(lats, latbins) - 1
lon_indices = np.digitize(lons, lonbins) - 1

# 一次性计算每个网格的观测数和温度总和
occ_temp, _, _ = np.histogram2d(lats, lons, bins=[latbins, lonbins])
sum_temp, _, _ = np.histogram2d(lats, lons, bins=[latbins, lonbins], weights=temps)

# 计算均值并处理空网格
FILL_VALUE = -9999.
mean_temp = np.full(occ_temp.shape, FILL_VALUE)
valid_mask = occ_temp != 0
mean_temp[valid_mask] = sum_temp[valid_mask] / occ_temp[valid_mask]
mean_temp = np.ma.masked_where(mean_temp == FILL_VALUE, mean_temp)

# 绘图逻辑与原代码一致
plt.figure(figsize=(5, 10))

# 绘制温度均值图
ax = plt.subplot(211)
ax.set_facecolor('0.8')
plt.pcolormesh(lonbins, latbins, mean_temp, cmap=cm.RdBu_r)
plt.colorbar(label='温度 (°C)')
plt.xlim(1, 4)
plt.ylim(40, 44)
plt.xlabel('经度')
plt.ylabel('纬度')

# 绘制观测数图
ax = plt.subplot(212)
plt.pcolormesh(lonbins, latbins, occ_temp)
plt.colorbar(label='观测数量')
plt.xlim(1, 4)
plt.ylim(40, 44)
plt.xlabel('经度')
plt.ylabel('纬度')

plt.show()

核心优化点 #

• 替换手动索引循环：用np.digitize批量生成所有点的网格索引，避免逐点循环判断，效率提升显著。
• 向量化统计计算：通过np.histogram2d同时完成网格内的计数和加权求和，一步到位替代逐点累加，充分利用numpy的底层优化。
• 简化空值处理：用np.full初始化填充值，结合掩码数组让空网格的处理逻辑更清晰简洁。

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