Python中基于时间维度统计网格单元内地理点数的方法

场景与现有代码 #

用户拥有包含经纬度、日期和置信度的CSV数据，希望创建地理网格并按时间维度统计每个网格内的数据点数。现有GeoPandas代码实现了基础的网格计数，但无法完成时间维度的分组统计。

示例CSV数据 #

latitude  longitude    acq_date  confidence
0        -8.1135   112.9281  2001-01-01          99
1        -6.4586   143.2235  2001-01-03          86
2         6.6564   125.0055  2001-01-03          85
3         6.6545   124.9990  2001-01-03          84
4         9.7481   107.9814  2001-01-03          96
...          ...        ...         ...         ...
456844   -4.0529   143.4047  2020-12-28          89
456845   -8.1128   112.9365  2020-12-30         100
456846   -2.5768   121.3746  2020-12-31         100
456847   -2.5754   121.3848  2020-12-31          84
456848   -1.4573   127.4369  2020-12-31          90

用户现有代码 #

# convert df into a geopandas geodataframe
gdf = geopandas.GeoDataFrame(df,
                             geometry=geopandas.points_from_xy(
                                 df.longitude, df.latitude),
                             crs='epsg:4326')
# gdf.head()
# drop lon lat
gdf = gdf.drop(columns=['longitude', 'latitude'])

# total area for the grid
xmin=93
ymin=-11
xmax=141
ymax=8
# how many cells across and down
n_cells = 104
cell_size = (xmax-xmin)/n_cells
# projection of the grid
crs = 'epsg:4326'

# create the cells in a loop
grid_cells = []
for x0 in np.arange(xmin, xmax+cell_size, cell_size):
    for y0 in np.arange(ymin, ymax+cell_size, cell_size):
        # bounds
        x1 = x0-cell_size
        y1 = y0+cell_size
        grid_cells.append(box(x0, y0, x1, y1))
cell = geopandas.GeoDataFrame(grid_cells, columns=['geometry'], crs=crs)

merged = geopandas.sjoin(gdf, cell, how='left', predicate='within')
# make a simple count variable that we can sum
merged['n_fires'] = 1
# Compute stats per grid cell -- aggregate fires to grid cells with dissolve
dissolve = merged.dissolve(by="index_right", aggfunc="count")
# put this into cell
cell.loc[dissolve.index, 'n_fires'] = dissolve.n_fires.values

解决方案 #

核心思路是同时按网格索引和时间维度分组聚合，以下是具体修改步骤：

1. 预处理时间列 #

先将日期列转为datetime类型，再根据需求提取时间粒度（年、月、日等）：

import pandas as pd
import geopandas as gpd
from shapely.geometry import box
import numpy as np

# 转换时间列为datetime类型
gdf['acq_date'] = pd.to_datetime(gdf['acq_date'])
# 按年份统计（可替换为其他粒度：如按年月用dt.to_period('M')，按日用dt.date）
gdf['year'] = gdf['acq_date'].dt.year

2. 修改聚合逻辑 #

在聚合时，同时按网格索引（index_right）和时间列分组，替代原有的单一网格索引分组：

merged = gpd.sjoin(gdf, cell, how='left', predicate='within')
merged['n_fires'] = 1

# 按网格索引和时间维度分组统计数据点数
dissolve = merged.dissolve(by=['index_right', 'year'], aggfunc={'n_fires': 'sum'}).reset_index()

# 数据量较大时，用groupby性能更优
# dissolve = merged.groupby(['index_right', 'year'])['n_fires'].sum().reset_index()

3. 关联统计结果到网格 #

根据需求选择长格式或宽格式输出：

方式1：长格式（每个网格+时间组合为一行）

# 合并统计结果到网格GeoDataFrame
cell_time_stats = cell.merge(dissolve, left_index=True, right_on='index_right', how='left')
# 填充空值为0（表示该网格对应时间内无数据点）
cell_time_stats['n_fires'] = cell_time_stats['n_fires'].fillna(0)

方式2：宽格式（每个网格为一行，不同时间为列）

# 透视表转换为宽格式
pivot_df = dissolve.pivot(index='index_right', columns='year', values='n_fires').fillna(0)
# 合并到网格
cell_wide_stats = cell.merge(pivot_df, left_index=True, right_index=True, how='left')

4. 完整示例代码 #

import pandas as pd
import geopandas as gpd
from shapely.geometry import box
import numpy as np

# 读取CSV数据
df = pd.read_csv('your_data.csv')

# 转换为GeoDataFrame
gdf = gpd.GeoDataFrame(df,
                       geometry=gpd.points_from_xy(df.longitude, df.latitude),
                       crs='epsg:4326')
gdf = gdf.drop(columns=['longitude', 'latitude'])

# 预处理时间列
gdf['acq_date'] = pd.to_datetime(gdf['acq_date'])
gdf['year'] = gdf['acq_date'].dt.year

# 创建网格
xmin=93
ymin=-11
xmax=141
ymax=8
n_cells = 104
cell_size = (xmax-xmin)/n_cells
crs = 'epsg:4326'

grid_cells = []
for x0 in np.arange(xmin, xmax+cell_size, cell_size):
    for y0 in np.arange(ymin, ymax+cell_size, cell_size):
        x1 = x0-cell_size
        y1 = y0+cell_size
        grid_cells.append(box(x0, y0, x1, y1))
cell = gpd.GeoDataFrame(grid_cells, columns=['geometry'], crs=crs)

# 空间连接与聚合
merged = gpd.sjoin(gdf, cell, how='left', predicate='within')
merged['n_fires'] = 1
dissolve = merged.groupby(['index_right', 'year'])['n_fires'].sum().reset_index()

# 生成长格式统计结果
cell_time_stats = cell.merge(dissolve, left_index=True, right_on='index_right', how='left')
cell_time_stats['n_fires'] = cell_time_stats['n_fires'].fillna(0)

# 查看结果
print(cell_time_stats.head())

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