基于GeoPandas的POI人口数赋值方法

这个方法是某篇文章中提到的，基于未知兴趣点和街道中心人口点进行的未知兴趣点人口赋值。

我们先来说一下数据，street是街道面数据，里面有一个population字段，用来记录街道总人口值。有一个字段Join_Count表示街道包含的原始点数量。

poi是渔网点数据，由渔网数据点获取邻近500米范围内的原始点(与人口斯皮尔曼系数最高的三类点：交通站点、学校、居民点)叠加而来，有一个属性Join_Count表示在邻接500米缓冲区内的原始点数量。

我们的目的是：

• 统计每一个街道内部的点
• 通过该点与街道人口的关系，计算该点的人口

import geopandas as gpd

path=r"\data"

import numpy as np
import pandas as pd
import math
import random

EDA阶段

• 查看数据
• 统一坐标系

poi=gpd.read_file(path+"/poi.json")

street=gpd.read_file(path+"/street.json")

poi.head()

	OBJECTID	FID_meshmap4_label	FID_meshmap_poi_jc1	Join_Count	TARGET_FID	BUFF_DIST	ORIG_FID	OBJECTID_1	NAME	KIND	WGS_Lon	WGS_Lat	geometry
0	1	249	125	1	216	500	249	216	呈贡县中卫小学	A702	102.776100	24.761585	POINT (102.77741 24.75828)
1	2	287	126	1	231	500	287	231	福宝宝幼儿园	A701	102.781719	24.768870	POINT (102.77741 24.76728)
2	3	365	127	1	265	500	365	265	野山药火腿鸡风干腊排骨美食城马金铺店	1380	102.796260	24.789305	POINT (102.79541 24.78528)
3	4	366	128	10	266	500	366	266	化成路口	808B	102.802564	24.787729	POINT (102.80441 24.78528)
4	5	367	129	1	267	500	367	267	兴盛园	1380	102.817426	24.782900	POINT (102.81341 24.78528)

street.head()

	OBJECTID_1	Join_Count	TARGET_FID	name	name_1	population	Ename	Shape_Leng	QUYU_NAME	ID	OBJECTID	NAME_12	KIND	WGS_Lon	WGS_Lat	Shape_Length	Shape_Area	geometry
0	1	174	0	阿拉街道	阿拉	87545	Ala	56728.703706	官渡	32	1907	羊甫车场	808F	102.793098	24.975722	56728.703396	7.544723e+07	POLYGON ((585551.541 2772498.155, 585514.492 2...
1	2	191	1	茨坝街道	茨坝	86047	Ciba	32123.282283	盘龙	1	3373	蓝龙潭	808B	102.726554	25.127394	32123.282373	3.581656e+07	POLYGON ((575100.675 2789697.506, 575084.789 2...
2	3	295	2	大观街道	大观	48914	Daguan	6201.573449	五华	15	1398	白药厂	808B	102.694308	25.034794	6201.573654	2.080762e+06	POLYGON ((569708.611 2771437.444, 569406.681 2...
3	4	18	3	大渔街道	大渔	23084	Dayu	38768.162785	呈贡	43	3947	关高路口	808B	102.805369	24.814496	38768.163019	2.470458e+07	POLYGON ((578718.634 2749972.017, 578715.300 2...
4	5	305	4	滇池国家旅游	滇池度假区	121943	Dianchi	21422.100871	西山	26	321	海埂	808B	102.673406	24.962374	21422.101000	2.186529e+07	POLYGON ((567755.987 2767210.931, 567731.801 2...

poi.plot()
street.plot()

<AxesSubplot: >

poi=poi.to_crs(epsg=4543)
street=street.to_crs(epsg=4543)1

poi.head(1)

	OBJECTID	FID_meshmap4_label	FID_meshmap_poi_jc1	Join_Count	TARGET_FID	BUFF_DIST	ORIG_FID	OBJECTID_1	NAME	KIND	WGS_Lon	WGS_Lat	geometry	pop
0	1	249	125	1	216	500	249	216	呈贡县中卫小学	A702	102.7761	24.761585	POINT (578633.905 2739501.580)	0

street.head(1)

	OBJECTID_1	Join_Count	TARGET_FID	name	name_1	population	Ename	Shape_Leng	QUYU_NAME	ID	OBJECTID	NAME_12	KIND	WGS_Lon	WGS_Lat	Shape_Length	Shape_Area	geometry
0	1	174	0	阿拉街道	阿拉	87545	Ala	56728.703706	官渡	32	1907	羊甫车场	808F	102.793098	24.975722	56728.703396	7.544723e+07	POLYGON ((585551.541 2772498.155, 585514.492 2...

按照街道选取与街道相交的点

• 遍历街道数据，选择每个与街道数据相交的poi
• 对该poi的人口数据进行更新，依据公式：
  $$PO{I}_p=STREE{T}_p\ast \frac{PO{I}_{counter}}{STREE{T}_{count}}$$

key=street['OBJECTID_1'].unique()

poi["pop"]=[0]*len(poi)

for i in key:
    n=street[street['OBJECTID_1']==i]
    node=n.geometry
    # 这边不知道怎么了，出了点问题，不能用广播，只能这样算了
    for j in range(len(poi)):
        P=poi.iloc[j,:]
        p=P.geometry
        if any(node.contains(p)):
            poi.iloc[j,-1]=int(math.ceil(n.population/n.Join_Count*P.Join_Count))

poi.head()

此时我们可以发现，已经按照人口比例进行赋值了。

	OBJECTID	FID_meshmap4_label	FID_meshmap_poi_jc1	Join_Count	TARGET_FID	BUFF_DIST	ORIG_FID	OBJECTID_1	NAME	KIND	WGS_Lon	WGS_Lat	geometry	pop
0	1	249	125	1	216	500	249	216	呈贡县中卫小学	A702	102.776100	24.761585	POINT (578633.905 2739501.580)	1554.0
1	2	287	126	1	231	500	287	231	福宝宝幼儿园	A701	102.781719	24.768870	POINT (578628.239 2740498.564)	1554.0
2	3	365	127	1	265	500	365	265	野山药火腿鸡风干腊排骨美食城马金铺店	1380	102.796260	24.789305	POINT (580437.259 2742503.011)	1554.0
3	4	366	128	10	266	500	366	266	化成路口	808B	102.802564	24.787729	POINT (581347.440 2742508.339)	15538.0
4	5	367	129	1	267	500	367	267	兴盛园	1380	102.817426	24.782900	POINT (582257.623 2742513.726)	1554.0

poi.to_file(path+"/newPoi.shp")

C:\Users\lenovo\AppData\Local\Temp\ipykernel_27348\2164432190.py:1: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.
  poi.to_file(path+"/newPoi.shp")