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​)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=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的高斯两步移动搜索法可达性分析相关推荐

  1. 基于高斯两步移动搜寻法(2SFCA)的城市绿地可达性分析

    [2SFCA的基本思路,可以略过] 对每个供给点j,搜索所有在j搜寻半径(d0)范围内的需求点(k),计算供需比Rj:对每个需求点i,搜索所有在i搜寻半径(d0)范围内的供给点(j),将所有的供需比R ...

  2. [续]基于高斯两步移动搜寻法(2SFCA)的城市绿地可达性分析[格网]

    原文链接 [续]基于高斯两步移动搜寻法(2SFCA)的城市绿地可达性分析[格网]https://mp.weixin.qq.com/s/BkjZVqCAORrKH5bsyOSTbg 之前我用立方学社公开 ...

  3. CV:基于keras利用cv2自带两步检测法对《跑男第六季第五期》之如花片段(或调用摄像头)进行实时性别脸部表情检测

    CV:基于keras利用cv2自带两步检测法对<跑男第六季第五期>之如花片段(或调用摄像头)进行实时性别&脸部表情检测 目录 输出结果 设计思路 核心代码 输出结果 设计思路 核心 ...

  4. 两步移动搜索法介绍与原理

    两步移动搜索法 两步移动搜寻法是公共服务设施空间可达性研究中的重要方法,在国内外公共服务设施布局研究中得到了广泛应用,本文介绍了两步移动搜索法的基本内容与原理. 一.两步移动搜索法是什么? 基本两步移 ...

  5. 2020FME博客大赛——基于FME利用高德路径规划AP实现公共服务设施可达性分析——以厦门山海健康步道为例

    作者:郭文义 单位:厦门市环境科学研究院 等时圈(siochrone),指从某点出发,以某种交通方式在特定时间内能到达的距离覆盖的范围(来自于网络).(An isochrones is an isol ...

  6. 岳翔南京大学计算机,基于组合IIS路径抽取的组合线性混成系统有界可达性分析-中国科学.PDF...

    基于组合IIS路径抽取的组合线性混成系统有界可达性分析-中国科学 中国科学: 信息科学 2017 年 第47 卷 第3 期: 288–309 论文 基于组合 路径抽取的组合线性混成系统有界 可达性分析 ...

  7. 两步路轨迹文件位置_#APP功能#两步路轨迹路网,户外导航不迷路!

    亲爱的驴友! 你是否曾在登山过程中,因为天气等原因导致前路受阻,绕道后却又迷路? 你是否为了尽快与前队队友集合,一路狂奔埋头往前冲,最终搞得筋疲力尽? 你又是否为了护送受伤的队员就医,必须半途返程,却 ...

  8. 量子:基于EPR块对的两步量子直接通信

    学习论文: 题目:Two-step quantum direct communication protocol using the Einstein-Podolsky-Rosen pair block ...

  9. java基于聚类的离群点检测_挑子学习笔记:基于两步聚类的离群点检测

    转载请标明出处:http://www.cnblogs.com/tiaozistudy/p/anomaly_detection.html 本文主要针对IBM SPSS Modeler 18.0中离群点检 ...

最新文章

  1. 清华贵系的期末大作业:奋战三周,造台计算机!
  2. 【ACL2020】这8份Tutorial不可错过!包括:常识推理、多模态信息抽取、对话、解释性等...
  3. SAP FI新手常用代码
  4. Spring查找方法示例
  5. 另一种声音:容器是不是未来?
  6. 3-11 Matplotlib数据可视化基础
  7. html5 动态 menuitem,利用HTML 5中的Menu和Menuitem元素快速创建菜单
  8. UBUNTU修改控制台语言
  9. 制造-销售”模式正在消亡,传统大型企业的上云之路要如何举步?
  10. Mac电脑上的Safari运行缓慢,卡的要死,该怎么解决?
  11. Java IO 系统(一)
  12. Linux 网络配置 修改DNS配置文件/etc/resolv.conf后,重启网络,DNS配置丢失
  13. PMP试题 | 每日一练,快速提分 9.1
  14. python判断信用卡号是否合法_三十八、练习、Python判断一个信用卡号是否合理
  15. java net php_atitit. web 在线文件管理器最佳实践(1)--- elFinder 的使用流程解决之道 。打开浏览服务器文件夹java .net php...
  16. MATLAB泰勒级数展开
  17. 最大化参数 火车头_火车头采集(LocoySpider)设置技巧
  18. GitLab Projects 2020 插件配置
  19. vga分屏2个显示器_VGA多分屏聚合器和有趣的问题
  20. 【读官方文档,学原味技术】SpringBoot-Staters和自定义Starter

热门文章

  1. 歪解还是正解的一个字....
  2. How to GROUD?
  3. 中文计数法亿兆京垓秭穰沟涧正载
  4. 王坚:我为什么反对有些企业的“去IOE”运动?
  5. 互动百科疑借“反垄断”进行炒作
  6. Android禁止EditText输入特殊字符
  7. Python课实例5:身体质量指数BMI
  8. 163电子邮箱能免费注册吗?163电子邮件注册移动办公解决方案
  9. 可口可乐公司推出全球第一款可加热饮用汽水,“可口可乐生姜+”上市
  10. 电脑被绑架开机自动装流氓软件