离散点数据的网格化统计

很久没有发文了，在破站上看到某up讲的经纬度网格计数的视频气小py-001：经纬度网格计数_2023.04.25_哔哩哔哩_bilibili，正好想延伸下思路，本人也是学气象的，经常会遇到很多类似的离散数据的网格化处理的问题，比如离散的降雨站点与数值预报的网格配对，比如雷达的基数据的空间网格统计等，记得多年前为如何处理这种插值或者统计到处咨询，后面还是自己解决了这个问题，现在分享给大家。

我的解决思路是也是先创建一个空白的网格，然后在网格中进行统计计算，但像视频中那样双重循环的话在数据量达几万或几十万后会很慢，所以换了另外一个方法，就是在创建空白网格的时候先对x，y坐标先进行标准化定位，作为新的xgrid/ygrid列，然后通过groupby函数对新的xgrid/ygrid聚合，聚合函数可以直接使用自带的，也可以自定义，最后透过透视表的方法形成新的标准网格，分解代码示例如下：

import numpy as np
import pandas as pd

# 创建随机数据
lon_list = []
lat_list = []
data_list = []
for i in range(1000):
    lon = np.random.uniform(101, 104)
    lat = np.random.uniform(28.5, 31)
    data = np.random.randint(1,50)
    lon_list.append(lon)
    lat_list.append(lat)
    data =data_list.append(data)
df = pd.DataFrame({"lon":lon_list,"lat":lat_list, "rain":data_list})

# 增加两列的网格坐标
xGrid = yGrid = 0.25 # 分辨率
dataXYV = df.values

# 对原始数据的经纬度进行规范化投射
dataGridX = np.floor(dataXYV[:, 0] / xGrid) * xGrid + xGrid / 2.0
dataGridY = np.floor(dataXYV[:, 1] / yGrid) * yGrid + yGrid / 2.0
datanew = np.c_[dataXYV, dataGridX, dataGridY]
# 创建表格型数据,分别为x,y,z,dbz,xgrid,ygrid
df_new = pd.DataFrame({
    'x': datanew[:, 0],
    'y': datanew[:, 1],   
    'data': datanew[:, 2],
    'xgrid': datanew[:, 3],
    'ygrid': datanew[:, 4]
})
# 聚合统计
dataGrid = df_new[:].groupby(['xgrid','ygrid']).count()
dataGrid 

# 生成透视表
dataGrid.reset_index().pivot(index="ygrid", columns="xgrid", values="data")

 注意：按照我之前的算法，生成的透视表网格坐标其实是在四个顶点的中心，所以使用的时候需要注意。当然也可以自己定义，还有如果该网格中无数据，会标记为Nan。

由于pandas的常用聚合函数运行速度很高，这样就算几十万行的数据，投射到网格速度也很快。用来处理雷达的基数据组网时很好用。

下面是我自己用来雷达数据组合反射率统计的函数：

def scatter2Grid(dataXYZV, xGrid, yGrid, method="max", z_range=None):
    """
    输入的XYZV格式的4维散点numpy格式数据，最后一列放数值，根据网格将散点投射到粗网格的中心点上
    xGrid为X方向的网格长度,yGrid为y方向的网格长度,z_range为垂直方向上的取值范围用来裁剪高度，默认全部
    method，默认取最大值
    """  
    # 对原始数据的经纬度进行规范化投射,默认投射到中心点
    dataGridX = np.floor(dataXYZV[:, 0] / xGrid) * xGrid + xGrid / 2.0
    dataGridY = np.floor(dataXYZV[:, 1] / yGrid) * yGrid + yGrid / 2.0
    #  创建表格型数据,分别为x,y,z,dbz,xgrid,ygrid
    df_new = pd.DataFrame({
        'x': dataXYZV[:, 0],
        'y': dataXYZV[:, 1],
        'z': dataXYZV[:, 2],  # 这一列如果没有的话，可以不用管它
        'data': dataXYZV[:, -1],
        'xgrid': dataGridX,
        'ygrid': dataGridY
    })

    if z_range:
        df_new.loc[(df_new['z'] < z_range[0]) | (df_new['z'] > z_range[1]),
                   "data"] = 0
    # 对表格数据聚类分组,如果都在一个网格内，根据mode的值采取不同的处理
    dataGrid = pd.DataFrame()
    if method == "mean":  # 取平均值
        dataGrid = df_new[:].groupby(['xgrid', 'ygrid']).mean()
    elif method == "min":  # 取最小值
        dataGrid = df_new[:].groupby(['xgrid', 'ygrid']).min()
    elif method == "max":  # 取最大值
        dataGrid = df_new[:].groupby(['xgrid', 'ygrid']).max()
    elif method == "median":  # 取最大值
        dataGrid = df_new[:].groupby(['xgrid', 'ygrid']).median()
    elif method == "count":  # 取最大值
        dataGrid = df_new[:].groupby(['xgrid', 'ygrid']).count()
    else:
        # 这一部分为自定义聚合函数，这个就自己写了
        pass
    df_Grid = dataGrid.reset_index().pivot(index="ygrid",
                                           columns="xgrid",
                                           values="data")
    return df_Grid

以下是测试，生成500万个离散点数据，花了18.9秒。

lon_list = []
lat_list = []
data_list = []
for i in range(5000000):
    lon = np.random.uniform(101, 104)
    lat = np.random.uniform(28.5, 31)
    data = np.random.randint(1,50)
    lon_list.append(lon)
    lat_list.append(lat)
    data =data_list.append(data)
df = pd.DataFrame({"lon":lon_list,"lat":lat_list, "rain":data_list})
df

 以下是统计花时，412毫秒，算很快了。

xGrid = yGrid = 0.01  # 分辨率
dataXYV = df.values
method = "count"  # 统计方法标识

df_Grid = scatter2Grid(dataXYZV=df.values,
                       xGrid=xGrid,
                       yGrid=yGrid,
                       method=method)
df_Grid