基于路网和GeoPandas的高斯两步移动搜索法可达性分析

import warnings
warnings.filterwarnings("ignore")
import geopandas as gpd
import pandas as pd
import os
import collections
import heapq
import re
import math

# 文件路径
path=r"C:\Users\落花雨\Desktop\data"
LIST=["\\"+i for i in os.listdir(path)]

KM_Poi=gpd.read_file(path+LIST[5]) # 昆明兴趣点
YJBNS_Ent=gpd.read_file(path+LIST[2]) # 应急避难所入口数据
YJBNS_Polgon=gpd.read_file(path+LIST[3]) # 应急避难所面数据
ZCQ_Roa=gpd.read_file(path+LIST[1]) # 主城区道路
OD=gpd.read_file(path+LIST[0])# OD 矩阵

数据探查

KM_Poi.head(2)

	FID	OBJECTID	osm_id	name	type	population	QUYU_NAME	geometry
0	0	3	474552927	官渡街道	town	183624	官渡	POINT (576290.931 2761828.261)
1	1	12	1483174846	海口街道	town	69019	西山	POINT (561304.355 2742770.016)

YJBNS_Ent.head(2)

	FID	OBJECTID	NAME	QUYU_NAME	geometry
0	0	6	云南大学呈贡校区-西门	呈贡区	POINT (585001.518 2747248.457)
1	1	7	云南大学呈贡校区-北门	呈贡区	POINT (585528.681 2748211.540)

YJBNS_Polgon.head(2)

	FID	NAME	Shape_Leng	Shape_Area	geometry
0	0	昆明呈贡新区第一小学广电苑校区	610.369700	19710.134256	POLYGON ((11444231.072 2852271.929, 11444275.0...
1	1	四川师范大学附属昆明实验学校（天	966.366133	52802.239035	POLYGON ((11445552.439 2854525.840, 11445615.3...

ZCQ_Roa.head(2)

	FID	FID_roads	osm_id	name	ref	type	oneway	bridge	maxspeed	FID_昆明	...	name_1	center	centroid	childrenNu	level	parent	subFeature	acroutes	length	geometry
0	0	0	23412703.0	Kunming-Shilin Expressway		motorway	1	1	0	2	...	官渡区	0	0	0	district	0	2	0	1152.168205	MULTILINESTRING ((102.73576 25.01994, 102.7363...
1	1	1	24044953.0	白塔路		secondary	1	0	0	1	...	盘龙区	0	0	0	district	0	1	0	1994.046225	MULTILINESTRING ((102.72596 25.05293, 102.7260...

2 rows × 21 columns

数据预处理

# 查看坐标系
KM_Poi.crs # 兴趣点是CGCS2000 / 3-degree Gauss-Kruger CM 102E
# 其坐标系编码为4543

<Projected CRS: EPSG:4543>
Name: CGCS2000 / 3-degree Gauss-Kruger CM 102E
Axis Info [cartesian]:
- X[north]: Northing (metre)
- Y[east]: Easting (metre)
Area of Use:
- name: China - between 100°30'E and 103°30'E.
- bounds: (100.5, 21.13, 103.5, 42.69)
Coordinate Operation:
- name: 3-degree Gauss-Kruger CM 102E
- method: Transverse Mercator
Datum: China 2000
- Ellipsoid: CGCS2000
- Prime Meridian: Greenwich

ZCQ_Roa.crs
# Road的坐标系需要转化为投影坐标

<Geographic 2D CRS: GEOGCS["CGCS_2000",DATUM["D_2000",SPHEROID["S_2000 ...>
Name: CGCS_2000
Axis Info [ellipsoidal]:
- lon[east]: Longitude (Degree)
- lat[north]: Latitude (Degree)
Area of Use:
- undefined
Datum: D_2000
- Ellipsoid: S_2000
- Prime Meridian: Greenwich

YJBNS_Polgon.crs  # 应急避难所场所倒是WGS84
# 需要坐标转换

<Projected CRS: EPSG:3857>
Name: WGS 84 / Pseudo-Mercator
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: World between 85.06°S and 85.06°N.
- bounds: (-180.0, -85.06, 180.0, 85.06)
Coordinate Operation:
- name: Popular Visualisation Pseudo-Mercator
- method: Popular Visualisation Pseudo Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

OD.crs

<Projected CRS: EPSG:4543>
Name: CGCS2000 / 3-degree Gauss-Kruger CM 102E
Axis Info [cartesian]:
- X[north]: Northing (metre)
- Y[east]: Easting (metre)
Area of Use:
- name: China - between 100°30'E and 103°30'E.
- bounds: (100.5, 21.13, 103.5, 42.69)
Coordinate Operation:
- name: 3-degree Gauss-Kruger CM 102E
- method: Transverse Mercator
Datum: China 2000
- Ellipsoid: CGCS2000
- Prime Meridian: Greenwich

YJBNS_Ent.crs

<Projected CRS: EPSG:4543>
Name: CGCS2000 / 3-degree Gauss-Kruger CM 102E
Axis Info [cartesian]:
- X[north]: Northing (metre)
- Y[east]: Easting (metre)
Area of Use:
- name: China - between 100°30'E and 103°30'E.
- bounds: (100.5, 21.13, 103.5, 42.69)
Coordinate Operation:
- name: 3-degree Gauss-Kruger CM 102E
- method: Transverse Mercator
Datum: China 2000
- Ellipsoid: CGCS2000
- Prime Meridian: Greenwich

坐标系转换

YJBNS_Polgon=YJBNS_Polgon.to_crs(epsg=4543)
ZCQ_Roa=ZCQ_Roa.to_crs(epsg=4543)
YJBNS_Ent=YJBNS_Ent.to_crs(epsg=4543)

此时已经完成了统一坐标系的工作

第一步计算

本步骤需要找出应急避难所搜索半径内的所有POI，并计算供需比

构建缓冲区

缓冲区用于搜索供给点

YJBNS_Polgon_C=YJBNS_Polgon.copy()

# 构建缓冲区示例
YJBNS_Polgon_C["geometry"][0].buffer(50)

# 我们不需要中间
YJBNS_Polgon_C["geometry"][0].buffer(50).difference(YJBNS_Polgon_C["geometry"][0])

# 查看应急避难所分布情况
YJBNS_Polgon_C['geometry'].plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c7f17706d8>

# 设置缓冲半径
BUFFER=3000

YJBNS_Polgon_C['geometry']=YJBNS_Polgon_C["geometry"].buffer(BUFFER)

YJBNS_Polgon_C['geometry'].plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c7f178ab00>

曼哈顿距离

相较于欧式距离度量，曼哈顿距离能更反应城市POI之间的实际连通性能

OD.head(2)

	FID	ObjectID	Name	OriginID	Destinatio	Destinat_1	Total_长	geometry
0	0	12851	官渡街道 - 云师大附属世纪金源学校-东门	103	422	1	1296.528731	LINESTRING (576290.931 2761828.261, 575577.741...
1	1	12852	官渡街道 - 云大附中星耀校区-北门	103	420	2	1490.718704	LINESTRING (576290.931 2761828.261, 577144.590...

对于OD成本矩阵,我们希望对每个应急避难所都有一份居民点

在构建的OD矩阵中，Name属性以 街道 - 应急避难所入口 做划分

# 先获取名称列表
Ent_name_list=YJBNS_Ent["NAME"].unique()
Str_name_list=KM_Poi["name"].unique()

S2E=collections.defaultdict(list)# 街道对应入口字典
E2S=collections.defaultdict(list) # 入口对应街道字典

for i in range(OD.shape[0]):
    # 获取并分割名称
    Name=list(map(lambda x:x.strip(),OD.iloc[i,2].split("-",1)))
    # 记录对应的值
    dis=OD.iloc[i,6]
    # 街道
    # 发现数据会被重复读写三次！！
    # 不知道为什么！！
    heapq.heappush(S2E[Name[0]],(dis,Name[1]))
    heapq.heappush(E2S[Name[1]],(dis,Name[0]))

# 此时发现所有记录都已经记载了
for i in Ent_name_list:
    if i not in E2S.keys():
        print(E2S[i])
        
# 马街街道是个神仙
# 因为找不到目的地而被舍弃了
for j in Str_name_list:
    if j not in S2E.keys():
        print(j)

马街街道

# 小根堆会将元素剔除
# 若是要保留元素，可以留下副本
S2E_c=S2E.copy()
E2S_c=E2S.copy()

到这一步我们已经处理了曼哈顿距离，现在能够从避难所入口出发，寻找给定阈值的居民点了

也能够从居民点出发，寻找给定阈值的服务区

下一步我们将获取避难所的面积信息

应急避难所面积

在本步骤中，我们需要获取应急避难所的面积信息，并将其与入口点对应起来

我们先查看避难所面数据

YJBNS_Polgon_C1=YJBNS_Polgon.copy()
YJBNS_Polgon_C1.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1c7f17a6128>

YJBNS_Polgon_C1.head()

	FID	NAME	Shape_Leng	Shape_Area	geometry
0	0	昆明呈贡新区第一小学广电苑校区	610.369700	19710.134256	POLYGON ((581419.907 2745109.067, 581460.202 2...
1	1	四川师范大学附属昆明实验学校（天	966.366133	52802.239035	POLYGON ((582607.927 2747151.992, 582666.062 2...
2	2	云南中医药大学应急避难场所	3778.904004	607319.343216	POLYGON ((582979.894 2749463.449, 582932.370 2...
3	3	捞鱼河湿地公园	3832.693417	513661.330361	POLYGON ((577548.739 2746297.639, 577654.965 2...
4	4	云师大附中	1255.955110	107052.706438	POLYGON ((585882.850 2749957.852, 585741.543 2...

经过数据探查，发现该项数据中的 NAME 属性列能够跟入口名称对应起来，并且 Shape_Area 为几何计算得到的面积

但是两张表之间的名称并非完全匹配。我们只需要匹配关键字即可

# 需要做映射的数组
E_N=E2S.keys()
P_N=YJBNS_Polgon_C1["NAME"]

Name_Area=[]
for i in P_N:
    for j in E_N:
        if re.search(i,j):
            Name_Area.append([j,(YJBNS_Polgon_C1[YJBNS_Polgon_C1["NAME"]==i]["Shape_Area"]).tolist()[0]])

# 创建一张查询表
Area=dict(Name_Area)

现在需要设计一个函数，用于取出给定阈值内的数据

def getData(dataset,threshold):
    tem_dataset=dataset.copy()
    if len(tem_dataset)==0:
        # 空表
        return []
    T=heapq.heappop(tem_dataset)
    ans=[]
    while T[0]<=threshold:
        ans.append(T)
        T=heapq.heappop(tem_dataset)
    return ans 

人口数据

人口数据存放在昆明POI的Population字段中

KM_Poi.head()

	FID	OBJECTID	osm_id	name	type	population	QUYU_NAME	geometry
0	0	3	474552927	官渡街道	town	183624	官渡	POINT (576290.931 2761828.261)
1	1	12	1483174846	海口街道	town	69019	西山	POINT (561304.355 2742770.016)
2	2	14	-1995559001	马金铺街道	town	54798	呈贡	POINT (581339.943 2742863.779)
3	3	15	-1659466282	阿拉街道	town	87545	官渡	POINT (579586.509 2766116.964)
4	5	22	-782694876	双龙街道	town	11841	盘龙	POINT (585541.983 2778974.166)

若要查看其人口数据，只需要以下语句即可

KM_Poi[KM_Poi["name"]=="官渡街道"]["population"].tolist()[0]

183624

# 封装成函数
def getPop(name):
    try:
        return KM_Poi[KM_Poi["name"]==name]["population"].tolist()[0]
    except:
        return -1
    
# 但是这样时间复杂度太高了
# 所以我决定打表
Name_Pop={}
for i in range(KM_Poi.shape[0]):
    T=KM_Poi.iloc[i]
    Name_Pop[T['name']]=T['population']

供需比计算

根据公式：
R j = S j ∑ k ∈ { d k j ≤ d 0 } G ( d k j , d 0 ) D k R_j =\frac{S_j}{\sum_{k\in \{{d_{kj}\le d_0}\}}G(d_{kj},d_0)D_k} Rj​=∑k∈{dkj​≤d0​}​G(dkj​,d0​)Dk​Sj​​
进行计算

# 首先确定搜索半径
# 本demo采用给定阈值的搜索半径
Threshold=3000
# test
getData(S2E_c["官渡街道"],Threshold)

[(1296.52873128, '云师大附属世纪金源学校-东门'),
 (1490.71870426, '云大附中星耀校区-北门'),
 (1600.79115752, '云大附中星耀校区-南门'),
 (1708.15483322, '云大附中星耀校区-1号门'),
 (1782.80973124, '云师大附属世纪金源学校-西南门'),
 (1850.44174205, '云师大附属世纪金源学校-西门'),
 (1955.38161, '云师大附属世纪金源学校-西北门'),
 (2780.42288513, '新亚洲体育城体育中心-2号门'),
 (2960.49871649, '新亚洲体育城体育中心-1号门')]

# test
getData(E2S_c["昆明呈贡新区第一小学广电苑校区"],12000)

[(3219.24822336, '马金铺街道'),
 (3354.50862843, '大渔街道'),
 (5100.55004948, '雨花街道'),
 (8610.87021572, '乌龙街道'),
 (9965.30536041, '龙城街道'),
 (10175.1377117, '洛龙街道')]

我们现在需要处理的就是将这些数据按照公式累加，先确定高斯方法：

def Gauss(dis,dis_D):
    f1=((math.e)**(-0.5))*((dis/dis_D)**2)-math.e*(-0.5)
    f2=1-math.e**(-0.5)
    return f1/f2

# 对每个供给点进行计算
# 保留个副本
YJBNS_Ent_c=YJBNS_Ent.copy()
R_list=[] # 用来存数据的

for i in range(len(Ent_name_list)):
    # 找到在搜索半径内的每个需求点
    name=YJBNS_Ent_c.iloc[i,2]
    points=getData(E2S[name],Threshold)
    # 当然，有些根本就没有
    # 这部分我们会在之后的算法进行改进
    # 我们需要获取这些点的人口数
    ans=0
    for poi in points:
        Name=poi[1]
        dis=poi[0]
        population=Name_Pop[Name]
        ans+=population*Gauss(dis,Threshold)
    # 对于找不到的就需要处理
    # 实际上我们的数据也不是完全一一对应的
    try:
        R_list.append(Area[name]/ans)
    # 目前对于除零(该地区搜索半径没有点)
    # 以及面积不对应(缺失避难所面积)
    # 统一归零处理
    except:
        R_list.append(0)

YJBNS_Ent_c['R_val']=R_list
YJBNS_Ent_c.head()

	FID	OBJECTID	NAME	QUYU_NAME	geometry	R_val
0	0	6	云南大学呈贡校区-西门	呈贡区	POINT (585001.518 2747248.457)	0.0
1	1	7	云南大学呈贡校区-北门	呈贡区	POINT (585528.681 2748211.540)	0.0
2	2	8	云南大学呈贡校区-东门	呈贡区	POINT (586966.101 2747785.881)	0.0
3	3	9	云南大学呈贡校区-南门	呈贡区	POINT (586146.480 2746968.095)	0.0
4	4	10	西南联大研究院附属学校-3号门	呈贡区	POINT (587599.978 2748137.044)	0.0

第二步计算

我们在上一步中，已经获取了R值

在本步中，我们将从居民点出发，寻找搜索半径内所有的应急避难所入口

具体公式如下:

A i = ∑ f ∈ { d i j ≤ d 0 } R j A_i=\sum_{f\in \{{d_{ij}\le d_0}\}}R_j Ai​=f∈{dij​≤d0​}∑​Rj​

KM_Poi_c=KM_Poi.copy()

# 出于线程优化考虑，我们依旧打表
R_N=dict(zip(YJBNS_Ent_c['NAME'],YJBNS_Ent_c['R_val']))

A_list=[]
for i in range(KM_Poi_c.shape[0]):
    name=KM_Poi_c.iloc[i,3]
    points=getData(S2E[name],Threshold)
    # 此时仅需要计算R_val的和
    ans=0
    for poi  in points:
        # 依旧存在对应不上的问题
        try:
            ans+=R_N[poi[1]]
        except:
            ans+=0
    A_list.append(ans)

KM_Poi_c['A_val']=A_list
KM_Poi_c.head()

	FID	OBJECTID	osm_id	name	type	population	QUYU_NAME	geometry	A_val
0	0	3	474552927	官渡街道	town	183624	官渡	POINT (576290.931 2761828.261)	0.245758
1	1	12	1483174846	海口街道	town	69019	西山	POINT (561304.355 2742770.016)	0.000000
2	2	14	-1995559001	马金铺街道	town	54798	呈贡	POINT (581339.943 2742863.779)	0.000000
3	3	15	-1659466282	阿拉街道	town	87545	官渡	POINT (579586.509 2766116.964)	0.000000
4	5	22	-782694876	双龙街道	town	11841	盘龙	POINT (585541.983 2778974.166)	0.000000

存储与显示

KM_Poi_c.to_file(r"C:\Users\落花雨\Documents\Tencent Files\651421775\FileRecv\新建文件夹\Data\Accessibility.shp")

KM_Poi_c.plot()

结果显示如下：

Demo改进

目前考虑的改进方向为：

• 自适应搜索半径
• 更加合理的泊松车流阻尼
• 更加复杂的人地关系模型
• 数据源的正确性修正