当前位置：首页 > article >正文

Python散点统计栅格化

article 2025/1/24 16:46:47

上周在气象家园论坛有一篇《如何用 python 对离散数据做二维分箱？》的帖子，刚好我有一些相关经验就参与了讨论，在此做一个详细的分享。

在该帖中，作者灭火器描述

二维分箱（2d binning）指规划好经纬度的正交网格，根据每个格子里落入的离散点计算格子的统计值，最后得到各个变量的格点数据……将非规则网格的L2卫星数据处理成格点数据时，经常会利用二维分箱的方法计算格子里的平均值。

我的相关实践有：

极轨卫星的甲烷数据栅格化为分辨率为1°格点
闪电定位仪数据绘制密度分布地图

我更愿意把这个“二维分箱”称为“散点统计栅格化”，和“散点插值栅格化”类似，都是从不规则离散点数据处理获得规则栅格数据。

插值栅格化后格点值可能受格点外散点影响，而统计栅格化后格点值只受处于该格点内散点的影响。可以看出统计栅格化对距离的处理比较粗糙，应当更好实现，也许正是因为简单所以并没有像插值工具那样繁荣的工具生态。

接下来介绍一下我总结的散点统计栅格化Python方案，看看有多简单。核心依赖只有pandas和xarray包。

示例

有离散点数据df（此处为构造的假数据）：

import numpy as np
import pandas as pd

θ = np.random.uniform(0, 2 * np.pi, size=10000)
r = np.random.uniform(0, 1, size=10000) ** 0.6
x1 = r * np.cos(θ) * 10 + 125
y1 = r * np.sin(θ) * 20 + 30
θ = np.random.uniform(0, 2 * np.pi, size=10000)
r = np.random.uniform(0, 1, size=10000) ** 0.6
x2 = r * np.cos(θ) * 15 + 85
y2 = r * np.sin(θ) * 10 + 30
df = pd.DataFrame({
    'lon': np.concatenate([x1, x2]),
    'lat': np.concatenate([y1, y2])
})
df['value'] = np.hypot(df['lon'] - 105, df['lat'] - 30)

若希望以1°分辨率计算散点平均值栅格化：

df_mean = df.copy()
df_mean[['lat', 'lon']] = df_mean[['lat', 'lon']].round().astype(np.int32)
da_mean = df_mean.groupby(['lat', 'lon'])['value'].mean()\
    .to_xarray().reindex(
        lat=range(0, 61), lon=range(70, 141)
    )

详解

df_mean = df.copy()

复制一份是为了保留原始数据，以防后续需要，可看情况省略。

df_mean[['lat', 'lon']] = df_mean[['lat', 'lon']].round().astype(np.int32)

把坐标化为整数，这一步就已经把每个散点属于哪个栅格标记完成。

.astype(np.int32)是为了把数值类型变成整型，以回避浮点数应当相等却不等的问题。

da_mean = df_mean.groupby(['lat', 'lon'])['value'].mean()\
    .to_xarray().reindex(
        lat=range(0, 61), lon=range(70, 141)
    )

其中df_mean.groupby(['lat', 'lon'])完成了对标记为相同格点散点的分组。

['value'].mean()是对分完组后散点的value求平均值。

.to_xarray()是将统计完的一行一行的DataFrame按分组的lat和lon变形为二维的DataArray。

.reindex(lat=range(0, 61), lon=range(70, 141))是为了补齐栅格，因为如果某整行或整列都没有散点，该行或列是不会出现在.to_xarray()的DataArray中的。

扩展

统计方式

上面这个例子统计的是平均值，若df是闪电定位仪之类的散点数据，统计栅格内散点数量（即密度，单位为个/度^2），只要把.mean()方法改为.count()即可：

df_count = df.copy()
df_count[['lat', 'lon']] = df_count[['lat', 'lon']].round().astype(np.int32)
da_count = df_int.groupby(['lat', 'lon'])['value'].count()\
    .to_xarray().reindex(
        lat=range(0, 61), lon=range(70, 141)
    )

同样可以按需换成其他统计方法比如标准差.std()等，非常规统计量可以使用.apply()方法自定义。

缺失值

不包含任何散点的格点最后都被标记为缺失值np.nan，想填补缺失值的话加上.fiilna()方法即可，以补零为例：

df_count = df.copy()
df_count[['lat', 'lon']] = df_count[['lat', 'lon']].round().astype(np.int32)
da_count = df_int.groupby(['lat', 'lon'])['value'].count()\
    .to_xarray().reindex(
        lat=range(0, 61), lon=range(70, 141)
    ).fillna(0)

非整数分辨率

如果希望分辨率是0.3度，稍作修改：

w, e, s, n = 70, 140, 0, 60
res = 0.3
df_count = df.copy()
df_count['lat'] = ((df_count['lat'] - s - res / 2) / res).round().astype(np.int32)
df_count['lon'] = ((df_count['lon'] - w - res / 2) / res).round().astype(np.int32)
da_count = df_count.groupby(['lat', 'lon'])['value'].count()\
    .to_xarray().reindex(
        lat=range(int((n - s) // res + 1)), lon=range(int((e - w) // res + 1))
    )
da_count['lat'] = da_count['lat'] * res + s
da_count['lon'] = da_count['lon'] * res + w

思路仍然是整数化标记散点所属栅格，最后将整数标记恢复成坐标而已。

非均匀格点

虽然不常见，但是比非整数分辨率更泛化的情况是非均匀栅格格点，需要用到pandas.cut方法：

lon = np.array([70, 71, 90.5, 95, 100, 112, 130, 136, 140])
lat = np.array([0, 3, 19.5, 20.8, 35, 50, 60])
df_count = df.copy()
df_count['lat'] = pd.cut(df_count['lat'], bins=lat, labels=range(len(lat) - 1))
df_count['lon'] = pd.cut(df_count['lon'], bins=lon, labels=range(len(lon) - 1))
da_count = df_count.groupby(['lat', 'lon'], observed=True)['value'].count()\
    .to_xarray().reindex(
        lat=range(len(lat) - 1), lon=range(len(lon) - 1)
    )
da_count['lat'] = (lat[:-1] + lat[1:]) / 2
da_count['lon'] = (lon[:-1] + lon[1:]) / 2