应用场景

现有以csv格式存储的千万级别点位信息，每个点包含数值类属性值，需将这些点按一定网格规格进行统计，计算各网格内点的统计结果。

解决思路

这是一个典型的空间统计问题，主要思路是通过空间连接建立点和网格的一对一连接，然后按照网格进行分组统计计算即可。对于小数据量的输入，可在ArcGIS中通过Identity-Spatial join工具连用解决，但对于海量数据，在ArcMap中处理的效率就较低了。

由于点、面关联的空间计算相对简单，因此对于大量数据，考虑采用GeoPandas或PostGIS内置函数解决，本文简单阐述以上两种方法的思路及实现。

实现方案

GeoPandas

数据准备

点位数据，格式：csv
应用ArcMap的Fishnet工具生成的网格，格式：polygon类shp，可参考之前文章。

点位数据约1000w条，网格数为23w。

环境准备

在Python环境中安装GeoPandas包，采用官方推荐的conda install --channel conda-forge geopandas方式容易出现各种问题，因此用pip进行安装。
安装过程如下^[1]

根据Python环境版本，在网站下载GeoPandas的四个依赖包gdal、Shapely、Fiona、pyproj的whl文件。
切换到以上文件的下载目录，通过pip install xxx.whl方式依照gdal、Shapely、Fiona、pyproj的顺序进行安装，安装完成后pip install geopandas即可

代码实现

首先通过GeoPandas提供的points_from_xy函数构造geometry并生成点要素的GeoDataFrame
读取网格要素载入GeoDataFrame
通过sjoin函数实现空间连接^[2]，空间拓扑关系参考
通过groupby进行属性统计计算

import pandas as pd
import geopandas as gpd
import pyproj
import time
import os

# define dir
root_dir = r'./Dir/'
file_list = os.listdir(root_dir)

# Month Iterator
for file in file_list:
    # get filename 
    filename = file.split('.')[0]

    time_start=time.time()
    # load point data
    raw_point= pd.read_csv(root_dir + file,encoding='utf-8')
    # numeric format 
    raw_point[['longitude','latitude']] = raw_point[['longitude','latitude']].apply(pd.to_numeric)
    # subset
    raw_point = raw_point[['longitude', 'latitude','SOx_emission', 'NOx_emission', 'PM10_emission','PM2P5_emission', 'hc_emission', 'CO_emission']]
    print('load data done')
    
    # 通过points_from_xy将经纬度信息生成点要素
    vessel_point = gpd.GeoDataFrame(raw_point, geometry=gpd.points_from_xy(raw_point.longitude, raw_point.latitude))
    # define projection
    vessel_point.crs = pyproj.CRS.from_user_input('EPSG:4326') 
    # 导出点shp文件
    # vessel_point.to_file(root_dir +  filename + '.shp', driver='ESRI Shapefile',encoding='utf-8')
    print('make gdf done')
    
    # load grid data
    vessel_grid = gpd.read_file(r'./result/grid.shp')
    vessel_grid.crs = pyproj.CRS.from_user_input('EPSG:4326')
    # 空间连接
    vessel_merged = gpd.sjoin(vessel_grid, vessel_point, predicate='contains')
    print('merge done')
    
    # subset
    vessel_merged2 = vessel_merged[['Id', 'SOx_emission', 'NOx_emission', 'PM10_emission', 'PM2P5_emission',
           'hc_emission', 'CO_emission']]
    # 按网格求和
    vessel_grid_merged = vessel_merged2.groupby('Id').agg('sum')
    # 将计算结果与geometry网格进行关联，用于输出
    vessel_grid_shp = vessel_grid.merge(vessel_grid_merged, on='Id')
    print('stat done')
    
    # output shp file
    vessel_grid_shp.to_file(root_dir +  filename + '.shp', driver='ESRI Shapefile',encoding='utf-8')
    vessel_grid_merged.to_csv(root_dir +  filename + 'grid.csv')
    
    # calculate running time
    time_end=time.time()
    print('totally cost',time_end-time_start)

最终效果如下：

PostGIS

数据准备

点位数据，格式：csv

环境搭建

采用Docker部署PostGIS及PgAdmin环境

实现思路

使用PostGIS时，不需要实际生成网格并存储，而是通过生成器动态创建网格然后与其他空间数据做叠加统计分析^[3]。

将点位shp导入Pg数据库
根据vessel_point范围使用ST_SquareGrid生成动态网格
INNER JOIN进行空间连接汇总

SELECT sum(NOx_emission) as NOx_emission, squares.geom
FROM
    ST_SquareGrid(
        ST_SetSRID(ST_EstimatedExtent('vessel_point', 'geom'), 4326)
    ) AS squares
    INNER JOIN
    vessel_point AS p
    ON ST_Contains(p.geom, squares.geom)
GROUP BY squares.geom;