R의 좌표계 변환하기

좌표계 변환을 도와달라고 해서 잠깐 찾아보았다. 좌표계가 이렇게 다양한지도 몰랐고, 생각보다 정리되어 있는 곳이 많지 않다.

PROJ.4라는 좌표계 변환을 제공하는 라이브러리가 있으며 다양한 언어로 포팅되어 있다고 하는데, 비개발자에겐 접근이 쉬워 보이지 않는다. 재미삼아 R로도 되지 않을까 검색해봤다. 역시나 있다. 아니 자바스크립트로도 포팅되어 있는데 R에 없는게 이상하잖아;

아무튼 대충 구글링을 해서 예제를 만들어보았다.

1. 필요 라이브러리

library(sp)
library(rgdal)

2. 좌표 변환 함수

• parameter ** long : 경도 ** lat : 위도 ** from.crs : 갖고 있는 데이터 좌표계의 proj.4 인자 ** to.crs : 변환할 좌표계의 proj.4 인자
• return ** from.crs 좌표계 좌표를 to.crs 좌표계 좌표로 반환하여 data.frame 형식으로 반환한다.

좌표계에 해당하는 proj.4 인자는 이 곳에서 참고했다. * OSGeo 한국어 지부 한국 주요 좌표계 EPSG 코드 및 proj4 인자 정리

convertCoordSystem <- function(long, lat, from.crs, to.crs){
  xy <- data.frame(long=long, lat=lat)
  coordinates(xy) <- ~long+lat
  
  from.crs <- CRS(from.crs)
  from.coordinates <- SpatialPoints(xy, proj4string=from.crs)
  
  to.crs <- CRS(to.crs)
  changed <- as.data.frame(SpatialPoints(spTransform(from.coordinates, to.crs)))
  names(changed) <- c("long", "lat")
  
  return(changed)
}

3. 테스트

UTM-K(GRS-80) 좌표계에서 WGS84 경위도 좌표계로 변환한다.

coord <- data.frame(utmk.long=c(954677.6, 958869.4), utmk.lat=c(1951583, 1945669))
from.crs = "+proj=tmerc +lat_0=38 +lon_0=127.5 +k=0.9996 +x_0=1000000 +y_0=2000000 +ellps=GRS80 +units=m +no_defs"

to.crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
coord <- cbind(coord, convertCoordSystem(coord$utmk.long, coord$utmk.lat, from.crs, to.crs))

4. 결과

for(i in 1:nrow(coord)){
  p <- coord[i,]
  str <- paste("(", p$utmk.long, ",", p$utmk.lat, ") -> (", p$long, ",", p$lat, ")", sep="")
  print(str, quote=FALSE)
}

[1] (954677.6,1951583) -> (126.986819396889,37.5624908864062)
[1] (958869.4,1945669) -> (127.034612566292,37.5093829143664)

대충 구글링을 해보니 맞는 것으로 보인다… 만들고 나서 나도 확신을 못하겠지만 :( 아무튼 전문가분들이 보고서 고칠 점이 있으면 말해달라. 반영하겠음.