VISSIM二次开发（Python）&大作业总结2

写在前面

前一小节已经概括了基本所有这次使用的二次开发的知识，在这一小节，我们主要介绍本文对.att文件的读取和绘图以及分析的工作。

在这一部分我们将展开对这一部分的介绍，由于这部分的内容专门性比较强，所以就与前面的二次开发部分分开撰写，供有需要的同学们参考。

.att文件示例

在这里插入图片描述

其中表头会介绍数据字节含义，这里需要点出第一列就是仿真运行次数。

需求1：依据仿真结果绘制contour图

#开发环境
matplotlib                         3.3.4
numpy                              1.19.5
pandas                             1.2.4

coutour图是指等高线图，直接上代码

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

S=8  #设置需要读取的仿真次数

df1=pd.read_csv('file.att',sep=';',skiprows=21)#att中数据分隔符是“;”，然后需要跳过前面不需要转换成df的部分。
#读取数据需要跳过前面的不需要的行
#%%
# df2=df1.sort_values(by=['DATACOLLECTIONMEASUREMENT', 'TIMEINT'])

dfn1=df1[df1['$DATACOLLECTIONMEASUREMENTEVALUATION:SIMRUN'].isin([S])]#这里其实在前面更改clomun后会更加方便。选择仅仅包含这第八次运行的结果
T=int(672*5/8)#这是其中非匝道段的线圈有效数据的行数，因为后面会有很多不需要的统计值，可能需要自己更改
#%%
dfn=dfn1[dfn1['DATACOLLECTIONMEASUREMENT'].isin([1,2,4,7,9])]#这里选取了其中需要分析的截面
pre_df=dfn[['TIMEINT','DATACOLLECTIONMEASUREMENT','SPEED(ALL)']].iloc[:T,:]#这里主要是选择需要分析的几列
pre_df1=pre_df[['SPEED(ALL)']]#最后发现还是只需要一列
pre_df1.index = range(T)#为方便，对index重新编号

#下面就是把上面的一列数据转换
for i in range(84):
    if i ==0:
        A= pre_df1.iloc[0:5,0]
        C=pd.DataFrame(A)
        C.index = range(len(C))
    else:
        B=pd.DataFrame(pre_df1.iloc[i*5:i*5+5,0])
        B.index = range(len(B))
        C=pd.concat([C,B],axis=1,join='inner')#建立一个9*64的df作为等高线绘图的输入

t=C.shape[1]#下面就是等高线图的绘制
C=C.iloc[::-1,:]
x=np.linspace(0,t-1,t)
y=np.linspace(37,41,5)
X,Y=np.meshgrid(x,y)   #上面三句话就是建立等高线图的网格 
plt.figure()
plt.contourf(X,Y,np.array(C),7,cmap=plt.cm.RdYlGn,levels=np.linspace(0,90,50)) #绘制等高线图，中间需要以矩阵的形式，并且满足x列，y行
plt.colorbar()#显示色条
plt.xticks(range(0,84,12),['6:00','7:00','8:00','9:00','10:00','11:00','12:00'])#更改横纵坐标
plt.yticks([37,38,39,40,41],['NHNX41','NHNX40','NHNX39','NHNX38','NHNX37'])
# plt.savefig('C:\\Users\\13968\\Desktop\\0819-'+str(S)+'.jpg', dpi=1000)
plt.show()

其中控制绘图效果的主要就是levels=np.linspace(0,90,50)这个代码其中0，90是控制显示的纵坐标色彩范围（对应颜色），50是表述色彩过渡的趋势，可以参考下面链接中的方法进行选用。以及色板的使用都可以参考，当然官方文档写的也很清楚~

contourf的colorbar如何设置显示范围_solmn的博客-CSDN博客_colorbar设置想要的范围

【Python】绘制热力图seaborn.heatmap，cmap设置颜色的参数_小白兔de窝-CSDN博客_cmap

部分同学可能也有绘制真实场景的contour的需求，因为方法一致这里也不进行探讨了，主要就是构建这个绘图速度时空矩阵C，大家可以自己尝试，这部分中间的df传递有些混乱_{目前看来能达到效果就好}

最终绘制的效果如下~是不是还可以

在这里插入图片描述

需求2：依据结果进行分析

结果分析主要是参考这两篇文献来做的

Bottleneck Identification and Calibration for Corridor Management Planning

Calibration of a micro-traffic simulation model with respect to the spatial-temporal evolution of expressway on-ramp bottlenecks

评价指标	含义	要求
C1(Bottleneck Area Matching)	实际与仿真在瓶颈范围（时间和空间上）的匹配程度	>0.75
C2(Actual Speed Matching)	反映范围匹配以及实际速度的匹配	>0.70
GEH	由Geoffrey E. Havers提出，旨在比较两组流量数据的匹配程度	<5
速度相对误差		<15%

也简单罗列一下公式，其中C1和Devs还有一些其他的形式不过差别不大，具体怎么算可以参考前面推荐的文献~
$\\Devs= \frac { | S _ { r } - S _ { s } j } { S _ { r } }$

$\sqrt { G H } = \sqrt { \frac { ( E - v ) ^ { 2 } } { ( E + v ) / 2 } }$

$\\C_1=\frac{\sum _ { i = 1 } ^ { N } \{ ( \sum _ { t = 1 } ^ { T } [ B S _ { s } ( i , t ) A B S _ { r } ( i , t ) ] ) \cdot ( x _ { i } + 1 - x _ { i } ) \}}{\sum _ { i = 1 } ^ { N } ( \sum _ { t = 1 } ^ { T } [ B S _ { s } ( i , t ) V B S _ { r } ( i , t ) ] ) \cdot ( x _ { i } + 1 - x _ { i } ) \}}$

$\\C_2=\frac{2 \sum _ { i = 1 } ^ { N } \sum _ { i = 1 } ^ { T } \{[ B S _ { s } ( i , t ) V B S _ { r } ( i , t ) ] \cdot| S _ { s } ( i , t ) - S _ { r } ( i , t ) | ) \cdot ( x _ { i } + 1 - x _ { i } ) \}}{\sum _ { i = 1 } ^ { N } \sum _ { t = 1 } ^ { T } \{ [ B S _ { s } ( i , t ) V B S _ { r } ( i , t ) ] \cdot ( S _ { s } ( i , t ) + S _ { r } ( i , t ) ) ) \cdot ( x _ { i } + 1 - x _ { i } ) \}}$

这部分代码为好兄弟xjt今天看论文了嘛的博客_CSDN博客提供

#文件导入部分前面有讲解，在此不再赘述，主要介绍主要其功能性模块
#计算BS矩阵（45代表阈值，可自行调节，索引自行依据矩阵调节）
    BS_wn[f'speed_{data}'] = BS_wn[f'speed_{data}'].apply(lambda x: 0 if x >= 45 else 1)   #这里直接利用pandas读取后的Dataframe类型来进行操作，其他类型也可以操作，只要索引对应上
#修正BS矩阵 怕出错这里可不带入，对结果影响不大
BS_s_1 = BS_wn.values   #转换成矩阵形式，遍历一个线圈的历史数据进行判断
for spa_ii in range(1,BS_s_1.shape[1]):
    for tt_3 in range(BS_s_1.shape[0]):
        if BS_s_1[tt_3][spa_ii] == 0:  #主要关注瓶颈区的设置是否合理，因此不合理的瓶颈区设置的异常值就是0，所以找到所有0（非瓶颈区位置）
            for tt_4 in range(tt_3-3,tt_3):
                data_listwn = []
                if tt_4 <0:
                    pass
                else:
                    for tt_5 in range(tt_4,tt_4+5):
                        if tt_5 == tt_3 or tt_5 >=len(BS_s_1.shape[0])-1: #空间位置索引最大值
                            pass
                        else:
                            data_listwn.append(BS_s_1[tt_5][spa_ii])
        if data_listwn == [1, 1, 1, 1]:
            BS_s_1[tt_3][spa_ii] = 1
#指标各参数含义可详见公式
#计算C1
def Ca_C1(Bs_s, Bs_r):
    C1_up = np.sum(np.logical_and(Bs_s,Bs_r))
    C1_down = np.sum(Bs_s+Bs_r)
    C1 = 2*C1_up/C1_down
    return C1
# 计算C2指标：
def Ca_C2(Bs_s1, Bs_r1, S_s, S_r):
    Bs_or = np.logical_or(Bs_s1,Bs_r1)
    S_abs = np.abs(S_r-S_s)
    S_add =S_r+S_s
    C2_up = np.sum(np.multiply(Bs_or, S_abs))
    C2_down = np.sum(np.multiply(Bs_or,S_add))
    C2 =1 - 2*C2_up/C2_down
    return C2
# 计算GEH
def Ca_GEH(V_ss,V_r):
    Geh_up = np.multiply((V_ss-V_r), (V_ss-V_r))
    Geh_down = (V_ss+V_r)/2
    Geh = np.sqrt(Geh_up/Geh_down)
    return np.percentile(Geh, 0.85) #返回85分位值
#计算速度相对偏差
def Ca_DevS(S_s1, S_r1):
    S_abs = np.abs(S_r1-S_s1)
    idxnonzeros = np.where(S_r1 != 0)
    Devs = S_abs[idxnonzeros]/S_r1[idxnonzeros]
    Devs_re = Devs[np.where(Devs <= np.percentile(Devs, 0.85))]  #剔除15%的较大的值后进行计算，也可以不剔除
    Devs_avr = np.average(Devs_re)
    return Devs_avr