如何用 python 对离散数据做二维分箱？

灭火器

已知一系列离散点，用经纬度和变量值描述，可以保存为如下的二维表格：



二维分箱（2d binning）指规划好经纬度的正交网格，根据每个格子里落入的离散点计算格子的统计值，最后得到各个变量的格点数据，然后可以保存为 netCDF 文件。将非规则网格的 L2 卫星数据处理成格点数据时，经常会利用二维分箱的方法计算格子里的平均值。大概类似 NCL 官网这张图的效果：



NCL 里直接有现成的 bin_avg 函数，但 python 里要实现就更麻烦一些。

这里举出我常用的两种方法，就当抛砖引玉，想看看大家都是怎么实现的。简单粗暴的二维循环就忽略了。

方法一：scipy.stats.binned_statistic_2d，注意为了得到形如 (nlat, nlon) 的结果，需要颠倒 x 和 y 参数。

1. from scipy.stats import binned_statistic_2d
   
2. 
   
3. 
   
4. def binning2d(
   
5.     x: ArrayLike,
   
6.     y: ArrayLike,
   
7.     values: ArrayLike | list[ArrayLike],
   
8.     xbins: int | ArrayLike,
   
9.     ybins: int | ArrayLike,
   
10. ) -> tuple[NDArray, NDArray, NDArray]:
    
11.     def nanmean(arr: NDArray) -> float:
    
12.         arr = arr[~np.isnan(arr)]
    
13.         if arr.size == 0:
    
14.             return np.nan
    
15.         return arr.mean()
    
16. 
    
17.     binned, ybins, xbins, _ = binned_statistic_2d(
    
18.         y, x, values, bins=[ybins, xbins], statistic=nanmean
    
19.     )
    
20.     xlabels = get_bin_centers(xbins)
    
21.     ylabels = get_bin_centers(ybins)
    
22. 
    
23.     return xlabels, ylabels, binned
    
24. 
    
25. 
    
26. df = pd.read_csv("test.csv")
    
27. 
    
28. lon_bins, lon_labels = make_evenly_bins(-180, 180, 0.1)
    
29. lat_bins, lat_labels = make_evenly_bins(-90, 90, 0.1)
    
30. 
    
31. varnames = df.columns[2:]
    
32. _, _, binned = binning2d(
    
33.     x=df["lon"],
    
34.     y=df["lat"],
    
35.     values=[df[varname] for varname in varnames],
    
36.     xbins=lon_bins,
    
37.     ybins=lat_bins,
    
38. )
    
39. 
    
40. ds = xr.Dataset(
    
41.     {varname: (["lat", "lon"], data) for varname, data in zip(varnames, binned)},
    
42.     coords={"lon": lon_labels, "lat": lat_labels},
    
43. )

复制代码

方法二：pd.cut，更加简洁，速度也比 scipy 快很多。
更新：根据墨家大宝的回复修改了 observed 参数，并使用 reindex。

1. import numpy as np
   
2. import pandas as pd
   
3. import xarray as xr
   
4. from numpy.typing import ArrayLike, NDArray
   
5. 
   
6. 
   
7. def get_bin_centers(bins: ArrayLike) -> NDArray:
   
8.     bins = np.asarray(bins)
   
9.     assert len(bins) >= 2
   
10.     return (bins[1:] + bins[:-1]) / 2
    
11. 
    
12. 
    
13. def linspace2(x0: float, x1: float, dx: float) -> NDArray:
    
14.     assert dx > 0
    
15.     nx = int(abs(x1 - x0) / dx) + 1
    
16.     x = np.linspace(x0, x1, nx)
    
17. 
    
18.     return x
    
19. 
    
20. 
    
21. def make_evenly_bins(x0: float, x1: float, dx: float) -> tuple[NDArray, NDArray]:
    
22.     bins = linspace2(x0, x1, dx)
    
23.     labels = get_bin_centers(bins)
    
24. 
    
25.     return bins, labels
    
26. 
    
27. 
    
28. lon_bins, lon_labels = make_evenly_bins(-180, 180, 0.1)
    
29. lat_bins, lat_labels = make_evenly_bins(-90, 90, 0.1)
    
30. 
    
31. df = pd.read_csv("test.csv")
    
32. 
    
33. # 实际运算部分
    
34. ds = (
    
35.     df.assign(
    
36.         lon=pd.cut(df["lon"], lon_bins, labels=lon_labels, include_lowest=True),
    
37.         lat=pd.cut(df["lat"], lat_bins, labels=lat_labels, include_lowest=True),
    
38.     )
    
39.     .groupby(["lat", "lon"], observed=True)
    
40.     .mean()
    
41.     .to_xarray()
    
42.     .reindex(lon=lon_labels, lat=lat_labels)
    
43. )
    
44. 
    
45. ds.to_netcdf("test.nc")

复制代码

另外还有两种我尝试失败的 xarray 方法，看有没有坛友知道原因：

1. lon_cat = pd.cut(df['lon'], lon_bins, labels=lons, include_lowest=True)
   
2. lat_cat = pd.cut(df['lat'], lat_bins, labels=lats, include_lowest=True)
   
3. ds = df.to_xarray().assign(lon=('index', lon_cat), lat=('index', lat_cat)).groupby(['lat', 'lon']).mean()
   
4. 
   
5. # v2024.07.0 版本以上
   
6. from xarray.groupers import BinGrouper
   
7. lon_grouper = BinGrouper(lon_bins, labels=lons, include_lowest=True)
   
8. lat_grouper = BinGrouper(lat_bins, labels=lats, include_lowest=True)
   
9. ds = df.to_xarray().groupby(lon=lon_grouper, lat=lat_grouper).mean()

复制代码

墨家大宝

试试我的秘方：

df[['lon', 'lat']] = df[['lon', 'lat']] // 0.1 * 0.1 +0.05
da = df.groupby(['lon', 'lat']).mean().to_xarray().reindex(lat=lats, lon=lons)

内马尔

{:5_213:}{:5_213:}

edwardli

groupby即可。

edwardli

groupby即可。pandas.dataframe和xarray.dataset都有这个功能。

道理同此：  dailyTP.groupby('latitude').mean(dim="longitude")

【气小Py-007：中期预报降水趋势_2023.06.13-哔哩哔哩】【视频标记点 04:00】 https://b23.tv/p81Hsp8

灭火器

墨家大宝 发表于 2025-1-15 14:30
试试我的秘方：

df[['lon', 'lat']] = df[['lon', 'lat']] // 0.1 * 0.1 +0.05

学到了，确实是秘方，量化到 0.1 分辨率再加一半的分辨率，相当于得到了分箱的标签，然后做 groupby，虽然有些格子没有出现在结果里，但是最后通过 reindex 把这些格子（行和列）又补回来了，更快也更省内存。
不过秘方有两个假设：一是格子是等间距的，二是格子的左右边缘数值应该是分辨率的整数倍，比如说 70、70.1、70.2，而 70.05 就会产生一些误差。不过估计减去 minlon 再加上 minlon 就能解决。

灭火器

edwardli 发表于 2025-1-15 15:42
groupby即可。pandas.dataframe和xarray.dataset都有这个功能。

道理同此：  dailyTP.groupby('latitude ...

看了一下视频，dailyTP 是时间、经度和纬度三个维度互相独立的网格数据吧，这个帖子想讨论的是经度和纬度不独立，二者会同时变化的一维散点数据，例如 flatten 后的卫星 L2 网格那种。
我试了下连续两次 groupby_bins，结果会是一维的，不符合预期。

墨家大宝

灭火器 发表于 2025-1-15 21:56
学到了，确实是秘方，量化到 0.1 分辨率再加一半的分辨率，相当于得到了分箱的标签，然后做 groupby，虽 ...

我凭印象盲打的，第二行里lons和lats一般我都会先处理成每个bin的中心便于画图。浮点数和非等距问题可以通过把第一行换成pd.cut解决，只不过对于均匀格点，尤其是1度等整数均匀格点，手动整除取整比cut更方便。

楼主有写公众号吗？方便的话加上我这个方法让我转载一下呗，提到“墨大宝”就行，我的公众号已经空了很久了。。。

任性小王子

墨家大宝 发表于 2025-1-16 09:13
我凭印象盲打的，第二行里lons和lats一般我都会先处理成每个bin的中心便于画图。浮点数和非等距问题可以 ...

学习到了，非常感谢大佬们

灭火器

墨家大宝 发表于 2025-1-16 09:13
我凭印象盲打的，第二行里lons和lats一般我都会先处理成每个bin的中心便于画图。浮点数和非等距问题可以 ...

无公众号，可以给你投稿