COMPAS 공모전에 진심인 한 사람으로서.. 공모전 시즌마다 지리데이터를 많이 다루게 되었는데
나는 지리데이터 분석이 재밌고 시각화 해놓으면 보기도 예뻐서 좋다 ✌✌
지리데이터 분석에서 자주쓰는 방법을 잘 정리해두면 앞으로의 공모전에 좋을 것 같아서 공모전 하나를 마친 뒤 이 글을 쓰기 시작했다 ! (최우수상 언제 한번 타보냐 ㅠㅠ)
지도 시각화 툴은 folium을 자주 사용하고, 그 베이스는 당연히 geopandas가 된다. (사실 이전 공모전에서 plotly 지도도 그려 보려고 했었는데 잘 안됐다 😥)
그래서 이제부터 이 두 가지를 주로 활용하는 분석 방법을 정리해보고자 한다.
1. 위도, 경도가 있는 csv 데이터를 geopandas 형식으로 변경하기
위도, 경도가 컬럼으로 주어진 csv 데이터를 받을 때가 있는데, 이 csv를 geopandas 형식으로 변경해서
위도, 경도 기반으로 geometry를 만들어 볼 수 있다.
이러한 변경방식은 하도 자주 사용하다보니 함수화해서 편하게 사용하고 있다 :)
from shapely.geometry import Point, Polygon, LineString
import geopandas as gpd
def geo_transform(DataFrame) :
# csv to geopandas
# lon, lat data를 geometry로 변경
DataFrame['lat'] = DataFrame['lat'].astype(float)
DataFrame['lon'] = DataFrame['lon'].astype(float)
DataFrame['geometry'] = DataFrame.apply(lambda row : Point([row['lon'], row['lat']]), axis=1)
DataFrame = gpd.GeoDataFrame(DataFrame, geometry='geometry')
DataFrame.crs = {'init':'epsg:4326'}
DataFrame = DataFrame.to_crs({'init':'epsg:4326'}) # 좌표계 epsg : 4326
return DataFrame위도 경도에 맞는 Point geometry를 생성하고
우리나라의 표준 좌표계(crs)인 epsg:4326에 맞게 crs를 설정해주는 과정이다.
이 함수를 적용하면 그림과 같은 결과가 나온다.
2. 특정 point에 대한 일정 반경 지도에 표현하기
나는 공모전에서 cctv 커버 면적을 표현할 때 이 반경 표현법을 사용했는데, 지도에 표현할 때는
folium의 Circle을 사용해서 원하는 크기(반경)의 원을 그리면 된다.
# CCTV 반경 50m 그리기
# map 생성
m = f.Map(location=[lat, lon], zoom_start=11, tiles='CartoDB positron')
for _, row in df.iterrows() :
f.Circle(location=(row['lat'], row['lon']), popup=row['type'], radius=50, color='#FF580B',
fill='#FF580B').add_to(m)
m # 지도 표현folium Map의 location에는 다음과 같이 원하는 중심 좌표 (lat, lon)를 쓰면 된다.
location = [35.2721639, 128.8479674]
그리고 Map의 타일 형식에는 여러가지가 있는데, default는 openstreetmap이고, 다른 타입의 타일을 적용하고 싶으면 다음 링크를 참고하면 된다.
[참고 링크] folium-map-tiles
본 코드는 geometry 정보가 있는 데이터프레임(df)의 반복문을 돌면서 지도에 lat, lon 위치에 대한 50m 반경(radius)의 원을 표현하는 것이다.
popup 옵션을 설정하면 해당 위치에 마우스를 대면 설정한 값의 팝업이 뜨며, color와 fill 옵션은 각각 원의 색상, 원을 채우는 색상이다. 지도에 다음 그림과 같은 방식으로 표현된다.
지도 표현 뿐만 아니라 반경 영역을 아예 geometry로 만들어서 사용하고 싶다면, buffer를 활용해서 중심좌표 기준 일정 반경의 원을 만들 수 있다.
temp = df.to_crs({'init':'epsg:5179'}).buffer(100).to_crs({'init':'epsg:4326'})
buffer_geometry = gpd.GeoDataFrame({'geometry' : temp})buffer를 생성할 때는 좌표계 epsg:4326 대신 epgs:5179를 사용한다. 그럼 다음 그림과 같이 우리가 원하는 미터 단위의 buffer가 잘 생성된다. (왜인지는 모르겠고 반경 지오메트리 생성법 JONNA 열심히 구글링 하다 찾았다)
3. 특정 영역(폴리곤)의 면적(area) 계산하기
multipolygon과 같은 geometry 데이터의 면적을 구하고 싶을 때가 있다. geopandas에서 .area를 하면 면적이 구해지긴 하는데, 이게 좌표계 epsg:4326 기준으로 하면 우리가 생각하는 그 정상적인 미터, 키로미터 단위의 수가 나오지 않으므로, 좌표계를 epsg:6933으로 변경해주어야 한다. (왜인지는 모르겠고 이것도 멀티폴리곤 면적 계산법 JONNA 열심히 구글링 하다 찾았다)
# 제곱 킬로미터 기준(100*100m은 0.01km^2, 10000m^2) : 10 ** 6
df['area'] = df['geometry'].to_crs({'init' : 'epsg:6933'}).map(lambda x : round(x.area, 2) ) 기본은 제곱 미터 기준인거 같은데, 제곱 키로미터 단위를 하고자 한다면 여기서 나누기 10 ** 6을 하면 된다.
참고로 지구는~ 둥그니까 계산된 면적의 수치에 약간씩 오차가 있을 수 있다.
4. 특정 영역들에 대하여 합집합, 교집합, 차집합 구하기
두 개의 geopandas 데이터프레임이 존재할 때, 이 두 개의 geometry에 대한 합집합, 교집합, 차집합을 구하고 싶다면?
geopandas의 overlay가 그 기능을 지원해 준다.
# 합집합
result = gpd.overlay(df1, df2, how='union')
# 교집합
result = gpd.overlay(df1, df2, how='intersection')
# 차집합
result = gpd.overlay(df1, df2, how='difference')
# 대칭 차집합
result = gpd.overlay(df1, df2, how='symmetric_difference')추가로 이런 작업을 하면 geometry가 여러개로 쪼개지는 경우가 있는데..(예를 들면 교집합을 해서 중간중간의 geometry가 겹치는 경우 결과 데이터프레임의 row가 여러개 만들어진다)
이걸 다시 합치려면 .dissolve() 를 사용하면 된다.
dissolve() 내 파라미터로는 묶고자 하는 기준 컬럼(예를 들면 격자 넘버 컬럼?)을 설정해주면 그걸 기준으로 묶이게 된다.
5. 히트맵 그리기
데이터에 위도 경도 포인트별 수치로 표현되어 있다면, 폴리움의 히트맵으로 지도에 표현해보는것도 좋다.
- 1단계
from folium.plugins import HeatMap
import branca.colormap as cm
from collections import defaultdict
steps = 100
color_map= cm.LinearColormap(
colors=['blue', 'green', 'yellow', 'orange', 'red'], # heatmap 색상
vmin=0.4, vmax=0.9
).to_step(n=10)
gradient_map=defaultdict(dict)
for i in range(steps):
gradient_map[1/steps*i] = color_map.rgb_hex_str(1/steps*i)
1단계는 히트맵을 그릴 색상과 비율 범위를 선정하고, 비율에 따라 그라데이션 맵을 설정해 놓는 것이다.
vmin, vmax의 설정을 건드리면 색상을 변경할 최소 비율과 최대 비율을 변경할 수 있다.
- 2단계
# Draw a basemap
m = f.Map(location=[lat, lon], tiles='openstreetmap', zoom_start=11)
# Add heatmap
HeatMap(
data=df[['lat', 'lon']],
radius=10,
gradient=gradient_map
).add_to(m)
color_map.caption='ADD TITLE' # colorbar 문구 설정
m.add_child(color_map)
# Display the map
m2단계는 그라데이션 맵과 데이터를 기반으로 히트맵을 그리는 과정이다.
radius는 맵에 표시할 데이터들의 지점 반경이고, 이게 커질수록 히트맵이 촘촘하게(진하게?) 그려진다고 보면 된다. 이 과정을 거치면 다음과 같이 지도에 히트맵이 그려지게 될 것이다.
6. 폴리곤 그리기
원, 점이 아닌 다각형 모양의 폴리곤을 그리고 싶다면 ! 다음과 같은 코드를 사용해보자.
m = f.Map(location=[lat, lon], zoom_start=11, tiles="cartodbdark_matter")
for _, row in df.iterrows() :
sim_geo = gpd.GeoSeries(row['geometry']).simplify(tolerance=0.001)
geo_j = sim_geo.to_json()
geo_j = f.GeoJson(data=geo_j, style_function=lambda x : {'fillColor' : 'blue'})
# f.Popup(row['name']).add_to(geo_j)
geo_j.add_to(m)
m방법은 원이랑 유사해서 자세한 설명은 생략하는데, Circle 대신 GeoJson으로 goemetry 자체를 맵에 그린다고 생각하면 된다. 이 과정을 거치면 다음 그림과 같은 맵을 그릴 수 있다.
7. 코로플레스 맵 그리기
행정 구, 동, 격자 등의 단위 영역별로 수치에 따라 색상이 달라지는 코로플레스 맵을 폴리움으로 그려보자.
# 격자 단위 데이터 시각화
m = f.Map(location=[lat, lon], zoom_start=11)
c = f.Choropleth(geo_data = grid.to_json(),
data=df[['gid','value']],
columns=['gid','value'],
fill_color='YlOrRd',
key_on='properties.gid', # 격자번호 기준
highlight=True,
fill_opacity=0.8, line_opacity=0.3,
legend_name='value'
).add_to(m)
c.geojson.add_child(f.features.GeoJsonTooltip(['gid'])) # 구역 이름 표시
m나는 격자단위로 값을 표현했기 때문에 gid와, 표현할 값 컬럼(value)을 설정해서 그려줬다. 색 채우는 방식은 포맷이 여러개 있는데, 포맷을 변경하고 싶다면 다음 링크의 표현법을 참고하길 바란다.
이 과정을 거치면 다음과 같은 예쁜(?) 코로플레스 맵을 그릴 수 있다.
(이거 격자가 작아서 잘 안보이는데.. 확대하면 코로플레스가 맞아요..)
공모전 한번 치룰 때마다(몇번 안한거 같지만...?) 이런거 어떻게 했더라~~ 하고 한참 전에 돌려본 코드 뒤적뒤적 찾아보는 내가 너무 답답해서 정리했는데 이젠 그만 뒤적일 수 있을것 같다..(뭉클)
지리데이터 분석이 뽀개졌는지 정리하는 내 머리가 뽀개졌는지는 모를 일이지만,, 정리할 게 더 생기면 추가로 적어봐야지 🤔
