GIS 기반의 프로젝트를 수행중에 있는 어느 날
어떤 특정한 위경도 좌표가 입력될 때 대한민국의 시군구 중 어디인지를
알려주는 함수를 필요로 하게 되었다.
API 기반의 서비스를 이용하지 않고 순수 로직에 의해 만드는것을 목표로
하고 있어 다양한 요소에 개발을 필요로 하고 있는데
큼지막한 요소로는 2가지 미션들이 있다.
- 좌표계의 통일을 위한 변환방법
- 특정 (x, y) 좌표가 시군구 경계에 해당하는 컨백스홀 안에있는지 밖에있는지를 판단하는 논리 - 추후 포스팅 예정
이 중 첫번째 좌표계의 통일을 위한 변환방법을 R의 sp 패키지를 기반으로 알아보고자 한다.
우리나라에서 사용하는 좌표계의 종류 - 출처 : biz-gis.com
GIS에 대해서 1도 몰랐던 시절 좌표계의 종류가 이렇게 많은지 모르고
있었다.
더군다나 실제 GIS 데이터들은 이런 좌표단위들을 다양하게 사용하고 있어서
혼란이 많았다.
예를 들면 이렇다.
서울의 광화문 중앙의 좌표를 (위도,경도)1
라고 할 때
A 테이블에선 (37.5760336, 126.9769142) 로 사용중이고
B 테이블에선 (37.5749854, 126.9746067) 로 관리가 되고
있다.
미세한 차이로 보이지만 실제로 이 위경도값을 그대로 구글지도에 입력하여
찾아보면 멀리 떨어진 지점이다.
Google Map 에서 A 테이블의
(37.5760336, 126.9769142)를 찍어 보았을 때
Google Map 에서 B 테이블의
(37.5749854, 126.9746067)를 찍어 보았을 때
이런 예제처럼 A테이블 정보와 B테이블의 정보를 매쉬업 하기 위해선
단위를 통일하여야 매쉬업 할 때 왜곡되어지지 않는다.
자연스럽게 A, B중 하나는 좌표계를 변환할 필요가 생기는데 이
시점에서
R에서 sp package 를 이용한 좌표계 변환 성공예제를 정리하고 공유하고자
한다.
행정경계 자료 획득 - 시군구 별 경계
시군구를 구분짓는 경계에 대한 데이터 베이스가 필요했는데
이 데이터는 구했지만 여기서 부터 문제가 있었다.
필자에게 익숙한 경위도 지리좌표계가 아닌 다른 단위의 좌표계로 관리가
되고 있다는 점이었다.
문제의 시군구별 행정경계 데이터는 http://www.biz-gis.com/GISDB/ 에서
제공하고 있다.
다운로드 하는 방법은 http://www.biz-gis.com/GISDB/ 로 들어간 후
“행정경계-시군구” 항목을 클릭하여 다운로드 버튼을 누르면 된다.
“vADM2_TM.zip” 이라는 압축파일이 다운로드가 되고, 압축을 풀면 아래와 같은 파일들이 떨어진다.
## [1] "v시군구_TM.dbf" "v시군구_TM.sbn" "v시군구_TM.sbx" "v시군구_TM.shp"
## [5] "v시군구_TM.shx"행정경계정보가 담긴 이 5개 파일을 R에서 불러와 분석을 하기 위해
rgdal::readOGR() 함수를 이용한다. 2
library(rgdal)
d <- readOGR("data")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/lovetoken/Dropbox/02_Study/01_Statistic/31_R/03_Posting/02_Package/sp/proj4string/data", layer: "v시군구_TM"
## with 232 features
## It has 4 fieldsreadOGR() 함수의 사용법은 .shp 포맷을 포함한
기하정보들을 모두 가진 디렉토리를 첫번째 인자에 입력하면 된다.
만약 입력한 디렉토리 안에 “vADM2_TM.zip” 파일을 압축해제한 5개 파일 외에
다른 레이어 파일들이 있을때는 무엇을 기준으로 읽어들일지 알 수 없으므로
두번째 layer 인자에 명시하여 알려 줄 수 있다.
d <- readOGR("data", "v시군구_TM")readOGR() 함수를 통해 불러온 객체 d 를
axes = T 를 통해 축을 표현시켜 무작정 플롯팅을 시켜보면
아래와 같이 뜬다.
plot(d, axes = T)
시군구 경계가 구분되어 있는 대한민국 지도가 그려지는 것을 확인할 수
있다.
plot() 함수를 실행할 때 x, y 점들을 포인팅 하는 원리로써
이런 지도가 그려지는것은
d 객체 안에 경계선에 대해서 수많은 경계점들의 좌표정보를
가지고 있기 때문이란 것을 추측할 수 있고
직접 까보면 그 실체를 알 수 있게 된다. str(d) 출력의
일부(tail)를 확인해 보면
str(d) .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 179310 418281
.. .. .. .. .. .. ..@ area : num 4775
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:26, 1:2] 179312 179320 179328 179337 179341 ...
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 179400 418308
.. .. .. .. .. .. ..@ area : num 1161
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:24, 1:2] 179420 179419 179420 179414 179410 ...
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 179184 418418
.. .. .. .. .. .. ..@ area : num 3849
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:25, 1:2] 179181 179168 179160 179159 179152 ...
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 176061 418565
.. .. .. .. .. .. ..@ area : num 9792
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:39, 1:2] 176124 176126 176117 176108 176099 ...
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 178754 418658
.. .. .. .. .. .. ..@ area : num 9316
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:37, 1:2] 178730 178736 178751 178760 178761 ...
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 171224 419615
.. .. .. .. .. .. ..@ area : num 1090438
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:31, 1:2] 170917 170743 170551 170426 170368 ...
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 176879 420421
.. .. .. .. .. .. ..@ area : num 558583
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:47, 1:2] 177012 177056 177119 177237 177352 ...
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 188820 407056
.. .. .. .. .. .. ..@ area : num 692251527
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:5345, 1:2] 191196 191278 191338 191358 191399 ...
.. .. .. ..@ plotOrder: int [1:29] 29 14 1 27 28 17 20 7 16 2 ...
.. .. .. ..@ labpt : num [1:2] 188820 407056
.. .. .. ..@ ID : chr "97"
.. .. .. ..@ area : num 697667674
.. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
.. .. .. ..@ Polygons :List of 1
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 226596 433472
.. .. .. .. .. .. ..@ area : num 428700798
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:5130, 1:2] 233437 233541 233595 233829 233846 ...
.. .. .. ..@ plotOrder: int 1
.. .. .. ..@ labpt : num [1:2] 226596 433472
.. .. .. ..@ ID : chr "98"
.. .. .. ..@ area : num 428700798
.. .. [list output truncated]
..@ plotOrder : int [1:232] 113 121 191 116 189 195 117 208 112 199 ...
..@ bbox : num [1:2, 1:2] -11603 -42102 548505 569052
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "x" "y"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr NA“sp” 라는 클래스로 관리되는 d 객체 안에는 coords 라는
속성에 뭔지 모를 수많은 값들이 모여져서 다각형(Polygon)을 구성하게 되는
식이라 추측할 수 있다.
그런데 플롯팅 된 지도의 x, y 축 단위를 우리에게 낯익은 경위도 지리
좌표계(Geographic Coordinate)로 바꾸려면 어떻게 해야할까?
지금의 x, y축의 수치를 보면 십만단위 범위값까지 나오는데 이는 구 형의
지구를 평면으로 펼친 후 포인팅하기위한 평면직각좌표계(Projected
Coordinate)를 기준으로 측정할 때 나올 수 있는 단위라고 한다.
필자가 개발하려는 것은 경위도 값(ex. 37.5760336, 126.9769142)이 들어올
때 어떤 시군구에 속하는지를 맵핑하는 함수를 개발하려다 보니 지금 보는
d 객체의 평면좌표계 값을 경위도 지리좌표계로 변형할 필요가
생긴 것이다.
Coordinate Reference System
좌표변형을 위해선 두가지를 사전에 알아야 한다.
from Coordinate Reference System(fromCRS) 과 to Coordinate Reference
System(toCRS) 이다.
좌표계 시스템(CRS)에 대해선 전문적인 지식들이 필요로 하지만 이에 대한
언급을 자세히 하지 않겠다.3
재량을 통해4 fromCRS 를 알아낸 후 그에 맞게
약속된 문법에 맞추어 준비한다.
본 예제에서는 중부원점 기준 좌표로 명시되어 있어 그에 맞는 CRS 를
준비했다. 참고로 CRS() 함수를 이용하여 CRS 임을
명시한다.
from_crs = CRS("+proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000 +ellps=bessel +units=m")문법에 대해서 간단히 아는선에서 이야기하면
“+proj” 은 투영방법을 명시하는 부분이고
“+lat_0”, “+lon_0” 은 어떤 포인트를 중심으로 투영했는지
“+k” 는 구의 이심률
“ellps” 는 타원체의 규격을 정의하는 파라미터들이라고 한다.5
fromCRS 를 명시했으니 toCRS 도 필요하다.
to_crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")경위도 좌표로 바꾸기 위한 toCRS 를 알아낸 후 위처럼 명시했다. fromCRS 과 toCRS 가 준비되었으니 이제 변형만 하면된다.
좌표변형
R의 sp package 에서 proj4string(),
spTransform() 를 통해 좌표변형을 성공할 수 있었다.
우선 아까 준비한 fromCRS 정보 from_crs 를 아래처럼
d 객체에 반영할 것이다.
d 는 sp class 로 coordinate reference 값을 속성으로 따로
관리하는데 “coord. ref.” 부분을 보면 이전에는 NA 로 알 수
없는 상태이다.
d## class : SpatialPolygonsDataFrame
## features : 232
## extent : -11603.22, 548504.6, -42102, 569052 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## variables : 4
## names : 광역시코드, 광역시명, 시군구코드, 시군구명
## min values : 11, 강원도, 1111, 가평군
## max values : 50, 충청북도, 5013, 횡성군이를 fromCRS 로 변형한다.
coordinate reference 값만 변형하기 위해선 proj4string()
함수를 이용하면 된다.
proj4string(d) <- from_crs다시 d를 조회해보면
d## class : SpatialPolygonsDataFrame
## features : 232
## extent : -11603.22, 548504.6, -42102, 569052 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000 +ellps=bessel +units=m
## variables : 4
## names : 광역시코드, 광역시명, 시군구코드, 시군구명
## min values : 11, 강원도, 1111, 가평군
## max values : 50, 충청북도, 5013, 횡성군coordinate reference 값이 from_crs 에 담겨있는 string
으로 바뀐것을 확인할 수 있다.
이 상태에서 spTransform() 함수를 통해 toCRS 형식으로
변경한다.
d2 <- spTransform(d, to_crs)
d2## class : SpatialPolygonsDataFrame
## features : 232
## extent : 124.5918, 130.9408, 33.11139, 38.61427 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
## variables : 4
## names : 광역시코드, 광역시명, 시군구코드, 시군구명
## min values : 11, 강원도, 1111, 가평군
## max values : 50, 충청북도, 5013, 횡성군CRS 를 변경시킨 d2 객체를 확인해 보면
coordinate reference 값이 to_crs string 으로 바뀐것 뿐만
아니라
extent 에 명시된 xmin, xmax, ymin, ymax 값도 바뀐것을 눈치 챌 수
있을것이다.
d2를 이번엔 무작정 플랏팅 해 볼까?
plot(d2, axes = T)
이로써 경위도 지리좌표계로 변형된 것을 확인할 수 있다.
sp class 로 이루어진 GIS 테이블들을 좌표통일 후 매쉬업 할 수 있겠다는
생각이 든다.