R) 전처리 - 지점간 거리 계산

R) 전처리 - 지점간 거리 계산

실제 지도상의 위치간 거리는 유클리드 거리로 계산하면 안된다. 그 이유와 계산방법을 서울 2호선 지하철역 데이터와 함께 알아보자.

개요

두 지점간의 거리 계산에서 왜 유클리드 거리 계산법을 쓰면 안되는 것일까? 정확하게는 지표면의 두 점의 거리계산에서 유클리드 거리 계산법을 사용하면 안된다. 아무튼 일상에서 사용하는 데카르트의 직교좌표계 위의 두 점의 최단거리 계산의 경우 유클리드 거리를 계산한다. 하지만 지표면의 경우 평면이 아니라 구면상의 두 점의 최단거리를 계산하는 문제이기 때문에 무작정 유클리드 거리를 계산하면 실제 최단거리 보다 더 짧은 값을 산출하게 된다.
거리 계산

그래서 지리좌표계를 고려하여 계산을 해야하며 R을 활용한 연산을 해보도록 하자.


준비

여기에서 사용하는 데이터는 서울 열린데이터 광장 에서 제공하는 서울특별시 노선별 지하철역 정보(신규)를 추가 가공한 데이터를 사용한다.
서울특별시 노선별 지하철역 정보
seoul_subway_line_2_geocoded.csv 다운받기 [클릭]

거리 계산을 정확하게 하기 위해서는 geosphere 패키지를 사용하는 것이 좋다. 미리 설치하고 불러오도록 하자.

1
library("geosphere")

기본 연산

geosphere 함수에는 지표면상 거리 계산 관련 함수가 많지만 그 중에서 두 지점 사이의 최단거리를 계산하는 distGeo() 함수와 다중 지점의 중심점을 산출할 수 있는 centroid() 함수를 알아본다.

distGeo()

다음은 distGeo() 함수의 사용 예시이며 경도(longitude), 위도(latitude) 순서대로 입력해야 한다. 한 지점의 위경도 좌표는 c() 함수로 묶어주어야 하며 지점간 위경도 좌표는 쉼표로 구분해준다.

1
2
3
distGeo(c(128.1, 38.1),
c(128.2, 38.2))
## [1] 14143.53

결과는 미터단위로 산출되며 앞의 결과의 경우 두 지점간의 겨리가 약 14km 이상임을 알 수 있다. 그리고 기본 연산은 WGS84(WGS1984와 동일) 기준이며 distGeo() 함수의 장반경(張半徑) 인자 a와 편평률(扁平率, flattening) f는 각각 6378137(m)과 1/298.25223563 으로 설정되어있으며 필요시 입력 자료의 지리좌표계에 맞는 장반경과 편평률을 기입하도록 하자.

distGeo() 함수는 친절하게 벡터화 연산을 지원하기 때문에 다중 지점간의 거리를 데이터프레임으로 정리한 경우 다음과 같이 편리하게 사용할 수 있다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
df_pos = data.frame(lon = c(128.0, 128.1, 128.2),
lat = c( 38.0, 38.1, 38.2))
df_pos
## lon lat
## 1 128.0 38.0
## 2 128.1 38.1
## 3 128.2 38.2

distGeo(df_pos)
## [1] 14150.79 14143.53 NA

df_pos[, "dist"] = distGeo(df_pos)
df_pos
## lon lat dist
## 1 128.0 38.0 14150.79
## 2 128.1 38.1 14143.53
## 3 128.2 38.2 NA

centroid()

centroid() 함수는 다중 지점간의 좌표를 종합하여 중심 지점 좌표를 반환한다. 이 함수도 distGeo() 함수처럼 각 지점을 별도로 입력하는 방법과 데이터프레임 형식으로 입력하는 방법 둘 다 지원하기 때문에 다음 예시와 같이 편리하게 사용할 수 있다.

1
2
3
4
5
6
7
8
9
df_pos[, 1:2]
## lon lat
## 1 128.0 38.0
## 2 128.1 38.1
## 3 128.2 38.2

centroid(df_pos[, 1:2])
## lon lat
## [1,] 128.1 38.10005

함수에 입력된 좌표의 숫자만 보면 세 지점의 좌표가 일직선 상에 있는 것으로 보이지만 지표면 상의 지점이기 때문에 일직선 위에 각 지점이 놓여있다고 보기 어렵다. 그래서 결과를 보면 위도(latitude)의 값이 38.1 보다 약간 더 큰 38.10005가 산출된 것을 알 수 있다.


응용

이제 seoul_subway_line_2_geocoded.csv 파일을 다뤄보자.

1
2
3
4
5
6
7
8
9
df = read.csv("data/seoul_subway_line_2_geocoded.csv")
head(df)
## CD stn_nm_kr stn_nm_en line_no CD_ext loop order lon lat
## 1 200 까치산 Kkachisan 02호선 234-4 0 0 126.8467 37.53177
## 2 201 시청 City Hall 02호선 201 1 1 126.9769 37.56570
## 3 202 을지로입구 Euljiro 1(il)-ga 02호선 202 1 2 126.9827 37.56607
## 4 203 을지로3가 Euljiro 3(sam)-ga 02호선 203 1 3 126.9911 37.56629
## 5 204 을지로4가 Euljiro 4(sa)-ga 02호선 204 1 4 126.9977 37.56662
## 6 205 동대문역사문화공원 Dongdaemun History Culture Park 02호선 205 1 5 127.0090 37.56568

각 변수 설명은 다음과 같다.

● CD: 지하철역 코드
● stn_nm_kr: 지하철역명(국문)
● stn_nm_en: 지하철역명(영문)
● line_no: 지하철 노선명
● CD_ext: 지하철역 외부 코드
● loop: 내선순환역 여부
● order: 시청역을 기점으로한 내선순환역 배열 순서
● lon: 경도(longitude)
● lat: 위도(latitude)

2호선 내선순환역 정보와 필요한 변수만 남겨보자.

1
2
3
4
5
6
7
8
9
df_sub = df[df$loop == 1, -(3:6)]
head(df_sub)
## CD stn_nm_kr order lon lat
## 2 201 시청 1 126.9769 37.56570
## 3 202 을지로입구 2 126.9827 37.56607
## 4 203 을지로3가 3 126.9911 37.56629
## 5 204 을지로4가 4 126.9977 37.56662
## 6 205 동대문역사문화공원 5 127.0090 37.56568
## 7 206 신당 6 127.0195 37.56567

거리 계산

각 지점의 위경도를 한 번에 distGeo() 함수에 할당하여 지점간 거리 계산을 할 수 있다.

1
2
3
4
5
6
7
8
9
df_sub[, "dist"] = distGeo(df_sub[, c("lon", "lat")])
tail(df_sub)
CD stn_nm_kr order lon lat dist
39 238 합정 38 126.9145 37.55001 1213.0519
40 239 홍대입구 39 126.9245 37.55753 1097.9280
41 240 신촌 40 126.9365 37.55506 848.0974
42 241 이대 41 126.9459 37.55676 906.1749
43 242 아현 42 126.9561 37.55738 743.0818
44 243 충정로 43 126.9640 37.55965 NA

함수의 특성상 마지막 데이터는 결측이 되는데 다루는 지하철역은 순환 연결 구조이기 때문에 마지막의 결측치를 메꿔주자.

1
2
3
4
5
6
7
8
9
nrow(df_sub)
## [1] 43

df_sub[43, "dist"] = distGeo(df_sub[43, c("lon", "lat")],
df_sub[1 , c("lon", "lat")])
tail(df_sub, 2)
## CD stn_nm_kr order lon lat dist
## 43 242 아현 42 126.9561 37.55738 743.0818
## 44 243 충정로 43 126.9640 37.55965 1320.3468

중심 계산

2호선 내선순환 지하철역의 중심 지점 좌표는 다음과 같이 계산할 수 있다.

1
2
3
4
cent = centroid(df_sub[, c("lon", "lat")])
cent
## lon lat
## [1,] 126.9849 37.5229

시각화

앞의 내용을 종합하여 Google Map 기반 시각화를 해보자. 기초가 없는 사람은 앞의 포스팅을 참고하도록 하자.

GCP 가입 및 지오코딩 API 신청
R) 시각화 - Google 지도

앞에서 산출한 중심 좌표를 중심으로 그린 Google 지도는 다음과 같다.

1
2
3
4
5
6
7
8
9
mykey = "AIzaSyBpkJubAqtwhJx0YRpDXEVR8Y9puaE6zf8" # 본인 API 키로 변경 필요
library("ggmap")
register_google(key = mykey)

map = get_googlemap(center = as.numeric(cent),
maptype = "roadmap",
zoom = 12,
scale = 2)
ggmap(map)

서울 지도

이제 지도 위에 2호선 내선순환역 위치를 표기해보자.

1
2
3
4
5
ggmap(map) + 
geom_point(data = df_sub,
aes(x = lon, y = lat),
color = "#00AA00",
size = 5)

2호선 내선순환역 위치 표기
일단 그림은 그려졌는데 모든 역이 표기가 다 되지 않은 느낌을 받을 것이다. 위 코드를 실행했을 때 다음과 같은 경고 메세지를 확인했을텐데 이는 어떤 이유로 geom_point() 함수로 표기하려던 데이터 두 개(2 rows)가 유실되었다는 것을 말해준다.

경고메시지(들):
Removed 2 rows containing missing values (geom_point).

중심점을 조금 옮겨서 그려보자.

1
2
3
4
5
6
7
8
9
map2 = get_googlemap(center = c(cent[1, 1] + 0.012, cent[1, 2]),
maptype = "roadmap",
zoom = 12,
scale = 2)
ggmap(map2) +
geom_point(data = df_sub,
aes(x = lon, y = lat),
color = "#00AA00",
size = 5)

2호선 내선순환역 위치 표기(2차)
온전하게 표기된 것을 확인할 수 있다.

이제 각 역을 이어보자.

1
2
3
4
5
6
7
8
9
ggmap(map2) + 
geom_path(data = df_sub,
aes(x = lon, y = lat),
color = "#FF0000",
size = 2) +
geom_point(data = df_sub,
aes(x = lon, y = lat),
color = "#00AA00",
size = 5)

2호선 내선순환역 연결
시청역과 충정로역이 이어지지 않은 것을 볼 수 있다.

임시방편으로 데이터 마지막에 시청역 정보를 rbind() 함수로 붙여주면 올바르게 시각화 되는 것을 볼 수 있다.

1
2
3
4
5
6
7
8
9
10
ggmap(map2) + 
geom_path(data = rbind(df_sub, df_sub[1, ]),
aes(x = lon, y = lat),
color = "#FF0000",
size = 2) +
geom_point(data = df_sub,
aes(x = lon, y = lat),
color = "#00AA00",
size = 5) +
theme_void()

완전히 이어진 2호선 내선순환역

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×