고속도로 DSRC 데이터 분석
0. 목표
RSE 단말기 정보를 활용한 네트워크 자료(shp) 제작
DSRC 원시자료를 가공하여 경로데이터 형성
1. EDA
DSRC 하루치 데이터가 커서 대용량의 DF를 만들면 메모리에러가 뜸
read_csv 시 날짜 시간데이터를 하나의 컬럼으로 합치고, 컬럼별 자료를 보고 불필요한건 날릴 필요가 있음
import pandas as pd
from tqdm import tqdm
df_dsrc_raw = pd.read_csv('dsrc_raw.csv',
encoding='cp949',
parse_dates={'Datetime' : [0,2]},
dtype={'노선번호':str}
)
df_dsrc_raw.drop('Unnamed: 9', axis=1, inplace=True)
df_dsrc_raw.info()
df_dsrc_raw.isnull().sum()
df_dsrc_raw.sort_values(by=['가상OBU_ID','Datetime'],inplace=True)
'DSRC차종구분명'의 countplot을 보니 차종오류가 많아(16%)정도 제외하는 것이 좋다고 생각됨
* 교통량조사 자료를 통해 보정하는것이 괜찮을 것 같음. 전국기준으로
* 한글 표출을 잊지말자..
import matplotlib
matplotlib.rcParams['font.family'] ='Malgun Gothic'
matplotlib.rcParams['axes.unicode_minus'] =False
df_dsrc_raw.drop(['DSRC차종구분명','DSRC차종구분코드'], axis=1, inplace=True)
나머지 변수 요일명, RSE_ID, 가상OBU_ID, 노선번호, 도로이정은 딱히 살펴볼만한것이 없지만
노선번호의 경우 int값이 아닌 string으로 부르는것이 좋음, 4자리 숫자로 맞춰져있음(ex. 0010 경부선)
## 네트워크 파일을 만들려고 보니 RSE, LINK에 문제가 많음
- RSE 좌표가 없는 값이 있음 - > 도로중심선 좌표를 보간법을 사용하여 대체하였음
- RSE 에서의 남해선 노선번호과 도로중심선 좌표의 남해선 노선번호가 상이함 -> 좌표값 보고 수정함
df_station_lat_lon = pd.read_csv('station_lat_lon.csv', encoding='cp949', dtype={'노선번호':str})
df_station_lat_lon.drop(['GRS80X좌표값','GRS80Y좌표값','도로명'],axis=1,inplace=True)
df_station_lat_lon.rename(columns={'노선번호':'line_num',
'이정':'station',
'X좌표값':'lat',
'Y좌표값':'lon'
},
inplace=True
)
df_station_lat_lon.loc[df_station_lat_lon['line_num']=='0102','line_num'] = '0100'
df_rse_null = df_dsrc_rse[df_dsrc_rse['lat'].isnull()]
df_rse_null_1 = df_rse_null[df_rse_null['line_num']=='0100']
df_rse_null_2 = df_rse_null[df_rse_null['line_num']!='0100']
df_station_lat_lon.sort_values(by=['station'],inplace=True)
from scipy.interpolate import interp1d
def fill_na(df_na,line_num,df_sta):
f_linaer_lat = interp1d(x = df_sta[df_sta['line_num']==line_num]['station'],
y = df_sta[df_sta['line_num']==line_num]['lat']
)
f_linaer_lon = interp1d(x = df_sta[df_sta['line_num']==line_num]['station'],
y = df_sta[df_sta['line_num']==line_num]['lon']
)
df_na.loc[:,'lat'] = df_na.apply(lambda x : f_linear_lat(x['station']),axis=1)
df_na.loc[:,'lon'] = df_na.apply(lambda x : f_linear_lon(x['station']),axis=1)
return df_na
df_rse_null_1 = fill_na(df_rse_null_1, '0100', df_station_lat_lon)
df_rse_null_2 = fill_na(df_rse_null_2, '6000', df_station_lat_lon)
df_rse_null_filled = pd.concat([df_rse_null_1,df_rse_null_2])
df_dsrc_rse.loc[df_dsrc_rse['lat'].isnull(),:] = df_rse_null_filled
df_dsrc_rse = df_dsrc_rse.astype({'lat':'float64',
'lon':'float64'
}
)
df_dsrc_rse.info()
del df_rse_null, df_rse_null_1, df_rse_null_2, df_rse_null_filled
2. 네트워크 파일 제작
geopandas와 fiona를 통해 제작할 수 있는 것으로 파악됨
df_dsrc_link = pd.read_csv('dsrc_link.csv', encoding='cp949')
df_dsrc_rse = pd.read_csv('rse.csv', encoding='cp949', dtype={'위도값' : 'float64',
'경도값' : 'float64'
}
)
위도 경도는 소수점에 의해 큰차이가 있을것같아서 dtype을 지정해줬음
df_dsrc_rse.rename(columns={'위치명' : 'locatio',
'도로이정' : 'station',
'RSE위치구분코드' : 'RSE_location_code',
'노선번호' : 'line_num',
'위치단축명' : 'location_short_name',
'도로등급구분코드' : 'road_class_code',
'위도값' : 'lat',
'경도값' : 'lon',
},
inplace=True)
후에 shp파일 제작시 한글이 있으면 인코딩 문제가 발생하기 때문에 컬럼명을 모두 영어로 변경
위도 경도를 살펴보니 잘못된 데이터가 있었음
위도경도가 바뀐걸로 보임
sns.scatterplot(df_dsrc_rse['lat'],df_dsrc_rse['lon'])
df_dsrc_rse.loc[df_dsrc_rse['RSE_ID'] == 2100701100334, 'lat'] = 37.3083
df_dsrc_rse.loc[df_dsrc_rse['RSE_ID'] == 2100701100334, 'lon'] = 127.019
수정 후 정상화 되었음. data.ex.co.kr 에 수정요청해놨음
제거할까도 생각했는데 그러면 링크 생성시 오류가 생길 것 같아서 수정해서 사용하는 방안을 택함
먼저 노드를 생성하는 것 먼저 배웠음
node_gdf = gpd.GeoDataFrame(df_dsrc_rse,
geometry = gpd.points_from_xy(df_dsrc_rse['lon'], df_dsrc_rse['lat']),
crs='EPSG:4326',
)
type(node_gdf)
node_gdf.plot()
node_gdf.to_file(filename='node.shp',driver='ESRI Shapefile')
노드 shp 생성은 생각보다 간단함
자료와 x,y지정하고 좌표계 입력하면됨
여러가지 옵션이 있을것같지만
필요한것만..
표준노드링크위에 내가 만든 노드를 올려보니 잘 맞는듯함
허나, RSE 좌표계는 EPSG:4326 이고, 표준노드링크의 좌표계는 EPSG:5187 임을 유의해야함
링크를 노드의 station 순서대로 줄세우고, 방향 선택하고, 본선노드만 선택해서 튜플의 리스트로 Linestring을 형성하려했으나,
튀는 노드가 발생함.. 다른 방식을 구상해봐야겠음
import geopandas as gpd
from shapely.geometry import LineString
df_dsrc_link.isnull().sum()
df_dsrc_link_drop = df_dsrc_link.dropna()
df_dsrc_link_drop.isnull().sum()
df_S = df_dsrc_link_drop[df_dsrc_link_drop['start_end']=='S']
link_dict ={}
for line_num in df_S['line_num'].unique():
df_temp = df_S[df_S['line_num']==line_num]
df_temp.sort_values(by=['link_order'],inplace=True)
node_list = []
for index, row in df_temp.iterrows():
node_list.append((row.org_lon,row.org_lat))
link_dict[line_num] = node_list
df_dsrc_link_drop.to_csv('sample.csv')
line = LineString(link_dict['0010'])
gdf = gpd.GeoDataFrame(geometry=[line], crs="EPSG:4326")
gdf.to_file("link.shp", driver="ESRI Shapefile")
그래서 표준노드링크를 geopandas로 불러보고 형태를 똑같이 따라해서 만들어봤음
표준노드링크 shp파일을 불러오면 맨마지막에 'geometry' 라는 컬럼에 linestring이 하나씩 담겨있음
이를 모방해서 내가 가진 dsrc 데이터에 geometry 컬럼을 추가하고 linestring을 넣고 shp를 작성하였음
df_dsrc_link_drop.loc[:,'geometry'] = df_dsrc_link_drop.apply(lambda x : LineString([[x['org_lat'],x['org_lon']],[x['dst_lat'],x['dst_lon']]]),axis=1)
df = gpd.GeoDataFrame(df_dsrc_link_drop, geometry='geometry')
df.to_file('MyGeometries.shp', driver='ESRI Shapefile')
첫째로, 기존 노드파일 생성한거와 맞지 않음.. 좌표계가 다른가? 소수점이 다른가?
둘째로, 튀는 링크가 너무 많음.. 왜지?
도로중심선을 기준으로 모두 재생성해버릴까 고민중임
- 경도 위도를 바꿔서 입력해서 이상했던거였음 수정함
표준노드링크랑 비교(붉은색 : 표준노드링크, 파란색 : 만든 노드 및 링크) 해보니 어긋나는 부분이 있어서 도로중심선으로 수정할거임
3. RSE 좌표를 표준노드링크에 매칭
좌표가 튀는 것을 교정하고 싶기도 하고
RSE 좌표 중 표준링크과 가장 가까운 노드에 배정시켜봤다
from math import dist
coords_list = [list(i) for i in zip(gpd_node_ktdb['lat'],gpd_node_ktdb['lon'])]
df_dsrc_rse[['NODE_ID','NODE_lat','NODE_lon']]=0
for i in tqdm(range(len(df_dsrc_rse))):
RSE_coords = [df_dsrc_rse.iloc[i,7],df_dsrc_rse.iloc[i,8]]
distances = [dist(RSE_coords,NODE_coord) for NODE_coord in coords_list]
closest_idx = distances.index(min(distances))
df_dsrc_rse.iloc[i,9] = gpd_node_ktdb['NODE_ID'][closest_idx]
df_dsrc_rse.iloc[i,10] = gpd_node_ktdb['lat'][closest_idx]
df_dsrc_rse.iloc[i,11] = gpd_node_ktdb['lon'][closest_idx]
하지만 실패했다
흠....
RSE좌표를 도로중심선으로 바꿔줘야하나..
노선명을 통일시키긴했는데..
빨간점이 RSE의 좌표고, 파란점이 RSE에서 가장 가까운 노드
그 노드를 링크로 이었다
이걸로는 튀는걸 잡을 수 없음
튀는애가 가까운 표준노드로 갈뿐..
논리가 잘못됬음
중심선 매칭으로 다시 해봐야겠다
4. RSE 좌표를 도로중심선 좌표에 매칭
생각해보니 raw 파일 EDA는 너무 성급함
네트워크 파일 먼저
도로중심선을 보간법으로 생성하였음
이정 범위가 다르면 수정해주고
#%% match station
df_dsrc_rse[['sta_lat','sta_lon']] = None
cond_1 = (df_dsrc_rse['RSE_location_code'] == 1) | (df_dsrc_rse['RSE_location_code'] == 4)
for line_num in df_dsrc_rse['line_num'].unique():
cond_2 = df_dsrc_rse['line_num'] == line_num
cond_3 = df_station_lat_lon['line_num'] == line_num
if df_dsrc_rse[cond_1&cond_2]['station'].max() <= df_station_lat_lon[cond_3]['station'].max():
pass
else:
print('line_num :', line_num)
print('rse range :', df_dsrc_rse[cond_1&cond_2]['station'].min(),'~',df_dsrc_rse[cond_1&cond_2]['station'].max())
print('sta range :', df_station_lat_lon[cond_3]['station'].min(),'~',df_station_lat_lon[cond_3]['station'].max())
print('---')
'''
line_num : 0150
rse range : 8.22 ~ 340.8
sta range : 0.0 ~ 340.7
---
line_num : 1100
rse range : 0.1 ~ 26.6
sta range : 0.0 ~ 26.6
---
line_num : 0140
rse range : 0.91 ~ 41.7
sta range : nan ~ nan
---
line_num : 5000
rse range : 1.92 ~ 9.63
sta range : nan ~ nan
---
line_num : 0141
rse range : 100.1 ~ 145.1
sta range : nan ~ nan
---
'''
def interp_xy(df,line_num,df_sta,cond_1,cond_2,cond_3):
print(line_num)
if len(df_sta[cond_3])==0:
return
elif df[cond_1&cond_2]['station'].max() > df_sta[cond_3]['station'].max():
df.iloc[df[cond_1&cond_2]['station'].idxmax(),2] = df_sta[cond_3]['station'].max()
elif df[cond_1&cond_2]['station'].min() < df_sta[cond_3]['station'].min():
df.iloc[df[cond_1&cond_2]['station'].idxmin(),2] = df_sta[cond_3]['station'].min()
f_linear_lat = interp1d(x = df_sta[cond_3]['station'],
y = df_sta[cond_3]['lat']
)
f_linear_lon = interp1d(x = df_sta[cond_3]['station'],
y = df_sta[cond_3]['lon']
)
df.loc[cond_1&cond_2,'sta_lat'] = df[cond_1&cond_2].apply(lambda x : f_linear_lat(x['station']),axis=1)
df.loc[cond_1&cond_2,'sta_lon'] = df[cond_1&cond_2].apply(lambda x : f_linear_lon(x['station']),axis=1)
cond_1 = (df_dsrc_rse['RSE_location_code'] == 1) | (df_dsrc_rse['RSE_location_code'] == 4)
for line_num in df_dsrc_rse['line_num'].unique():
cond_2 = df_dsrc_rse['line_num'] == line_num
cond_3 = df_station_lat_lon['line_num'] == line_num
interp_xy(df_dsrc_rse, line_num, df_station_lat_lon, cond_1, cond_2, cond_3)
sum(df_dsrc_rse['sta_lat']=='None')
df_dsrc_rse.info()
df_dsrc_rse[['sta_lat','sta_lon']] = df_dsrc_rse.astype({'sta_lat':'float64','sta_lon':'float64'})[['sta_lat','sta_lon']]
df_dsrc_rse.loc[df_dsrc_rse.isnull()['sta_lat'],'sta_lat'] = df_dsrc_rse.loc[df_dsrc_rse.isnull()['sta_lat'],'lat']
df_dsrc_rse.loc[df_dsrc_rse.isnull()['sta_lon'],'sta_lon'] = df_dsrc_rse.loc[df_dsrc_rse.isnull()['sta_lon'],'lon']
del cond_1, cond_2, cond_3, line_num
첫째,
분홍색이 RSE 좌표, 파란색이 표준노드링크 좌표, 빨간색이 도로중심선 좌표인데
도로 그리기를 빨간색 기준으로해서그런지 제일 잘맞음
RSE도 도로변에 있는걸 생각하면 상당히 정확한듯?
표준노드링크는 잘모르겠음
둘째,
여전히 튀는걸 잡을 수 없음
링크데이터에 이상이 있는걸로 생각하고
RAW 경로 분석을 시작해야겠음
5. RAW 다시 살펴보기
어제 발견한 dataprep이라는 패키지를 활용해서 변수들을 살펴보겠음
오후 5시에 최대 204만대의 DSRC 데이터가 담김
차종을 보니 1종과 승용차가 나눠서 찍힘
의미가있나..
차종오류도 15.17%로 큼 1종은 77.52%
총 38개의 노선이 있으며 경부선 21.11%로 데이터가 가장 많음
missing value도 없는걸로
불필요한 차종데이터를 없애고 경로자료를 만들어봐야겠다
5. trip 구분해내기
다른 논문을 보니 같은 OBU_ID를 가진 통행에 대해서 구분을 하기위해
시간차이를 사용했다
시간차이가 중앙값+3*median(abs(Xi-median(X)) 보다 작으면 같은 통행이라는건데
연산을 해보니 좀 불합리한 것 같다.
다른 논문을 살펴보기 3MAD를 사용한건 이상치를 구별하는 방법이라고한다
평균에 비해 median이 로버스트하니 사용한다고 한다.
흠
관측시간, x좌표, y좌표, 차종으로 GMM 클러스터링으로 통행을 구분해야겠다
첫번째 난관은 몇개의 트립으로 군집화해야할지 정할 수 없다는 것
살펴보니 AIC, BIC 스코어로 정한다고 한다.
계산법이나 원리는 모른다
n_components = np.arange(1, 50)
models = [GaussianMixture(n, covariance_type='full', random_state=0).fit(data_np) for n in n_components]
plt.plot(n_components, [m.bic(data_np) for m in models], label='BIC')
plt.plot(n_components, [m.aic(data_np) for m in models], label='AIC')
plt.legend(loc='best')
plt.xlabel('n_components');
데이터(OBU_ID==1)를 살펴보면 17~18개의 트립이 발생할 것 같은데
BIC 점수를 보니 4개의 트립이 적당하다고 계산된다.
feature가 너무 적나?
분산이 너무 작나?
이유가 뭐지??
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
xs = data_np[:,0]
ys = data_np[:,1]
zs = data_np[:,2]
color = df_temp['group']
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(xs, ys, zs, c=color, marker='o', s=15)
ax.set_xlabel('sta', fontsize=20)
ax.set_ylabel('lon', fontsize=20)
ax.set_zlabel('time', fontsize=20)
plt.show()
