[GIS]국토지리원, EPSG:5179(UTM-K), QGIS

개요

지인이 정밀 지도(건물, 길)데이터를 쓰고 싶은데, EPSG:5179로 되어있는 데이터는 대체 뭐냐고 자문을 구해왔습니다. 저도 EPSG:5187은 자주 쓰는데 EPSG:5179는 처음 다뤄봐서 그 내용을 정리하고자 합니다.

개념

우리가 대표적으로 많이 쓰는 좌표계는 EPSG:4326이랑 EPSG:3857입니다. EPSG:4326은 우리에게 익숙한 위경도 좌표계라서 통일해서 쓰기 편하고, EPSG:3857은 지도 데이터를 보는 것 + 지도에 필요한 부가 기능(길찾기, 영역 계산 등의 공간 연산)을 구현하기 위해 많이 쓰입니다. 두 좌표계는 정의 범위가 전 세계입니다. 그렇다 보니 지역의 정밀 데이터를 표현하기에는 정확도가 비교적 떨어집니다. 이를 방지하기 위해 특정 구역 내에서 정의된 좌표계들이 생겨나는데요, EPSG:5187은 대한민국의 중부 지역의 어떤 점을 원점으로 두고 한반도를 범위로 삼는, 한국에서 대표적으로 쓰이는 좌표계 입니다. 물론 한국에서도 여러 상황과 정밀도, 표현을 위해 다양한 좌표계를 선언해서 쓰고 있습니다.

EPSG:5179는 UTM-K라고도 불립니다. 추후 다른 글에서 UTM이 뭔지, 좌표계의 큰 분류에는 어떤 분류가 있을 지 다룰 예정입니다. 좌표계는 특정 변수(계수라고도 부릅니다)에 의해 정의할 수 있습니다.
EPSG:5179의 경우로 살펴볼까요?

EPSG:5179+proj=tmerc +lat_0=38 +lon_0=127.5 +k=0.9996 +x_0=1000000 +y_0=2000000 +ellps=GRS80 +units=m +no_defs 

여기서 주목할 점은 아래와 같습니다.

• proj=tmerc : 투영법입니다. 상세 내용은 Traverse Mercartor에서 확인할 수 있습니다.
• lat_0, lon_0 : 원점입니다.
• k = scale factor
• x_0, y_0 : (0,0)은 lat_0, lon_0에서 얼마나 떨어져 있는지 나타냅니다.
• ellps=GRS80 : 투영체를 WGS84로 쓰는데 여기서는 측량에 좀 더 적합한 GRS80을 사용했습니다.
• units=m : 단위는 metre입니다.

응용

다른 데이터와의 호환성을 위해 EPSG:5179로 표현된 데이터를 EPSG:3857로 변환하겠습니다.
GeoPandas를 통해 변환하는 방법도 있지만, GeoPandas로 변환하려면 GeoDataFrame으로 데이터를 로드해서 변환하다 보니 좀 번거롭습니다.
아래 코드는 GDAL 내의 ogr(벡터를 다루는 라이브러리)을 사용해서 변환하는 코드입니다.

코드

from osgeo import ogr, osr

# 입력 및 출력 파일 경로 설정
input_file = "input_5179.shp"
output_file = "output_3857.shp"

# 드라이버 설정
driver = ogr.GetDriverByName("ESRI Shapefile")

# 입력 데이터셋 열기
input_ds = driver.Open(input_file, 0)
'''
벡터 데이터는 여러 layer를 포함할 수 있습니다.
'''
input_layer = input_ds.GetLayer()

# 좌표계 설정
source_srs = osr.SpatialReference()
source_srs.ImportFromEPSG(5179)

target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(3857)

# 좌표 변환 객체 생성
transform = osr.CoordinateTransformation(source_srs, target_srs)

# 출력 데이터셋 생성
'''
DataSource가 파일, Layer가 파일에 속한 layer라고 생각하시면 됩니다.
'''
output_ds = driver.CreateDataSource(output_file)
output_layer = output_ds.CreateLayer("transformed_layer", target_srs, input_layer.GetLayerDefn().GetGeomType())

# 필드 복사
'''
attributes 들을 복사하는 과정입니다. 
'''
input_layer_defn = input_layer.GetLayerDefn()
for i in range(input_layer_defn.GetFieldCount()):
    output_layer.CreateField(input_layer_defn.GetFieldDefn(i))

# 피처 변환 및 복사
'''
여기서 말하는 feature는 geometry 정보를 가지고 있는 객체를 의미합니다.
'''
for feature in input_layer:
    geom = feature.GetGeometryRef()
    geom.Transform(transform)

    output_feature = ogr.Feature(output_layer.GetLayerDefn())
    output_feature.SetGeometry(geom)

    for i in range(feature.GetFieldCount()):
        output_feature.SetField(i, feature[i])

    output_layer.CreateFeature(output_feature)

# 메모리 해제
input_ds = None
output_ds = None

실전에 적용하기

문제는 국토지리원에서 받은 정밀지도(shp)를 qgis에 띄우면 좌표계가 제대로 인식되지 않습니다.
수동으로 할당해주면 해결할 수 있지만 성가십니다. 코드로 어떻게 할까요?

코드

from osgeo import gdal, osr, ogr

input_file = "{shp파일 위치}"
#output_file = "{변환 결과 저장할 위치}"

# 드라이버 설정
driver = ogr.GetDriverByName("ESRI Shapefile")

# 입력 데이터셋 열기
input_ds = driver.Open(input_file, 0)

'''
원본 CRS 정보 복사
'''
input_layer = input_ds.GetLayer()
input_spatial_ref = input_layer.GetSpatialRef()
source_srs = osr.SpatialReference()
source_srs.ImportFromWkt(input_spatial_ref.ExportToWkt())

'''
변환할 EPSG:5179 정보 설정
'''
target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(5179)

'''
layer 내의 geometry를 순회하면서 각 geometry 별로 osr의 transform 함수를 이용해 좌표계 변환이 들어갑니다. 
'''
transform = osr.CoordinateTransformation(source_srs, target_srs)

output_file = "/Users/hmjeong/Downloads/example_data/input/N3L_A0020000_A002_000000_edit.shp"
output_ds = driver.CreateDataSource(output_file)
output_layer = output_ds.CreateLayer("transformed_layer", target_srs, input_layer.GetLayerDefn().GetGeomType())

input_layer_defn = input_layer.GetLayerDefn()
for i in range(input_layer_defn.GetFieldCount()):
    output_layer.CreateField(input_layer_defn.GetFieldDefn(i))

for feature in input_layer:
    geom = feature.GetGeometryRef()
    geom.Transform(transform)

    output_feature = ogr.Feature(output_layer.GetLayerDefn())
    output_feature.SetGeometry(geom)

    for i in range(feature.GetFieldCount()):
        output_feature.SetField(i, feature[i])

    output_layer.CreateFeature(output_feature)

# 메모리 해제
input_ds = None
output_ds = None

이렇게 해서 qgis에 다시 파일을 로드하면 EPSG:5179로 인식하는 모습을 볼 수 있습니다. 그런데 원래 qgis가 가지고 있는 EPSG:5179는 아래와 같은데요, 변환 결과는 조금 다릅니다.
아래가 qgis에서 보여주는 EPSG:5179이구요,

Korea 2000 / Unified CS
속성
단위: 미터
정적 (평면 고정된 기준점에 의존)
천체(天體): Earth
메소드: Transverse Mercator
WKT
PROJCRS["Korea 2000 / Unified CS",
    BASEGEOGCRS["Korea 2000",
        DATUM["Geocentric datum of Korea",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4737]],
    CONVERSION["Korea Unified Belt",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",38,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",127.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",1000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",2000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Topographic mapping (small scale)."],
        AREA["Republic of Korea (South Korea) - onshore and offshore."],
        BBOX[28.6,122.71,40.27,134.28]],
    ID["EPSG",5179]]
Proj4
+proj=tmerc +lat_0=38 +lon_0=127.5 +k=0.9996 +x_0=1000000 +y_0=2000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

우리가 얻은 결과는 아래와 같습니다.

PROJCS["KGD2002 / Unified CS",
    GEOGCS["KGD2002",
        DATUM["Korean_Geodetic_Datum_2002",
            SPHEROID["GRS 1980",6378137,298.257222101],
            TOWGS84[0,0,0,0,0,0,0]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4737"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",38],
    PARAMETER["central_meridian",127.5],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1000000],
    PARAMETER["false_northing",2000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","5179"]]

이 정보는 epsg.io 에서 가져온 정보입니다.
KGD2002와 'Korea 2000 / Unified CS'은 무슨 차이가 있을까요?
KGD2002는 전세계 차원으로 통일성을 위해, 한국 지역을 중심으로 하는 좌표계들이 참고하는 geodetic을 2002년에 갱신한 버전입니다.
UNSD에서 KGD2002를 정의하고 협약을 맺고자 하는 내용
'Korea 2000 / Unified CS'은 한국 내에서 통상적으로 쓰던 정의 방법으로 보입니다.
'korea #### / Unified CS'의 형태로 년도별로 시리즈가 있습니다. Geodetic / 정의 년도의 차이로 보입니다.
이래나 저래나 문제없이 qgis에서 열릴 것입니다.

결론

EPSG:5179를 이해하기 위한 개념들을 가볍게 알아보고, 데이터의 interoperability를 위해 웹서비스에서 자주 쓰는 EPSG:3857로 변환하는 코드를 알아봤습니다.

참고

• proj 7계수로 정의된 EPSG:5179 : https://www.osgeo.kr/17
• EPSG 정의 정보 : https://www.epsg.io