[Python] 파이썬으로 GIS하기 : GeoPandas로 간단한 지오프로세싱 수행하기 (좌표계 설정, join, spatial join, intersect, 지도 그리기)

0. 준비물

실습 시 사용하는 데이터는

(1) 집계구 영역 shp파일, (2) 행정동 shp파일 (3) 집계구 단위 서울 생활인구 데이터 (csv)

세 가지다.

이 shp 파일의 투영좌표계는 UTM-K (EPSG:5179)이다.

집계구.zip

6.54MB

행정동.zip

1.32MB

https://data.seoul.go.kr/dataList/OA-14979/F/1/datasetView.do

열린데이터광장 메인

생활인구 데이터는 여기에서 내려받는다.

 

1. GeoPandas로 shp 파일 불러오기

import geopandas as gpd
shape_oa = gpd.read_file("Geopandas/OutputArea.shp", encoding = "CP949")

지오판다스에서는 read_file이라는 함수를 사용해 shp 파일을 불러온다.

type 함수를 활용하여 방금 불러온 shp 파일을 geodataframe으로 인식하는지 확인한다.

 

19153개의 집계구로 이루어져있음을 확인할 수 있다.

간단히 해당 shp의 형상을 확인할 수도 있다.

 

2. GeoPandas로 각 폴리곤의 센트로이드 구하기

shape_oa["geometry"].centroid

 

를 이용하면 아래와 같이 각 집계구 피처 별 centroid를 구할 수 있다.

각 centroid의 x좌표와 y좌표를 구해서 원본 집계구 파일의 x, y열(새로 만듦)에 넣어주자.

이렇게 x, y열이 새롭게 만들어졌으며, 그 안에 각 피처의 centroid 좌표가 입력되어있음을 알 수 있다.

 

3. GeoPandas로 투영 좌표계 입력 및 투영 좌표계 변환하기

ArcGIS에서는 투영을 지정할 때는 'define projection' 툴을, 투영을 변환할 때는 'projection' 툴을 사용한다.

그와 마찬가지로 지오판다스에서는 투영을 지정할 때는 'set_crs'를, 투영을 변환할 때는 'to_crs'를 이용한다.

shape_oa.crs

를 입력하면 해당 데이터의 좌표를 알려준다.

만약 아무것도 나오지 않으면 투영이 아직 정의되지 않았다는 뜻이다.

 

shape_oa.geometry = shape_oa.geometry.set_crs("EPSG:5179")

set_crs를 이용해 좌표계를 정의한다.

본 집계구 shp의 기본 좌표는 EPSG 5179라고 한다.

좌표계 코드에 대해 자세한 건 아래 참조

https://kimhongsi.tistory.com/entry/GIS-%ED%95%9C%EA%B5%AD%EC%97%90%EC%84%9C-%EC%9E%90%EC%A3%BC-%EC%93%B0%EB%8A%94-%EC%A2%8C%ED%91%9C%EA%B3%84-%EC%A2%8C%ED%91%9C-%EC%BD%94%EB%93%9C-%EC%A0%95%EB%A6%AC 

[GIS] 한국에서 자주 쓰는 좌표계, 좌표 코드 정리

좌표계를 지정해주었으니 이렇게 crs를 입력했을 때 좌표계 정보가 나온다.

shape_oa.geometry = shape_oa.geometry.to_crs("EPSG:4326") #WGS84 경위도 좌표

to_crs를 이용해 좌표계를 변환한다.

EPSG : 4326으로 바꾼다. 이는 WGS84 좌표계이므로, 경위도로 바뀔 것이다.

 

crs를 확인하니 WGS84로 바뀌었음을 알 수 있다.

 

4. Pandas로 데이터 가공을 해보자

df=pd.read_csv("Geopandas/LOCAL_PEOPLE_20221120.csv", encoding = "CP949")

위에서 내려받은 생활인구 csv 파일을 불러온다.

Pandas에서는 read_csv를 이용해 표를 불러올 수 있다.

생활인구에 대한 정보가 표시된다.

참고로 *은 3명 이하라는 뜻이다. 특정이 쉬우므로 3명 이하면 0명/1명/2명/3명 관계 없이 모두 *로 처리한다.

 

이 파일은 모두 11월 20일에 대한 정보를 담고 있으므로 기준일 열을 삭제한다.

Pandas에서 열 삭제는 del을 활용한다.

 

del df['?"기준일ID"']

df = df.replace("*", "0")  # *을 0으로 치환

*이 모두 0으로 바꾸었다.

Pandas에서 A를 찾아 B로 바꾸어주는 함수는 replace이다.

 

df.dtypes

 

*이 있었다는 것에서 알 수 있듯, 현재 본 데이터는 문자열(object)로 되어있다.

우리는 이 생활인구수를 숫자로 바꾸어야 한다.

따라서 시간대구분, 행정동코드, 집계구코드를 제외한 나머지를 다 float로 바꾸어 줄 것이다.

 

시간대구분, 행정동코드, 집계구코드를 제외하기 위해 이 세 열을 index로 바꾸어준다.

df = df.set_index(['시간대구분', '집계구코드', '행정동코드'])

그럼 이와 같이 세 열이 index로 바뀌었다.

이제는 type을 일괄 변경할 수 있다.

 

df = df.astype(float)

이제 생활인구수 데이터가 모두 float 64bit 형태로 바뀌었다.

다시 index 처리한 것을 풀어주자. 

df = df.reset_index()

다시 원래대로 돌아왔다. 물론 세 열은 여전히 object 파일 형식을 유지하고 있다.

집계구 코드처럼 수치를 나타내지 않는데 숫자 자체가 길어버리면 문자열로 유지하는 게 낫다.

왜냐하면 만약 int(32비트 기준)로 type설정을 해주게 되면,  overflow(범위 이탈)되어 집계구코드 숫자가 음수로 바뀌기 때문이다.

 

이제 우리가 필요로 하는 요소만 뽑아내보자 (pivot). 

 

우리는 젊은 여성 생활인구를 파악하기 위해, "여자25세부터29세생활인구수" 열을 사용하고자 한다.

 

이렇게 20대 후반 여성들이 각 집계구별로, 각 시간대 별로 몇 명이나 있는지 표시하는 표로 만드는 걸 목표로 한다.

 

df_subset= df[  ["시간대구분", "집계구코드","여자25세부터29세생활인구수"]   ] #대괄호 안에 list의 형태로 넣어주어 여러 열을 불러옴

먼저 내가 원하는 데이터 (시간대 구분, 집계구 코드, 20대 후반 여성 생활인구) 세 가지를 포함한 subset을 만든다.

 

*******************************************

만약 노인인구만 추리고 싶다면

 

df_subset= df[  ["시간대구분", "집계구코드","여자60세부터64세생활인구수","여자65세부터69세생활인구수" , "여자70세이상생활인구수",  
                "남자60세부터64세생활인구수","남자65세부터69세생활인구수" , "남자70세이상생활인구수",  ]   ] #대괄호 안에 list의 형태로 넣어주어 여러 열을 불러옴

df_subset['노인인구'] = df_subset['여자60세부터64세생활인구수'] + df_subset['여자65세부터69세생활인구수'] + df_subset['여자70세이상생활인구수'] + df_subset['남자60세부터64세생활인구수'] + df_subset['남자65세부터69세생활인구수'] + df_subset['남자70세이상생활인구수']

del df_subset['여자60세부터64세생활인구수']
del df_subset['여자65세부터69세생활인구수']
del df_subset['여자70세이상생활인구수']
del df_subset['남자60세부터64세생활인구수']
del df_subset['남자65세부터69세생활인구수']
del df_subset['남자70세이상생활인구수']

 

*******************************************

이를 이용해 pivot table을 만들어보자.

 

pd.pivot_table(df_subset, index = "집계구코드", columns="시간대구분")["여자25세부터29세생활인구수"]

그 결과 해당 일자에 시간대 별로 몇 명의 20대 후반 여성이 해당 집계구 위치에 있었는지 정리한 표가 완성된다.

 

이를 바로 저장해도 되지만, 우리는 shp파일과 join시킬 것이므로 칼럼명을 조정해준다.

칼럼명은 일반적으로 글자 수 제한이 있고 (13 characters), 숫자로 시작하면 안 된다.

그러니 위의 column 이름을 0은 T00to01,  23은 T23to24와 같이 바꾸어준다.

위의 코드로 0은 00으로, 1은 01로, 23은 23으로 바꿀 수 있으며

df_subset["시간대구분"].astype(int).astype(str).str.zfill(2)  +"to" +( df_subset["시간대구분"].astype(int)+1).astype(str).str.zfill(2)

를 해주면 

가 되고

df_subset["시간대"] = "T" + df_subset["시간대구분"].astype(str).str.zfill(2)  +"to" +( df_subset["시간대구분"].astype(int)+1).astype(str).str.zfill(2)

를 하면

 

로 시간대라는 열이 새로이 추가된다.

 

pivot = pd.pivot_table(df_subset, index = "집계구코드", columns="시간대")["여자25세부터29세생활인구수"]

pivot이라는 표를 이와 같이 만들었다.

이제 데이터 가공은 마쳤으니, shp파일과 join시켜 위의 데이터를 공간 상으로 올려보자.

 

5. Pandas의 merge 함수로 표와 shp파일을 join하기 

- 전처리 : 집계구코드의 type을 문자열로 바꾸기

집계구 shp의 집계구코드 데이터의 유형은 문자열(object)인데에 반해,

현재 pivot의 집계구코드 열 type은 int 64비트인 것으로 나온다.

이를 문자열로 고쳐주자.

 

pivot["집계구코드"] = pivot["집계구코드"].astype(np.int64).astype(str)   #너무 커서 int 32비트로는 안 됨

만약 pivot 집계구코드 열의 type이 int32로 나온 경우에는 위의 코드를 쓰고,

아닌 경우

pivot["집계구코드"] = pivot["집계구코드"].astype(str)

를 써서 문자열로 바꾸어준다.

- merge로 join하기

 

이제 집계구 shp 파일과 이 pivot dataframe을 join해보자.

Pandas에서는 merge라는 함수로 join할 수 있다.

shape_oa = pd.merge(shape_oa, pivot, left_on = "TOT_REG_CD", right_on = "집계구코드", how = "left")
#join의 기준이 되는 key 칼럼 이름이 같으면 그냥  on, 다르면 left right.
#적어도 shp 파일이 다 유지되도록 left join을 주로 함

merge ( 집계구 shp, 데이터프레임, left_on = 집계구 shp 기준 집계구코드 열, right_on = 데이터프레임 기준 집계구코드 열, how = "left")

key 칼럼이란 두 개의 데이터에서 공통적으로 겹쳐, 장차 join을 할 때 기준이 되는 열을 말한다.

여기에서는 집계구코드를 기준으로 join시키기 때문에 집계구코드 열이 key 칼럼이라 할 수 있다.

 

집계구 shp에서는 집계구코드가 TOT_REG_CD라는 이름으로, 데이터프레임에서는 집계구코드라는 이름으로 되어있으니,

left_on과 right_on에 각각 해당 열 이름을 입력한다. 

 

원래는 공간 위상만을 담고 있던 집계구 shp파일에 이와 같이 생활인구 데이터가 join되었음을 알 수 있다.

 

 

6. GeoPandas의 sjoin 함수로 spatial join하기

현재는 집계구별 생활인구를 알고 있는 상황인데,

문제는 우리가 집계구코드만으로는 이곳이 어디인지 잘 모른다는 것이다.

따라서 행정동 단위 shp파일과 집계구 단위 shp 파일을 spatial join하여, 각 집계구 행에다 해당 집계구가 속하는 행정동 정보를 붙여주려 한다.

 

shape_dong = gpd.read_file("Geopandas/ADM_DONG.shp", encoding="CP949")

행정동 shp 파일을 이와 같이 불러온다.

이렇게 시군구 코드 및 행정동 코드가 있는, 우리에게 익숙한 형태의 데이터이다.

 

shape_dong.geometry = shape_dong.geometry.set_crs("EPSG:5179")

shape_dong.geometry = shape_dong.geometry.to_crs("EPSG:4326")

위에서 했던 set_crs와 to_crs를 순서대로 해준다.

그럼 이와 같이 경위도 좌표로 바뀐 것을 알 수 있다.

 

이제 행정동 shp와 집계구 shp파일을 좌표계가 서로 같도록 통일해준 상태이다.

 

위의 코드로 spatial join에 대한 자세한 설명을 볼 수 있다.

 

shape_oa = gpd.sjoin(shape_oa, shape_dong, how = "left")

이와 같이 sjoin함수로 left 방식의 spatial join을 수행한다.

집계구에다가 해당 집계구가 속하는 행정동을 갖다 붙이는 형태이다.

 

이제 집계구 shp파일에 

집계구의 위상정보 + 집계구별 25~30 여성 생활인구 + 행정동 정보까지 붙은 상황임을 알 수 있다.

 

7. Pandas의 plot 기능으로 간단한 지도 만들기

간단하게 지도를 볼 수 있다.

물론 상세한 symbology 수정을 하거나 다양한 데이터를 추가하고 싶다면 QGIS나 ArcGIS를 이용하는 것이 좋다.

혹은 파이썬에서 mapclassify를 이용하면 단계구분도를 깔끔하게 그려낼 수 있다.

 

2022년 11월 20일 출근길 시간의 20대 후반의 여성은 마곡 일대에 많이 분포하는 것을 알 수 있다.

 

8. to_file하여 shp 파일 내보내기 

shape_oa.to_file("result.shp", encoding = "CP949")

를 하면 현 워킹디렉토리로 원하는 이름으로 저장이 되는데

 

이렇게 에러가 뜰 수 있다.

del shape_oa["집계구코드_x"]     #한글이 들어있으면 저장이 안 될 수도 있어서 집계구코드를 지워주었음
del shape_oa["집계구코드_y"]     #한글이 들어있으면 저장이 안 될 수도 있어서 집계구코드를 지워주었음

그래서 이렇게 집계구코드 열을 지워주었더니

warning만 뜨고 저장이 잘 되었다.

이 result 파일을 ArcMap에서 불러오면

 

이렇게 잘 불러와진다.

이렇게 속성테이블도 잘 열린다.

 

심볼을 입히면 이렇게 되고

 

지도의 요소를 추가한 결과는 아래와 같다.

 

 

 

개인적인 소감으로는, ArcMap 혹은 ArcGIS Pro에서 데이터 관리를 하기 까다로우므로

앞으로는 Python의 GeoPandas를 애용해야겠다는 생각이 들었다!

 

 

 

부록 1. GeoPandas의 overlay 함수로 intersect 하기

 intersect를 해주어 행정동을 집계구로 잘라낼 수 있다.

gpd.overlay(shape_dong, shape_oa, how = "intersection") #동을 집계구로 나눔

지오판다스의 overlay 함수를 이용해 intersection을 해주자.

그 결과, 이렇게 하나의 집계구임에도 여러 개의 동으로 이루어진 부분이 있음을 알 수 있다.

 

부록 2. GeoPandas에서 각 피처 별 면적 구하기

shape_oa.geometry.area

로 각 피처별 면적을 구할 수 있다.