基于路网和GeoPandas的高斯两步移动搜索法可达性分析
2024-06-16 02:17:55
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']
供需比计算
根据公式：
Rj=Sj∑k∈{dkj≤d0}G(dkj,d0)DkR_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)DkSj
进行计算
# 首先确定搜索半径
# 本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=0for 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值
在本步中，我们将从居民点出发，寻找搜索半径内所有的应急避难所入口
具体公式如下:
Ai=∑f∈{dij≤d0}RjA_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=0for poi in points:# 依旧存在对应不上的问题try:ans+=R_N[poi[1]]except:ans+=0A_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改进
目前考虑的改进方向为：
- 自适应搜索半径
- 更加合理的泊松车流阻尼
- 更加复杂的人地关系模型
- 数据源的正确性修正
基于路网和GeoPandas的高斯两步移动搜索法可达性分析相关推荐
- 基于高斯两步移动搜寻法（2SFCA）的城市绿地可达性分析
[2SFCA的基本思路,可以略过] 对每个供给点j,搜索所有在j搜寻半径(d0)范围内的需求点(k),计算供需比Rj:对每个需求点i,搜索所有在i搜寻半径(d0)范围内的供给点(j),将所有的供需比R ...
- [续]基于高斯两步移动搜寻法（2SFCA）的城市绿地可达性分析[格网]
原文链接 [续]基于高斯两步移动搜寻法(2SFCA)的城市绿地可达性分析[格网]https://mp.weixin.qq.com/s/BkjZVqCAORrKH5bsyOSTbg 之前我用立方学社公开 ...
- CV：基于keras利用cv2自带两步检测法对《跑男第六季第五期》之如花片段(或调用摄像头)进行实时性别脸部表情检测
CV:基于keras利用cv2自带两步检测法对<跑男第六季第五期>之如花片段(或调用摄像头)进行实时性别&脸部表情检测 目录 输出结果 设计思路 核心代码 输出结果 设计思路 核心 ...
- 两步移动搜索法介绍与原理
两步移动搜索法 两步移动搜寻法是公共服务设施空间可达性研究中的重要方法,在国内外公共服务设施布局研究中得到了广泛应用,本文介绍了两步移动搜索法的基本内容与原理. 一.两步移动搜索法是什么? 基本两步移 ...
- 2020FME博客大赛——基于FME利用高德路径规划AP实现公共服务设施可达性分析——以厦门山海健康步道为例
作者:郭文义 单位:厦门市环境科学研究院 等时圈(siochrone),指从某点出发,以某种交通方式在特定时间内能到达的距离覆盖的范围(来自于网络).(An isochrones is an isol ...
- 岳翔南京大学计算机,基于组合IIS路径抽取的组合线性混成系统有界可达性分析-中国科学.PDF...
基于组合IIS路径抽取的组合线性混成系统有界可达性分析-中国科学 中国科学: 信息科学 2017 年 第47 卷 第3 期: 288–309 论文 基于组合 路径抽取的组合线性混成系统有界 可达性分析 ...
- 两步路轨迹文件位置_#APP功能#两步路轨迹路网，户外导航不迷路！
亲爱的驴友! 你是否曾在登山过程中,因为天气等原因导致前路受阻,绕道后却又迷路? 你是否为了尽快与前队队友集合,一路狂奔埋头往前冲,最终搞得筋疲力尽? 你又是否为了护送受伤的队员就医,必须半途返程,却 ...
- 量子：基于EPR块对的两步量子直接通信
学习论文: 题目:Two-step quantum direct communication protocol using the Einstein-Podolsky-Rosen pair block ...
- java基于聚类的离群点检测_挑子学习笔记：基于两步聚类的离群点检测
转载请标明出处:http://www.cnblogs.com/tiaozistudy/p/anomaly_detection.html 本文主要针对IBM SPSS Modeler 18.0中离群点检 ...
最新文章
- 清华贵系的期末大作业：奋战三周，造台计算机！
- 【ACL2020】这8份Tutorial不可错过！包括：常识推理、多模态信息抽取、对话、解释性等...
- SAP FI新手常用代码
- Spring查找方法示例
- 另一种声音：容器是不是未来？
- 3-11 Matplotlib数据可视化基础
- html5 动态 menuitem,利用HTML 5中的Menu和Menuitem元素快速创建菜单
- UBUNTU修改控制台语言
- 制造-销售”模式正在消亡，传统大型企业的上云之路要如何举步？
- Mac电脑上的Safari运行缓慢，卡的要死，该怎么解决？
- Java IO 系统（一）
- Linux 网络配置 修改DNS配置文件/etc/resolv.conf后，重启网络，DNS配置丢失
- PMP试题 | 每日一练，快速提分 9.1
- python判断信用卡号是否合法_三十八、练习、Python判断一个信用卡号是否合理
- java net php_atitit. web 在线文件管理器最佳实践(1)--- elFinder 的使用流程解决之道 。打开浏览服务器文件夹java .net php...
- MATLAB泰勒级数展开
- 最大化参数 火车头_火车头采集（LocoySpider）设置技巧
- GitLab Projects 2020 插件配置
- vga分屏2个显示器_VGA多分屏聚合器和有趣的问题
- 【读官方文档，学原味技术】SpringBoot-Staters和自定义Starter