본문 바로가기
데이터 분석

[Python]500M*500M 정사각형으로 서울 구역 나누어 시각화에 이용하기

HanYaung 2022. 6. 29. 13:18

참고 블로그 :

 

[Python] 좌표변환 예제

UTM-K좌표계에서 WGS84좌표계로 변환하거나 WGS84좌표계에서 UTM-K좌표계로 변환해주는 예제...

blog.naver.com

https://chang12.github.io/python-map-grid/

 

Python 으로 지도 위에 정사각 격자를 그리려면?

Python 으로 지도 위에 정사각 격자를 그리려고 합니다. 정확히는 서울시를 1km x 1km 크기의 정사각 격자로 나누려고 합니다.

chang12.github.io

개요

서울 지역을 행정동보다 더 작게 구역을 나누어서 시각화에 이용해보려고 합니다. 

우선 목표는 구역을 나누고 그 구역에 포함되는 아파트들의 총 세대 수를 구해 지도에 시각화하는게 목표입니다.

 

결과 미리보기 

사용 라이브러리 

  1. pyproj  - 도표 변환과 측지 연산을 수행하는 파이썬 패키지
  2. pandas - dataframe을 사용하기 위함
  3.  folium - 웹기반 지도 라이브러리 
  4. geopandas - 지리정보 데이터 처리의 기하하적 연산과 시각화 등을 돕는 패키지

 

사용 데이터  

 1. 서울시 공동주택 아파트 정보.csv  https://data.seoul.go.kr/dataList/OA-15818/S/1/datasetView.do

 아파트 좌표 및 세대 수 읽어오기 

apart = pd.read_csv('/content/drive/MyDrive/fsi금융데이터/서울시 공동주택 아파트 정보.csv',encoding = "cp949") # 해당 엑셀 파일 불러오기
apart = apart[["k-아파트코드","좌표X","좌표Y","k-전체세대수"]].dropna().reset_index()

 

구간을 나누기

import itertools
#좌표 전환 https://m.blog.naver.com/wideeyed/221243506770  https://chang12.github.io/python-map-grid/
import folium
from pyproj import Transformer, transform ,Proj
proj_UTMK = Proj(init='epsg:5178') # UTM-K(Bassel) 도로명주소 지도 사용 중

# WGS1984
proj_WGS84 = Proj(init='epsg:4326') # Wgs84 경도/위도, GPS사용 전지구 좌표
transformer_5178_4326 = Transformer.from_proj(proj_UTMK, proj_WGS84)

p_4326 = (126.7835993001915,37.426168024755576) # 구글 지도 보고 적당히 고른 서남단 좌표.
meters = 500
n = 72 # 50x50

p_x_5178, p_y_5178 = transform(proj_WGS84,proj_UTMK , p_4326[0], p_4326[1])

#geo_jsons 생성 50*50 사이즈 정사각형 좌표로 서울을 나눔 
geo_jsons = {"type": "FeatureCollection","features":[]}
dic = {} # 폴리건을 만들기 위한 dict "n/n":[(좌표),(좌표),(좌표),(좌표)] 구조 
for x, y in itertools.product(range(n), range(n)):
    grid_lower_l = transformer_5178_4326.transform(p_x_5178 + (x + 0) * meters, p_y_5178 + (y + 0) * meters)
    grid_lower_r = transformer_5178_4326.transform(p_x_5178 + (x + 1) * meters, p_y_5178 + (y + 0) * meters)
    grid_upper_r = transformer_5178_4326.transform(p_x_5178 + (x + 1) * meters, p_y_5178 + (y + 1) * meters)
    grid_upper_l = transformer_5178_4326.transform(p_x_5178 + (x + 0) * meters, p_y_5178 + (y + 1) * meters)
    dic[str(x)+"/"+str(y)] = [(grid_lower_l[1],grid_lower_l[0]),(grid_lower_r[1],grid_lower_r[0]),(grid_upper_r[1],grid_upper_r[0]),(grid_upper_l[1],grid_upper_l[0])]

    geo_json =  {
                "type": "Feature",
                "properties":{"code": str(x)+"/"+str(y)},
                "geometry": {
                    "type": "Polygon",
                    "coordinates": [
                        [
                            grid_lower_l,
                            grid_lower_r,
                            grid_upper_r,
                            grid_upper_l,
                            grid_lower_l
                        ]
                    ]
                }
            }
    geo_jsons["features"].append(geo_json)

좌표를 500m씩 더해주기 위해서 EPSG:4326 (= 위도/경도) 좌표를 변환해야합니다 

epsg:5178로 변환을 해서  x,y좌표에 500씩 더해서 각 구간의 정사각형 구간의 꼭짓점을 구해줍니다. 

그리고 각 꼭짓점을 이용해서 polygon을 만들기 위해 dict에 꼭짓점들을 저장해줍니다.

 

구간을 polygon으로 만들고 각 polygon에 속하는 아파트 단지 찾기 

#https://incastle-study.tistory.com/3 폴리건을 사용해서 이안에 좌표가 포함되는지 확인 
from shapely.geometry import Point, Polygon
# #예시
# engineer = [(37.245452, 127.079964), (37.247041, 127.079842), (37.247912, 127.081486),(37.245588, 127.081612)]
# engineer_poly = Polygon(engineer)
# test_code1 = Point(37.245434, 127.078640)
# test_code1.within(engineer_poly)
# dic에 정사각형 폴리건들을 저장 
polygonDic ={}
for key in dic.keys():
  engineer=dic[key]
  engineer_poly = Polygon(engineer)
  polygonDic[key] = engineer_poly
#polygonDic에 저장된 폴리건 중에 포함되는 걸 찾기 
def pl(df):
  test_code1 = Point(df["좌표Y"], df["좌표X"])
  for key in polygonDic.keys():
    engineer_poly=polygonDic[key]
    if test_code1.within(engineer_poly) == True :
      return key
apart["polygon"] = apart.apply(pl,axis=1)
#나누어둔 구역에 매칭할 구역별 세대수  "x/y"로 매칭 (0~50/0~50) 
apartData = apart.groupby("polygon").sum()["k-전체세대수"]
apartData

아파트 단지마다 속한 polygon 저장해주고 이름 grouby("polygon").sum()해서 구역별 총 세대수를 구합니다. 

구간별 세대수 지도에 시각화하기 

m = folium.Map(location=(37.584930062517, 126.980289697036), zoom_start=11) # folium.Map 은 (lat, lng) 로 location 을 받음에 유의!
#구역 나눈거 적용 
folium.GeoJson(
    geo_jsons,
    name='세대수'
).add_to(m)

#구역별 아파트 세대수를 시각화하기
m.choropleth(geo_data=geo_jsons,
             data=apartData, 
             fill_color='PuRd', # 색상 변경도 가능하다
             fill_opacity=0.5,
             line_opacity=0.2,
             key_on='properties.code',
             legend_name="구역별 아파트 세대 수"
            )
#아파트 위치 circle로 찍기 
for i in apart[["좌표Y","좌표X"]].values:
  folium.Circle(location=[i[0],i[1]],radius=1).add_to(m)

m.save('index.html')
m

결과 

배운점 

1. 위치 정보를 갖는 geojson을 만드는 방법을 배움 

2. geojson을 polygon으로 만들어서 각 구역(polygon)별로 포함되는 좌표들을 구할 수 있었음

 

전체 코드 

!pip install pyproj
!pip install geopandas 
!pip3 install --upgrade setuptools
!pip3 install shapely
from google.colab import drive
drive.mount('/content/drive')
import pandas as pd
from pyproj import Transformer
import pyproj
import numpy as np
from shapely.geometry import Point as point
import geopandas as gpd
import folium
from shapely.geometry import Point, Polygon

apart = pd.read_csv('/content/drive/MyDrive/fsi금융데이터/서울시 공동주택 아파트 정보.csv',encoding = "cp949") # 해당 엑셀 파일 불러오기
apart = apart[["k-아파트코드","좌표X","좌표Y","k-전체세대수"]].dropna().reset_index()
apartCode = apart["k-아파트코드"]
apartX = apart["좌표X"]
apartY = apart["좌표Y"]
proj_UTMK = Proj(init='epsg:5178') # UTM-K(Bassel) 도로명주소 지도 사용 중

# WGS1984
proj_WGS84 = Proj(init='epsg:4326') # Wgs84 경도/위도, GPS사용 전지구 좌표
transformer_5178_4326 = Transformer.from_proj(proj_UTMK, proj_WGS84)

p_4326 = (126.7835993001915,37.426168024755576) # 구글 지도 보고 적당히 고른 서남단 좌표.
meters = 500
n = 72 # 50x50

p_x_5178, p_y_5178 = transform(proj_WGS84,proj_UTMK , p_4326[0], p_4326[1])

#geo_jsons 생성 50*50 사이즈 정사각형 좌표로 서울을 나눔 
geo_jsons = {"type": "FeatureCollection","features":[]}
dic = {} # 폴리건을 만들기 위한 dict "n/n":[(좌표),(좌표),(좌표),(좌표)] 구조 
for x, y in itertools.product(range(n), range(n)):
    grid_lower_l = transformer_5178_4326.transform(p_x_5178 + (x + 0) * meters, p_y_5178 + (y + 0) * meters)
    grid_lower_r = transformer_5178_4326.transform(p_x_5178 + (x + 1) * meters, p_y_5178 + (y + 0) * meters)
    grid_upper_r = transformer_5178_4326.transform(p_x_5178 + (x + 1) * meters, p_y_5178 + (y + 1) * meters)
    grid_upper_l = transformer_5178_4326.transform(p_x_5178 + (x + 0) * meters, p_y_5178 + (y + 1) * meters)
    dic[str(x)+"/"+str(y)] = [(grid_lower_l[1],grid_lower_l[0]),(grid_lower_r[1],grid_lower_r[0]),(grid_upper_r[1],grid_upper_r[0]),(grid_upper_l[1],grid_upper_l[0])]

    geo_json =  {
                "type": "Feature",
                "properties":{"code": str(x)+"/"+str(y)},
                "geometry": {
                    "type": "Polygon",
                    "coordinates": [
                        [
                            grid_lower_l,
                            grid_lower_r,
                            grid_upper_r,
                            grid_upper_l,
                            grid_lower_l
                        ]
                    ]
                }
            }
    geo_jsons["features"].append(geo_json)

polygonDic ={}
for key in dic.keys():
  engineer=dic[key]
  engineer_poly = Polygon(engineer)
  polygonDic[key] = engineer_poly
#polygonDic에 저장된 폴리건 중에 포함되는 걸 찾기 
def pl(df):
  test_code1 = Point(df["좌표Y"], df["좌표X"])
  for key in polygonDic.keys():
    engineer_poly=polygonDic[key]
    if test_code1.within(engineer_poly) == True :
      return key
apart["polygon"] = apart.apply(pl,axis=1)
#나누어둔 구역에 매칭할 구역별 세대수  "x/y"로 매칭 (0~50/0~50) 
apartData = apart.groupby("polygon").sum()["k-전체세대수"]
m = folium.Map(location=(37.584930062517, 126.980289697036), zoom_start=11) # folium.Map 은 (lat, lng) 로 location 을 받음에 유의!
#구역 나눈거 적용 
folium.GeoJson(
    geo_jsons,
    name='세대수'
).add_to(m)

#구역별 아파트 세대수를 시각화하기
m.choropleth(geo_data=geo_jsons,
             data=apartData, 
             fill_color='PuRd', # 색상 변경도 가능하다
             fill_opacity=0.5,
             line_opacity=0.2,
             key_on='properties.code',
             legend_name="구역별 아파트 세대 수"
            )
# #아파트 위치 circle로 찍기 
# for i in apart[["좌표Y","좌표X"]].values:
#   folium.Circle(location=[i[0],i[1]],radius=1).add_to(m)

m.save('index.html')
m

'데이터 분석' 카테고리의 다른 글

selenium 대신 playwright를 써야되는 이유  (0) 2024.03.15

'데이터 분석' Related Articles

티스토리툴바