서울의 도로교통 패턴 #3

이제부터는 링크별 매시간 평균속도 데이터와 놀 차례다. 원시데이터의 신뢰성을 전제한다면 이것은 통계학적으로 놀 거리가 굉장히 풍부한 panel data다. 다만 나에게 어울려 놀 재주가 부족할 따름이다. 숫자가 내 눈앞에 보기 좋게 펼쳐지지 않으면 찌르고 들어갈 각이 보이지 않기에, 일단은 성실하게 펼쳐 볼 수밖에 없다.

panel data가 대부분 그렇듯 이 데이터도 이 빠진 데가 있다. 가령 이번 작업에 사용할 2017년 10월의 데이터는 4797개 전체 링크 중 스무 개가 빠진 4777개를 포함하고 있다. 매시 평균속도가 빠짐없이 계산되어 4777×24×31 = 3,568,968개의 값이 있으면 좋겠지만 실제로는 28,314개가 결측으로 남아 3,525,774개의 유효값이 있다. 원칙적으로 해당 시간 프로브가 없으면 이력 자료를 활용한다고 했지만, 실제 운영 과정에서는 명시해놓지 않은 기준이 있는 것 같다.

시간대별로는 오전 4-5시 결측값이 6204개이고 나머지 시간대에는 1000개 내외다. 오전 4-5시 결측값이 유난히 많은 이유는 택시의 주야간 교대시간이 4시이기 때문일 것이다. 오후 4시경에도 교대시간의 흔적을 볼 수 있다.

2017년 10월 시간대별 결측값 분포

QGIS로 넘어가기 전, 시간대별 속도의 전체적인 분포를 살펴보기 위해 데이터를 정규확률도 plot에 넣어보았다. 산포도 작성에는 matplotlib을 이용하였다. 시험 삼아 여러 날짜의 산포도를 뽑아봤지만 패턴은 아래와 같이 대동소이하게 나타났다.

2017년 10월 1일(일) / 16일(월) / 27일(금) 시간대별 평균속도의 정규확률도

# -*- coding:utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

def scatterplot(date):
	month = 201710
	csvfile = f"topis.linkid.avg_{month}.csv"
	df24 = pd.read_csv(csvfile)
    fig, ax = plt.subplots()
    for i in range(12):
        hr = f"H{i*2+1:02}"
        sect = df24.loc[lambda df: df.DATE == month*100 + date, hr].sort_values().reset_index()[hr]
        randnm = pd.Series(np.sort(np.random.randn(len(sect))))
        df = pd.DataFrame({'x' : sect, 'y' : randnm})
        ci = (.35+i*.03, .45+i*.02, .45+i*.02)
        ax = df.plot.scatter(x='x', y='y', marker='+', c=ci, alpha=.6, label=hr, ax=ax)
    ax.axis([0,120,-3.5,3.5])
    ax.plot([0,1],[0,1], transform=ax.transAxes, ls=':', c='#aaaaaa', linewidth=1)
    ax.set_title(f'Normal probability plot: {date:02} Oct 2017')
    ax.set_xlabel('Observed (km/h)')
    ax.set_ylabel(r'Expected ($\sigma$)')
    ax.legend(loc=4, edgecolor='w', markerscale=2)
    fig.tight_layout()
    plt.show()

matplotlib을 동원해 plot을 만들기까지 들어간 시간에 비해 그림은 심심한 편이지만, 그래도 몇 가지 기본적인 사실은 확인할 수 있었다.

전체적인 속도 분포는 휴일과 평일 사이에 뚜렷한 차이가 없어 보이는 가운데, 40km/h 이하의 링크가 전체의 80% 이상을 차지하고 있다. 특히 정체가 심한 시간대에는 10km/h 미만 링크의 비율이 60km/h 이상 링크 비율과 비슷하다. 앞에서 서울시 기준에 따라 0 - 15km/h 구간을 하나로 묶었지만 더 세분화할 여지가 있어 보인다. 

반대로 변곡점인 40km/h부터 80km/h까지 가파르게 꺾이는 곡선 패턴은 여러 방향의 해석을 가능하게 한다. 여기에는 제한속도 60km/h만으로 설명하기 힘든 여러 요인이 있을 것이다. 아무튼 이 구간은 상대적 비중이 10% 정도로 매우 낮아, 앞에서 설정했던 40-60km/h 구간이 적절한 것인지 다시 생각하게 한다.

그래프 위쪽 바람에 날리는 스카프 모양은 80km/h 이상의 속도를 보이는 링크로 이루어져 있다. 원칙상 전부 도시고속도로 링크에 해당해야 하는데, 실제 그런지는 지도에서 확인할 수 있을 것이다.

새벽 시간대와 출퇴근 시간대의 속도는 대략 10km/h 내외 차이가 나는데, 잔잔한 물에 떠내려가듯 패턴이 유사하다. 새벽 시간대와 출퇴근 시간대의 교통상황이 전혀 다름에도 전체적인 패턴이 비슷하다는 것은 흥미롭다. 그렇지만 아직 이 그림만으로는 뭔가 말하기 어렵다. 다음 단계로 넘어가 봐야 한다.

누적된 속도 데이터를 지도에 옮기는 두 번째로 단순무식한 방법은, 각 시간대별 지도를 만든 다음 각각의 지도를 1프레임으로 하는 영상을 만드는 것이다. (제일 무식한 방법은 링크의 모든 속도를 ‘산술평균’하여 한 장의 지도로 만드는 것일 테다) 24시간 × 31일 = 744장의 이미지를 가령 12fps 영상으로 만든다면 10월 한 달의 교통상황이 62초로 압축되어 그럴싸한 time-lapse 영상이 될 것이다. 처음에는 이것을 만들 생각이었으나, 몇 가지 기술적 장벽에 부딪혀* 그 작업은 당장은 힘들게 되었다.

각 시간대별로 접근하는 방식을 피한다면, 펼쳐진 데이터를 의미 있게 압축하는 방법이 필요하다. 그렇다면 앞서 그래프를 본 이상 가장 먼저 해보아야 하는 압축은 날짜를 축으로 하여 같은 시간대 속도를 압축하는 것이다. 다시 말해 하루하루 각 시간대별 속도를 조화평균**하여 24시간의 교통 패턴이 어떻게 변화하는지를 보는 것이다. 그리고 휴일 여부를 일종의 dummy 변수로 삼아 그 차이가 있는지 아울러 보았다. 비록 산포도에서는 뚜렷이 나타나지 않았지만 평일과 휴일의 교통 패턴이 다르다는 점은 주지의 사실이다. 마침 지난 10월은 추석연휴 등으로 주말 포함 휴일이 16일에 달해 평일과 양분하는 양상을 보였다.

그래서 2017년 10월의 31일간 시간대별 평균속도를 다시 휴일과 평일 두 그룹으로 나누어 평균하는 과정을 수행하였다. 이 코드를 짤 때까지만 해도 pandas를 제대로 알지 못하고 dict에 의존해 해결하느라 먹지 않아도 될 애를 먹었다.

# -*- coding: utf-8 -*-
import csv
import numpy as np
import numpy.ma as ma
import scipy.stats as stats
import sys
import io
sys.stdout = io.TextIOWrapper(sys.stdout.detach(), encoding = 'utf-8')

def dailyHmean():
    month = 201710
    holidays = [1,2,3,4,5,6,7,8,9,14,15,21,22,28,29]
    mHolidays = [ month*100 + i for i in holidays ]
    with open(f'topis.linkid.avg_{month}.csv', 'r', newline='', encoding='utf-8') as csvfile:
        rawArr = [row for row in csv.reader(csvfile)]
    colnames = rawArr.pop(0)
    linkidDict = {}
    for row in rawArr:
        spds = list(map(lambda x: float(x) if x else np.nan, row[8:32]))
        rid = f'H{row[3]}' if int(row[0]) in mHolidays else f'W{row[3]}'
        if rid not in linkidDict:
            linkidDict[rid] = [np.array(spds)]
        else:
            linkidDict[rid] = np.concatenate((linkidDict[rid], [np.array(spds)]), axis=0)
    with open(f'topis{month}_holidays.csv', 'w', newline='', encoding='utf-8') as csvoutH, \
    open(f'topis{month}_weekdays.csv', 'w', newline='', encoding='utf-8') as csvoutW:
        csvwriteH = csv.writer(csvoutH, delimiter=',')
        csvwriteH.writerow(['TOPIS_ID'] + [f'H{val+1:02}' for val in range(24)])
        csvwriteW = csv.writer(csvoutW, delimiter=',')
        csvwriteW.writerow(['TOPIS_ID'] + [f'H{val+1:02}' for val in range(24)])
        for k, spdsMaskedArr in linkidDict.items():
            spdsMean = np.around(stats.hmean(ma.masked_invalid(spdsMaskedArr), axis=0), 2).tolist()
            topisId = int(k[1:])
            if k[0] == 'H':
                csvwriteH.writerow([topisId] + spdsMean)
            else:
                csvwriteW.writerow([topisId] + spdsMean)
    print(len(linkidDict), 'lines into csv done.')

* 744장의 지도를 일일이 뽑아 저장하는 중노동을 피하려면 ArcGIS의 ArcPy와 같이 QGIS의 scripting language인 PyQGIS를 이용해야 하는데, 이것은 여러 dependency를 포괄하는 거대한 라이브러리여서 단시간에 소화하기 어려웠다. 게다가 QGIS 2.x는 Python 2를 기반으로 하고 있어, Python 3만을 경험한 나는 코딩 과정에서 발생할 여러 버그 - 특히 소위 인코딩 문제 - 에 쉽게 대응할 엄두가 나지 않는다. PyQGIS는 12월 중 Python 3를 기반으로 하는 QGIS 3가 배포되면 다시 들여다볼 생각이다.

** 이하 모든 평균은 조화평균을 말한다. 사실 원시데이터 가공 단계에서 이미 엄밀성을 잃어버린 것으로 보이는 속도 데이터를 산술평균 대신 조화평균으로 처리함으로써 얻는 실익이 얼마나 될지는 의문이지만.