앞선 단계에서 데이터에 대한 수집과 간단한 처리 부분은 완료했다. 이제 이번 포스트에서는 본 프로젝트의 핵심인 H3의 사용의 시작과 쓰인 함수에 대한 소개를 하고자 한다.

H3는 Uber에서 개발한 육각형 그리드 시스템으로 파이썬으로 쉽게 그리드 시스템을 다룰 수 있다.

카일님의 블로그의 말을 빌리자면 육각형 그리드의 장점은 인접 타일 간의 거리가 같다는 점이 기존의 정사각형 그리드보다 강점을 지니며, 육각형 타일의 경우 재구성이 사각형 타일처럼 정확히 이루어지진 않으나, 약간의 회전과 함께 재구성이 가능하다.

(출처: https://zzsza.github.io/data/2019/03/31/uber-h3/)

서론은 줄이고 바로 H3 설치부터 다뤄보도록 한다.

0. H3 및 기타 패키지 설치

H3 패키지와 관련 함수들을 사용하기 위해 설치해야한 패키지의 목록은 다음과 같다.

- h3 

- geopandas

- geojson

- simplejson

- branca (LinearColormap에 쓰임)

각 패키지가 pip를 통해 바로 설치가 되면 좋겠지만, geopandas의 경우 pip를 통해 설치를 하면 설치 오류가 발생했고 관련 문제가 자주 일어나는 듯 보였다.

이 곳(https://codedragon.tistory.com/9556)의 글을 따라서 pip 설치가 아닌 whl파일을 직접 내려받아 로컬에서 설치를 진행하는 방식으로 진행해야 하며, geopandas를 설치하기 앞서 설치해야하는 패키지가 몇개 존재한다...

1. pyproj

2. Shapely

3. GDAL

4. Fiona

5. geopandas

https://www.lfd.uci.edu/~gohlke/pythonlibs/

각각의 패키지에 대한 whl파일은 위의 링크에서 내려받을 수 있다. 각각의 운영체제와 python의 버전에 맞게 내려받는다.

(whl 파일을 통한 설치 방법)

명령 프롬프트에서 pip install ~~~.whl의 명령을 통해 설치 파일을 내려 받은 whl파일로 지정하여 설치를 진행한다. 나같은 비전공자의 경우 cmd에서의 디렉토리 변경이 익숙치 않아 번거로우므로 워킹 디렉토리에 파일을 옮긴 다음 바로 pip install을 진행하는것이 편할 것이다ㅎㅎ

[예시]

(User 폴더에 whl파일 이동) -> (명령 프롬프트 실행) -> pip install geopandas-0.7.0-py2.py3-none-any.whl

(설치를 마무리한 뒤에는 whl파일을 삭제해도 무관하다.)

from h3 import h3
import pandas as pd
from tqdm import tqdm_notebook
import folium
from folium import GeoJson
import branca.colormap as cm
from geojson import Feature, FeatureCollection
import simplejson as json

본 단계에서 사용할 패키지들을 전부 import 했다.

위의 패키지를 전부 문제 없이 설치했다면 이렇게 라이브러리를 import해서 오류가 없는지 확인하고 다음 단계로 넘어가자.

1. H3 제공 주요 메서드

*프로젝트 중 설치한 패키지와 버전에 따라 메서드의 이름과 인자 이름이 다른 것을 확인하였으나, 공식 Documentation이 없어 이를 명확히할 수 없었다..!

아쉽게도 H3 패키지 자체에 시각화나 분석에 바로 쓸수 있는 메서드는 찾기 어렵다. h3패키지가 자체적으로 제공하는 주요 메서드들을 보면 다음과 같다.

1. geo_to_h3 : 위경도 좌표와 resolution을 받아 그 좌표가 속한 육각 타일의 id를 제공하는 메서드

2. h3_to_geo : 육각 타일 id를 받아 위경도 좌표 값으로 변환해주는 메서드

3. hex_ring : 육각 타일과 거리 k를 받아 거리가 k에 해당하는 타일의 id의 집합(set)을 반환하는 메서드

    ex. k=2 -> 12개 타일 반환

4. k_ring : 육각 타일과 거리 k를 받아 거리가 k이하인 타일의 id의 집합(set)을 반환하는 메서드

    ex. k=2 -> 1(거리=0) + 6(거리 =1) +12(거리=2) -> 총 19개 타일 반환

이런 주요 메서드를 보면 알 수 있듯이 h3 패키지가 분석을 자체적으로 도와주기보다는, 기존 geocode를 해당 위치와 resolution에 해당하는 타일 id를 반환해주고, 타일의 위치나 인접 타일에 대한 정보 정도를 다뤄주는 정도까지의 기능만 제공해준다,

즉 우리가 직접 육각 타일에 값을 부여하고, 그 결과를 시각화하는 작업은 직접 수행해야 된다는 것이다.

이를 위해 다음 단계에서 분석을 위한 함수를 직접 정의하는 과정을 다뤄보자!

2. H3 패키지 메서드 활용 함수

상술한 H3 패키지의 메서드를 이용해 분석단계에 활용한 함수들은 아래의 글에서 참고하였으며, 일부 수정을 거쳤다.

medium.com/better-programming/playing-with-ubers-hexagonal-hierarchical-spatial-index-h3-ed8d5cd7739d

1. counts_by_hexagon

def counts_by_hexagon(df, resolution):
    df = df[["lat","lng"]] # 데이터프레임에서 lat, lng만 가져와 처리
    
    df["hex_id"] = df.apply(lambda row: h3.geo_to_h3(row["lat"], row["lng"], resolution), axis = 1) 
    # h3.geo_to_h3 메서드를 통해 위경도 좌표가 속한 그리드의 id를 입력받아 hex_id 칼럼에 저장
    
    df_aggreg = df.groupby(by = "hex_id").size().reset_index() # hex_id를 기준으로 groupby하여, 행의 개수를 저장
    df_aggreg.columns = ["hex_id", "value"] # 칼럼이름을 각각 hex_id와 value로 변경
    
    df_aggreg["geometry"] =  df_aggreg.hex_id.apply(lambda x: {"type" : "Polygon",
                                                               "coordinates":[h3.h3_to_geo_boundary(x,geo_json=True)]})
    return df_aggreg

가장 처음으로 활용하게 되는 함수이다.

우선 위도('lat')와 경도('lng')값을 갖고 있는 데이터프레임을 인자로 받아 위경도 좌표값이 위치하는 타일의 id를 변환하여 'hex_id' 칼럼에 저장한다.

그 다음 hex_id 칼럼으로 groupby하여 size메서드를 통해 각 hex_id의 카운트 값을 계산하여 'value'칼럼으로 저장한다.

또한 이 타일을 시각화하기 위하여 'geometry'칼럼을 만들어 그곳에 geojson형태의 자료형을 입력해주었다. type은 'Polygon', 좌표값은 h3의 h3_to_geo_boundary 메서드를 활용하여 해당 타일을 이루는 6개의 점의 좌표값 리스트를 사용하였다.

2. sum_by_hexagon

def sum_by_hexagon(df,column, resolution): #count_by_hexagon과 유사하나, 특정 칼럼값을 지정하여, 타일별 그 값의 합을 연산
    
    df["hex_id"] = df.apply(lambda row: h3.geo_to_h3(row["lat"], row["lng"], resolution), axis = 1) 
    
    df_aggreg = df.groupby('hex_id')[column].sum().reset_index()# hex_id로 groupby한 뒤 지정한 칼럼의 sum을 계산
    df_aggreg.columns = ['hex_id','value']
    df_aggreg["geometry"] =  df_aggreg.hex_id.apply(lambda x: {"type" : "Polygon",
                                                               "coordinates":[h3.h3_to_geo_boundary(x,geo_json=True)]})
    return df_aggreg

위의 count_by_hexagon을 응용한 함수로, 위경도 값이 있는 데이터프레임에서 특정 value가 담긴 값을 타일별로 총합하여 반환하는 함수이다. count_by_hexagon에서 groupby 집계 방식만 바꿔주었다.

3. mean_by_hexagon

def mean_by_hexagon(df,column, resolution): #sum_by_hexagon과 유사하나 총합이 아닌, 평균으로 집계함
    
    df["hex_id"] = df.apply(lambda row: h3.geo_to_h3(row["lat"], row["lng"], resolution), axis = 1) 
    
    df_aggreg = df.groupby('hex_id')[column].mean().reset_index()#hex_id로 groupby하여 지정한 칼럼값의 평균으로 집계
    df_aggreg.columns = ['hex_id','value']
    df_aggreg["geometry"] =  df_aggreg.hex_id.apply(lambda x: {"type" : "Polygon",
                                                               "coordinates":[h3.h3_to_geo_boundary(x,geo_json=True)]})
    return df_aggreg

mean_by_hexagon은 마찬가지로 count_by_hexagon에서 groupby 집계 방식만 평균으로 바꾸어서 만들었다,

4. hexagons_dataframe_to_geojson

#hex_id와 value로 이루어진 데이터프레임을 geojson형태로 전환하는 함수
def hexagons_dataframe_to_geojson(df_hex, id='hex_id', value="value", file_output = None):
    list_features = []
    
    for i,row in df_hex.iterrows():
        feature = geojson.Feature(geometry = row["geometry"] , id=row[id], properties = {"value" : row[value]})
        list_features.append(feature)
        
    feat_collection = geojson.FeatureCollection(list_features)
    
    geojson_result = json.dumps(feat_collection)
    
    #optionally write to file
    if file_output is not None:
        with open(file_output,"w") as f:
            json.dump(feat_collection,f)
    
    return geojson_result

5. choropleth_map

# hex_id별 value가 입력된 데이터 프레임을 코로플레스 시각화 지도로 바꾸어주는 함수
def choropleth_map(df_aggreg,geojson_data, value='val',location =[37.65,126.865] ,border_color = 'black', fill_opacity = 0.7, initial_map = None, with_legend = True, kind = "filled_nulls"):
    #colormap시각화를 위한 중간값 m도출
    min_value = df_aggreg[value].min()
    max_value = df_aggreg[value].max()
    m = round ((min_value + max_value ) / 2 , 0)
    
    if initial_map is None:
        initial_map = folium.Map(location= location, zoom_start=12)
    
    # linear - 선형적인 관계를 나타냄 많고 적음을 나타내게 함
    # outlier - outlier를 찾기 위한 choropleth - 
    # filled_nulls - linear와 같으나, 0인 값을 회색으로 칠함.
    if kind == "linear":
        custom_cm = cm.LinearColormap(['green','yellow','red'], vmin=min_value, vmax=max_value)
    elif kind == "outlier":
        custom_cm = cm.LinearColormap(['blue','white','red'], vmin=min_value, vmax=max_value)
    elif kind == "filled_nulls":
        custom_cm = cm.LinearColormap(['gray','green','yellow','red'], 
                                      index=[0,min_value,m,max_value],vmin=min_value,vmax=max_value)
    geojson_data = geojson_data
    folium.GeoJson(geojson_data,style_function=lambda feature: {
            'fillColor': custom_cm(feature['properties']['value']),
            'color': border_color,
            'weight': 1,
            'fillOpacity': fill_opacity}
    ).add_to(initial_map)
    #add legend (not recommended if multiple layers)
    if with_legend == True:
        custom_cm.add_to(initial_map)    
    
    return initial_map

chropleth_map함수는 4번 hexagons_dataframe_to_geojson 함수로 만든 geojson 자료에 담긴 타일 정보를 folium 지도에 시각화하는 함수이다. 이름에 있는 것과 달리 본 함수는 folium의 Choropleth 메서드를 활용하는 것이 아닌, GeoJson 메서드를 활용하여 지도 위에 Polygon을 올린다. 인자로 받아주는 geojson_data에는 각 타일의 id와 그 위치가 담긴 'geometry'값이 있는데, 이를 참고하여 지도 위에 육각형 타일을 그리고, 각 타일에 대한 색깔은 branca의 LinearColormap 메서드를 활용하여 값의 크기에 비례하여 색을 결정한다.

이때, kind인자로 colormap으로 linear, outlier, filled_nulls 세 가지를 선택할 수 있다. linear은 기본적인 형태로 값이 작을수록 초록색, 클수록 빨간색을 띠는 colormap이며, outlier는 outlier가 아닌 일반적인 값은 흰색이고 값이 큰 부분과 작은 부분을 각각 빨간색과 파란색으로 채워준다. filled_nulls는 기본적으로 linear과 같은 colormap이나, 0인 부분을 회색으로 채워넣었다. (원래 코드에서는 다른 색이었는데 수정해 주었음.)

이외의 인자로는 Map의 시작위치를 설정하는 location, 타일들의 모서리의 색깔을 설정하는 border_color(없애고 싶은 경우 None), 타일의 불투명도를 설정하는 fill_opacity 등이 있다.

그 외에 인자로 설정할 수 있는 부분들이 함수에 존재하나, 본 프로젝트에서는 동일하게 사용될 인자는 따로 인자로 빼놓지 않았다.

6. get_neighbors, affect_neighbors

#hex_id값을 입력하면 인근 타일을 반환하는 코드. ring_size값에 따라 범위가 늘어남
def get_neighbors(h3_code,ring_size):
    return list(h3.hex_ring(h3_code,k=ring_size))

#한 그리드의 value값을 인접타일에 특정 비율만큼 추가하는 함수
def affect_neighbors(df, ratio):
    df = df.set_index('hex_id')
    df_2= df.copy()
    
    for h_id in tqdm(list(df.index)):
        neighbor_list = get_neighbors(h_id,1)
        for neighbor in neighbor_list:
            if neighbor in list(df_2.index):
                df_2.loc[neighbor,'value'] += df.loc[h_id,'value']*ratio
            else:
                df_2.loc[neighbor] = {'value':df.loc[h_id,'value']*ratio,'geometry':{"type" : "Polygon","coordinates":[h3.h3_to_geo_boundary(neighbor,geo_json=True)]}}
    return df_2.reset_index()

get_neighbors 함수는 하나의 타일 id에 대하여 인접한 타일의 id의 리스트를 반환하는 함수이다. ring_size를 인자로 설정하긴 했으나, 프로젝트 내에선 k=1로만 사용하였다.

affect_neighbors 함수는 타일 데이터가 담긴 데이터 프레임(df)과 비율(ratio)을 인자로 받아 각각의 타일에 대하여 인접한 타일에게 본 타일에 매핑된 값을 분배해주는 함수이다. 

다음 포스트에서 위 함수들을 활용하여 앞서 수집한 데이터를 타일에 매핑하여 분석해 볼 것이다.

본 단계에서 부터는 실질적으로 지리데이터를 갖고 서울시 전기차 충전소 분석을 진행하고자 한다.

이 글에서는 간단하게 아래의 두 가지 데이터만 활용해보고자 한다.

1. 서울시 전기차 충전소 현황 데이터

2. 서울시 행정동별 생활인구 데이터

두 자료가 각각의 전기차 충전의 공급 현황과 수요를 나타내는 데이터라 가정하며, 1번의 경우 앞선 포스트에서 수집 방법에 대해서 다루었으므로, 이번에는 '서울시 행정동별 생활 인구 데이터'를 다루어 보고자 한다.

1. 서울시 행정동별 생활인구 데이터 다루기

서울 열린 데이터 광장(http://data.seoul.go.kr/)에는 서울시의 행정동별 생활인구 데이터를 월별로 제공하고 있다. 

(자료 : http://data.seoul.go.kr/dataList/OA-14991/S/1/datasetView.do)

본 분석에서는 내국인에 대한 20년 1월 자료를 활용하였다. 최근 데이터를 사용하지 않은 것은 코로나 바이러스로 인한 생활인구 자료에 노이즈가 발생하였을 것을 우려했기 때문인데, 사실 큰 차이는 느껴지지 않았다.

자료를 보면 행정동별, 일별, 시간대별, 연령대별, 성별별 생활인구가 자세하게 기록되어있는 것을 볼 수 있다. 또한 행정동의 경우 동 이름이 아닌, 행정동 코드로 입력된 것을 확인할 수 있는데, 이를 매핑하기 위해 행정동 코드 자료 또한 구해야 한다. 행정동 코드는 위의 주소와 같은 페이지에서 구할 수 있다.

행정동 코드 정보를 얻었으니 우리는 pandas를 활용하여 행정동 코드를 행정동 이름으로 변경하는 작업을 해보고자 한다.

import pandas as pd
pop_df = pd.read_csv('LOCAL_PEOPLE_DONG_202001.csv')
hd_code = pd.read_csv('행정동코드.csv')

#행정동 코드를 key, 행정동명을 value로 하는 dictionary 정의
hd_dict = dict(zip(hd_code['행자부행정동코드'],hd_code['행정동명'])) 

#위에서 정의한 hd_dict를 참고하여 '행정동'이라는 열 생성
pop_df['행정동'] = pop_df['행정동코드'].replace(hd_dict) 

위 코드를 통하여 '행정동'이라는 이름의 열로 행정동명 데이터를 생성하였다. 다음으로는 본 데이터를 분석에 맞게 reshaping하는 작업이 필요하다. 

1. 본 분석에서는 생활인구를 크게 주간(10시~17시)과 야간(19시 ~07시)의 두 시간대로 나누어 진행한다. (출퇴근 시간대인 07시~10시와 17시~18시는 생활 인구가 크게 이동하는  제외하였다.) 행정동별로 해당 시간대별 평균 생활인구를 구한다.

2. 두 시간대의 행정동별 평균 생활인구수를 구하면, 야간대비 주간 생활인구 증가율을 계산한다.

3. 성별 및 연령대별 데이터가 나뉘어져 있는데, 전기차 수요와 보다 관련이 있는 25~54세에 해당하는 남성 및 여성의 생활 인구에 대해서만 다룬다.

위의 세 가지 작업을 실행하는 코드는 아래와 같다.

pop_df_day = pop_df[(pop_df['시간대구분']>=10)&(pop_df['시간대구분']<=17)]
pop_df_night = pop_df[(pop_df['시간대구분']>=19)|(pop_df['시간대구분']<=7)]
col_list = ['남자25세부터29세생활인구수','남자30세부터34세생활인구수','남자35세부터39세생활인구수','남자40세부터44세생활인구수',
            '남자45세부터49세생활인구수','남자50세부터54세생활인구수','여자25세부터29세생활인구수','여자30세부터34세생활인구수',
            '여자35세부터39세생활인구수','여자40세부터44세생활인구수','여자45세부터49세생활인구수','여자50세부터54세생활인구수']
            
#행정동별 groupby 후 각 연령대 자료 총계 계산
night = pop_df_night.groupby('행정동').mean()[col_list].sum(axis=1)
day = pop_df_day.groupby('행정동').mean()[col_list].sum(axis=1)

#야간대비 주간 증가율 계산
night_day = day/night

#새로운 데이터프레임 pop_df 생성
pop_df = pd.DataFrame(night_day,columns=['야간대비 주간생활인구 증가율']).reset_index()
pop_df.loc[:,'주간'] = list(day)
pop_df['야간'] = list(night)

위 작업을 통해 얻은 데이터는 아래와 같다.

이제 우리는 두가지 데이터 모두 전처리가 완료되었다. 이제 H3 라이브러리를 활용한 분석을 할 수 있으면 좋겠지만..

우리가 가진 데이터는 모두 위경도 자료가 매핑되어있지 않다. 따라서 우리는 각 데이터에 있는 주소 데이터를 활용해 위경도 값을 먼저 찾아줘야 한다.

2. googlemaps 라이브러리를 활용한 위경도 값 찾기

import googlemaps
from tqdm import tqdm_notebook

def latlng_getter(df, address_col):
    gmaps_key = 'YOUR-API-KEY'
    gmaps = googlemaps.Client(key=gmaps_key)
    for i in tqdm_notebook(range(len(df))):
        address = df.loc[i,address_col]
        json_obj = gmaps.geocode(address,language='ko')
        if type(json_obj) == list: #json obj가 리스트로 불려올 떄도, 그대로 dictionary로 불려올 떄도 있어서 둔 코드
            json_obj = json_obj[0]
        df.loc[i,'lat'] = json_obj['geometry']['location']['lat']
        df.loc[i,'lng'] = json_obj['geometry']['location']['lng']
    return df
    
pop_df = latlng_getter(pop_df,'행정동')
sttn_df = latlng_getter(sttn_df,'address')

위경도 값을 매핑하기 위하여 위와 같은 함수 latlng_getter를 정의하였다. 

googlemaps 라이브러리를 활용하여 googlemaps API에서 행단위로 선택한 열값(주소)을 검색하여 얻은 위경도 값을 각각 'lat', 'lng' 칼럼을 생성하여 부여하였다.

이번 분석을 통해 얻은 것

1. 공공데이터를 활용해 행정동별 생활인구를 다루어보았고, 행정동 코드 자료도 받아 이를 매핑하는 과정을 경험함.