按照多边形（如格网/行政区划）边界汇总统计其内的点/栅格是空间分析中的常用操作（栅格实际上可以看做点数据）。
对于软件，如果点是矢量，则可以直接通过ArcMap/QGIS的spatial join计算，点为栅格，则称zonal statistic……实际中数据多即需要批处理，可以采用QGIS、ArcPy或者基于Python（如geopandas rasterstats等模块）……
具体的实现流程视数据情况而定，数据量大的时候可以考虑看能不能先分块再算或者多分几步，这里仅给出一些我的处理案例，数据量更大的可能就需要考虑进一步优化实现思路了。此外：ArcGIS的底层是32位貌似且Python还是2.7的，因此运行会更慢，当然现在慢慢迁移到ArcGIS Pro，但是Pro比较吃电脑配置很多操作我更加推荐开源的QGIS，一方面QGIS参数的控制和结果更加细致，且运行速度（Python3.6）和结果预览都更快！
∑栅格×矢量zonal
以统计格网(polygons)内栅格值为例：
Python
- ArcPy，实现Zonal Statistic功能，应该在空间分析工具下，可自行查看文档
- 一个专门用于做分区统计的模块rasterstats文档，较推荐，这里给一个简单的个人案例，统计区划内的高程信息。经过简单测试发现时间为50秒，而在QGIS中运行耗时150秒，ArcGIS理论上更慢，说明该工具还可以（毕竟底层C++）。
- 传入统计的矢量数据或geometry类型+栅格图层文件路径即可，输出一个字典，可以直接转成df
- 也可以自定义计算函数作为统计参数stats
- 对于分类的栅格数据可以添加
categorical = True参数
1
2
3
4
5
6
7
8
9
10
11
12
| from rasterstats import zonal_stats
import pandas as pd
import geopandas as gpd
gdf = gpd.read_file("./青海西藏县.shp", encoding = 'utf-8')
tiffile = "A:/Data/dem.tif"
stats = zonal_stats(gdf,demfile,stats=['min', 'max', 'median', 'majority', 'sum'])
stats2 = zonal_stats("./青海西藏县.shp",demfile,stats='count')
res = pd.concat([gdf, pd.DataFrame(stats)], axis=1)
|
- 其他暂不探索
软件
之前汇总过一次全国960w格网的栅格计算，分省计算的，基于的是QGIS的批处理（通过导入参数json文件实现）/也尝试了其Python脚本发现不好用，可能是我还没上手习惯
ArcGIS略
QGIS：在工具箱搜索zonal或者分区统计即可
QGIS批处理，QGIS的脚本不如ArcPy好用，所以批处理之前是用参数导入的方式实现的，但是其参数设置还需要经过一道转换，略显麻烦，但也还行，代码见：
参数设置代码参考，其要求比较严格需要标准的路径、json对象和json文件等，当时debug了好一会~，json的介绍可以我写过的这里，需要设置输入文件、输出文件、栅格波段、计算方式（数量、均值等）、前缀等等
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
| from glob import glob
import json
inpath = "A:/Data/ChinaGrid/grids_v1/shp"
infiles = glob(inpath+"/*.shp")
infiles = [inpath+"/"+os.path.basename(file) for file in infiles]
outfiles = [infile.replace('v1','v2') for infile in infiles]
paras = []
for i in range(len(outfiles)):
infile = json.dumps(infiles[i])
inraster = json.dumps("A:/Data/ChinaGrid/base/china-gaia.tif")
prefix = json.dumps('imp_')
statis = json.dumps([1])
para = {
"PARAMETERS":{
"INPUT":infile,
"INPUT_RASTER":inraster,
"RASTER_BAND":'1',
"COLUMN_PREFIX": prefix,
"STATISTICS": statis},
"OUTPUTS": {"OUTPUT": outfiles[i]}
}
paras.append(para)
print(paras[0])
with open("A:\Data\ChinaGrid\zonal_batch.json",'w',encoding='utf-8') as f:
json.dump(paras,f)
|
∑点×矢量sjoin
以统计格网(polygons)内点的数量为例，比如我这里统计8w个0.05°格网内微博签到点的数量：
软件
- ArcMap：右击多边形直接join汇总点，或者通过ArcToolbox~Tools.tbxJoin
- QGIS：按位置连接，在QGIS工具箱搜索join或者连接都行，其工具检索双语支持，也可以批处理
Python
- 一种是基于ArcPy实现Spatial Join工具的脚本即可，查看tool help的代码即可，支持批处理，略
- 使用geopandas sjoin，但是没法直接一步到位，以下是我的处理，先sjoin在merge，剔除无关字段
我这里需要按月统计所有格网内的微博点数量，列名等参数可以自行更改，count()也可以换成求和或求值等，或者同时计算更多的类型使用pivot table，此处最后的结果是：在原有的shapefile表上增加了12列，每一列对应当月格网内微博点的数量
主要是以下几步：空间连接→计数→连接到polygon表
1
2
3
| table_sjoin = gpd.sjoin(points_table, polygons[['gid', 'geometry']])
table_df = pd.DataFrame(table_sjoin.groupby('gid')['_id'].count())
polygons = pd.merge(left = polygons,right=table_df, on='gid', how='left').fillna(0)
|
完整示例：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
| outdir = r"A:\my_research\dataset\wb1720\2017m"
polygons = grids.copy(deep=True)
for month, table in wb2017.groupby('month'):
outfile = outdir+"\\m"+str(month)+".shp"
table.to_file(outfile, encoding='utf-8',index=False)
table_sjoin = gpd.sjoin(table, polygons[['gid', 'geometry']])
table_df = pd.DataFrame(table_sjoin.groupby('gid')['_id'].count())
table_df = table_df.reset_index()
new_name = "count_m"+str(month)
table_df = table_df.rename(columns={'_id':new_name})
polygons = pd.merge(left = polygons,right=table_df, on='gid', how='left').fillna(0)
polygons[new_name] = polygons[new_name].astype('int32')
print(month)
polygons.to_file("./data/results.shp", encoding = 'utf-8')
|