대한민국 행정동 경계 좌표 추출 #2 - python > GeoJSON

대한민국 행정동 경계 좌표 추출 #1 - GeoJSON

위 글에서 QGiS로 WGS84 GeoJSON형식의 좌표를 추출해봤다. (원본 좌표계 값 그대로인 상태(UTM-K))

 

hp파일을 java, python등을 이용해서 경계 좌표 변환 & 추출 해보려고 한다.

편의상 python을 선택했다.

 

결론부터 말하면 python만으로의 행정동 좌표계 편환 추출은 할게 못된다 -_-

일단 모든 경우의 수에 대한 시도는 못해봤다.

각종(?) geometry library연동 등의 방법이 존재할 수 있겠으나 ... 일단 모르겠다. 안할란다. -_-

 

 

ogs2ogr를 이용해서 좌표계 변환된 ESRI shapefile로 만들고,

이를 python으로 json추출하는 방향으로 해보니 이 정도는 뭐 쓸만하....다.

 

 

1. QGiS 설치 후 PATH경로 추가

  - ogr2org, cs2cs 등 명령어 실행 위함

  - 환경변수 PATH 경로 추가 : C:\OSGeo4W64\bin

  - 기타 : GDAL_DATA 환경변수 추가 & 값 설정 (선택)

    변수 이름 : GDAL_DATA 

    변수 값 : C:\OSGeo4W64\share\gdal

 

 

2. 좌표계 변환 shape 파일 생성

  - (소스) EPSG:5181 (국가공간정보포털 제공 데이터)

   +proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000 +ellps=GRS80 +units=m +no_defs

  - (대상타켓) WGS84 (EPSG:4326)

   +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

ogr2ogr -s_srs "+proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000 +ellps=GRS80 +units=m +no_defs" -t_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" -f "ESRI Shapefile" -skipfailures 변환된.shp Z_SOP_BND_ADM_DONG_PG.shp --config SHAPE_ENCODING "EUC-KR" (설명) ogr2ogr \   -s_srs "+proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000 +ellps=GRS80 +units=m +no_defs" \   -t_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" \   -f "ESRI Shapefile" \   -skipfailures  \   변환된.shp <--- 타켓   Z_SOP_BND_ADM_DONG_PG.shp <-- 원본소스   --config SHAPE_ENCODING "EUC-KR" <-- 원본이 EUC-KR (한글 깨짐 방지)

  [참고] proj.4 파라미터 관련

    http://proj.maptools.org/gen_parms.html

    https://www.osgeo.kr/17

 

 

3. python 코드 작성 & 실행 : 값 추출 & GeoJSON 파일 생성

  - transfromToWGS84() :

   국가공간정보포털이 제공하는 shp의 좌표계는 UTM-K이니 이를 WGS84로 변환해야 해서 해본건데

   절대 하면 안된다!!!! 너무 너어무 너어어어어무 느리다. 쓸게 못된다.

   허긴... 방대한 데이터를(?) 이렇게 루핑 처리한다는 것 자체가 느리긴 느릴테지... 

   다른 방식의, 효과적인 방식이 있는지 모르겠는데(geopandas, fiona, gdal 연동?),

   아직 깊게(?) 들어가진 말자. 암튼 의도는 이러하고(?) 쓰지 말자.

  - 프로퍼티에 OBJECTID등을 포함시키고자하면 zip으로 묶는 dict구조에 해당 파라미터를 포함시키면 된다.

# -*- coding: utf-8 -*- 
from pyproj import Proj, transform
import numpy as np
import pandas as pd
import shapefile
from json import dumps

# shapefile 읽기
reader = shapefile.Reader("E:\Geo\센서스경계)행정동경계\Z_SOP_BND_ADM_DONG_PG\변환된.shp", encoding="euc-kr")
fields = reader.fields[1:]
field_names = [field[0] for field in fields]
buffer = []

# 좌표계 변환 (skip)
proj_UTMK = Proj(init='epsg:5181')
proj_WGS84 = Proj(init='epsg:4326')
def transformToWGS84(shape):
    shape_type = shape.__geo_interface__["type"]
    coords_utmk = shape.__geo_interface__["coordinates"]
    coords_wgs84 = []
    i = 0
    for cds in coords_utmk:
        for p in cds:
            coords_wgs84.append([transform(proj_UTMK, proj_WGS84, p[0], p[1])])
    if shape_type == "Polygon":
        return dict(type=shape_type, coordinates=[coords_wgs84])
    elif shape_type == "MultiPolygon":
        return dict(type=shape_type, coordinates=[[coords_wgs84]])
    
# geojson 데이터 구성
for sr in reader.shapeRecords():
    print(sr.record[2])
    atr = dict(zip(field_names, sr.record))
    geom = sr.shape.__geo_interface__
    #geom = transformToWGS84(sr.shape)
    buffer.append(dict(type="Feature", properties=atr, geometry=geom)) 

# GeoJSON 생성 (ensure_ascii=False : 한글깨짐 방지)
geojson = open("filename.geojson", "w")
geojson.write(dumps({"type": "FeatureCollection", "features": buffer}, indent=2, ensure_ascii=False) + "\n")
geojson.close()

 

4. 추출 결과 (GeoJSON)

- 행정동 코드, 행정동명, 

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "BASE_YEAR": "2018",
        "ADM_DR_CD": "1101053",
        "ADM_DR_NM": "사직동"
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              126.96893892002558,
              37.57814630501477
            ],
            [
              126.96895512144755,
              37.577935265875205
            ],
            [
              126.96921284490557


            [ ....  
            ],
            [
              126.96893892002558,
              37.57814630501477
            ]
          ]
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "BASE_YEAR": "2018",
        "ADM_DR_CD": "1101054",
        "ADM_DR_NM": "삼청동"
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [
              126.97713556993105,
              37.597683560300666
            ],
            [....
            ....
            ],
            
....
    {
      "type": "Feature",
      "properties": {
        "BASE_YEAR": "2018",
        "ADM_DR_CD": "3902062",
        "ADM_DR_NM": "예래동"
      },
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [
              [
                126.37313696693445,
                33.22901600038364
              ],
              [
                126.37312154153095,
                33.22901295010255
              ],
              [
...
              ],
              [
                126.46328170099486,
                33.3596931869656
              ],
              [
                126.46396143346223,
                33.3598578961412
              ]
            ]
          ]
        ]
      }
    }
  ]
}