← 메인 화면프로젝트 자료·코드
← 자료 목록

sejong_area.py

프로그램·코드 · 11.7 KB원본 다운로드

원본 내용을 읽기 전용으로 표시합니다.

  1. # 세종시 전체 분석
  3. import os
  5. import numpy as np
  6. import pandas as pd
  7. import matplotlib.pyplot as plt
  8. from sklearn.preprocessing import StandardScaler
  9. from sklearn.cluster import KMeans
  10. from sklearn.metrics import silhouette_score
  11. from sklearn.decomposition import PCA
  13. from common import *
  16. # 탐색적 데이터 분석
  17. # 업종 분포, 행정동별 지표, 행정동, 업종 구성비, 원도심/신도심
  19. def 탐색분석(now, prev, CATS):
  20. # 세종 전체의 업종 구조를 숫자와 그림으로 정리
  21. log("\n" + "=" * 70)
  22. log("2. 탐색적 데이터 분석(EDA)")
  23. log("=" * 70)
  25. # 업종 대분류별 업소 수
  26. # 세종 전체에서 어떤 업종이 많은지 먼저 확인
  27. 업종수 = now["상권업종대분류명"].value_counts()
  28. 표 = pd.DataFrame({"업소수": 업종수})
  29. 표["비율"] = (표["업소수"] / 표["업소수"].sum()).round(3)
  30. 표.index.name = "상권업종대분류명"
  31. log("\n[업종 대분류 분포]")
  32. log(표.to_string())
  34. plt.figure(figsize=(7, 4))
  35. plt.barh(업종수.index[::-1], 업종수.values[::-1], color=COLORS[0])
  36. plt.title("세종시 업종 대분류별 업소 수 (2026.6)")
  37. plt.xlabel("업소 수")
  38. 그림저장("01_업종대분류_분포.png")
  40. # 행정동별 지표
  41. # 동마다 상권 크기(업소수), 활력(신규비율), 위험(소멸률), 성격(체인, 1층, 다양성)을 한 표로 만듦
  42. dong = pd.DataFrame()
  43. dong["업소수"] = now.groupby("행정동명").size()
  44. dong["권역"] = now.groupby("행정동명")["권역"].first()
  45. dong["신규비율"] = now.groupby("행정동명")["신규"].mean()
  46. dong["체인비율"] = now.groupby("행정동명")["체인"].mean()
  47. dong = dong.sort_values("업소수", ascending=False)
  48. dong["소멸률"] = prev.groupby("행정동명")["소멸"].mean()
  49. dong["순증가율"] = now.groupby("행정동명").size() / prev.groupby("행정동명").size() - 1
  50. 층있음 = now[now["1층여부"] >= 0] # 층정보 있는 업소만
  51. dong["1층비율"] = 층있음.groupby("행정동명")["1층여부"].mean()
  52. dong["업종다양성"] = now.groupby("행정동명")["상권업종중분류명"].apply(다양성)
  53. log("\n[행정동별 주요 지표]")
  54. log(dong.round(3).to_string())
  55. dong.round(4).to_csv(os.path.join(결과폴더, "table_행정동별_지표.csv"), encoding="utf-8-sig")
  57. # 행정동별 업소 수 그래프 (원도심은 주황, 신도심은 파랑)
  58. 색 = []
  59. for g in dong["권역"]:
  60. if g == "읍면(원도심)":
  61. 색.append(COLORS[1])
  62. else:
  63. 색.append(COLORS[0])
  64. plt.figure(figsize=(7, 6))
  65. plt.barh(dong.index[::-1], dong["업소수"][::-1], color=색[::-1])
  66. plt.title("행정동별 업소 수 (주황: 읍면 원도심, 파랑: 동 신도심)")
  67. 그림저장("02_행정동별_업소수.png")
  69. # 행정동별 히트맵
  70. # 동마다 업종 비율을 색으로 비교
  71. comp = pd.crosstab(now["행정동명"], now["상권업종대분류명"], normalize="index")
  72. comp = comp[CATS]
  73. comp = comp.loc[dong.index]
  74. plt.figure(figsize=(9, 8))
  75. plt.imshow(comp.values, cmap="Blues", aspect="auto")
  76. plt.xticks(range(len(CATS)), CATS, rotation=45, ha="right")
  77. plt.yticks(range(len(comp)), comp.index)
  78. for i in range(comp.shape[0]):
  79. for j in range(comp.shape[1]):
  80. 글자색 = "black"
  81. if comp.values[i, j] > 0.25: # 진한 칸은 글자 안 보여서 흰색으로
  82. 글자색 = "white"
  83. plt.text(j, i, round(comp.values[i, j] * 100), ha="center", va="center", fontsize=7, color=글자색)
  84. plt.title("행정동별 업종 대분류 구성비(%)")
  85. plt.colorbar()
  86. 그림저장("03_행정동_업종구성_히트맵.png")
  87. comp.round(4).to_csv(os.path.join(결과폴더, "table_행정동별_업종구성비.csv"), encoding="utf-8-sig")
  89. # 원도심 vs 신도심 업종 구성
  90. # 을 두 묶음으로 합쳐서 차이를 숫자로 확인
  91. grp = pd.crosstab(now["상권업종대분류명"], now["권역"], normalize="columns").loc[CATS]
  92. log("\n[권역별 업종 대분류 구성비]")
  93. log(grp.round(3).to_string())
  95. # 행정동별 소멸률, 신규비율 산점도
  96. # x = 얼마나 없어지나, y = 얼마나 새로 생기나
  97. plt.figure(figsize=(7, 5))
  98. 신도심 = dong[dong["권역"] == "동(신도심)"]
  99. 원도심 = dong[dong["권역"] == "읍면(원도심)"]
  100. plt.scatter(신도심["소멸률"], 신도심["신규비율"], s=np.sqrt(신도심["업소수"]) * 6, color=COLORS[0], label="동(신도심)")
  101. plt.scatter(원도심["소멸률"], 원도심["신규비율"], s=np.sqrt(원도심["업소수"]) * 6, color=COLORS[1], label="읍면(원도심)")
  102. for 이름 in dong.index:
  103. if dong.loc[이름, "업소수"] > 400: # 큰 동만 이름 표시
  104. plt.text(dong.loc[이름, "소멸률"], dong.loc[이름, "신규비율"], 이름, fontsize=8)
  105. plt.xlabel("1년 소멸률 (2025.6 업소 중 사라진 비율)")
  106. plt.ylabel("신규 비율 (2026.6 업소 중 새로 생긴 비율)")
  107. plt.title("행정동별 소멸률과 신규 비율 (원 크기 = 업소 수)")
  108. plt.legend()
  109. 그림저장("04_행정동_소멸률_신규비율.png")
  111. return dong, comp
  114. # 머신러닝 - K-Means 로 상권 유형 나누기
  115. # 세종 전체를 250m 격자로 나누고, 격자마다 업종 비율/업소 수/1층 비율/체인 비율을 구해서
  116. # 비슷한 격자끼리 묶기
  118. def 상권유형군집(now, prev, CATS):
  119. # 250m 격자 단위로 K-Means 군집화해서 상권 유형을 자동으로 나눔
  120. # 행정동은 너무 커서 한 동 안에 학원가, 먹자골목, 주택가가 섞임 -> 격자로 더 잘게 봄
  121. log("\n" + "=" * 70)
  122. log("3. 머신러닝 ①: 격자(250m) 단위 상권 유형 군집화 (K-Means)")
  123. log("=" * 70)
  125. # 격자 만들기
  126. now["gx"] = (now["x"] // GRID_M).astype(int)
  127. now["gy"] = (now["y"] // GRID_M).astype(int)
  128. prev["gx"] = (prev["x"] // GRID_M).astype(int)
  129. prev["gy"] = (prev["y"] // GRID_M).astype(int)
  131. # 격자별로 값 계산
  132. g = now.groupby(["gx", "gy"])
  133. cell = pd.DataFrame()
  134. cell["업소수"] = g.size()
  135. cell["체인비율"] = g["체인"].mean()
  136. cell["신규비율"] = g["신규"].mean()
  137. cell["위도"] = g["위도"].mean()
  138. cell["경도"] = g["경도"].mean()
  139. cell["대표동"] = g["행정동명"].agg(최빈값)
  140. 층있음 = now[now["1층여부"] >= 0]
  141. cell["1층비율"] = 층있음.groupby(["gx", "gy"])["1층여부"].mean()
  142. cell["1층비율"] = cell["1층비율"].fillna(cell["1층비율"].median()) # 층정보 없는 격자는 중앙값으로
  143. 업종비율 = pd.crosstab([now["gx"], now["gy"]], now["상권업종대분류명"], normalize="index")
  144. cell = cell.join(업종비율[CATS])
  145. cell = cell[cell["업소수"] >= MIN_STORES].copy() # 업소 적은 격자 빼기
  146. cell["log업소수"] = np.log1p(cell["업소수"]) # 업소 수 차이가 너무 커서 로그 씌움
  147. log(f"분석 대상 격자: {len(cell)}개 (업소 {MIN_STORES}개 이상), 포함 업소 {cell['업소수'].sum():,}개")
  149. # K-Means
  150. # 변수 단위가 다르니까 표준화한 뒤 군집
  151. # 군집 수 k는 엘보 + 실루엣으로 봄. 실루엣이 전체적으로 낮아서(0.1대) 경계가 뚜렷하진 않음
  152. FEATS = CATS + ["log업소수", "1층비율", "체인비율"]
  153. X = StandardScaler().fit_transform(cell[FEATS]) # 표준화 (단위 맞추기)
  155. # k 정하기 : k를 3~8까지 해보고 점수 확인
  156. ks = [3, 4, 5, 6, 7, 8]
  157. inertias = []
  158. sils = []
  159. for k in ks:
  160. km = KMeans(n_clusters=k, n_init=20, random_state=SEED)
  161. km.fit(X)
  162. inertias.append(km.inertia_)
  163. sils.append(silhouette_score(X, km.labels_))
  164. log("k별 실루엣 점수: " + ", ".join(f"k={k}:{s:.3f}" for k, s in zip(ks, sils)))
  166. # k=3은 유형이 너무 적어서 4 이상 중에 실루엣 제일 높은 걸로
  167. BEST_K = 4
  168. 최고점수 = -1
  169. for k, s in zip(ks, sils):
  170. if k >= 4 and s > 최고점수:
  171. 최고점수 = s
  172. BEST_K = k
  173. log(f"선택된 군집 수 K = {BEST_K} (k>=4 중 실루엣 최대)")
  175. plt.figure(figsize=(10, 3.6))
  176. plt.subplot(1, 2, 1)
  177. plt.plot(ks, inertias, marker="o")
  178. plt.title("엘보 곡선")
  179. plt.xlabel("군집 수 k")
  180. plt.subplot(1, 2, 2)
  181. plt.plot(ks, sils, marker="o")
  182. plt.axvline(BEST_K, color="gray", linestyle=":")
  183. plt.title("실루엣 점수")
  184. plt.xlabel("군집 수 k")
  185. 그림저장("05_군집수_선택.png")
  187. # 최종 K-Means
  188. km = KMeans(n_clusters=BEST_K, n_init=50, random_state=SEED)
  189. km.fit(X)
  190. cell["군집"] = km.labels_
  192. # 군집별 평균
  193. prof = cell.groupby("군집")[CATS + ["업소수", "1층비율", "체인비율", "신규비율"]].mean()
  194. prof["격자수"] = cell.groupby("군집").size()
  195. 전체평균 = cell[CATS].mean()
  197. # 군집 이름 붙이기
  198. # 손으로 붙이면 k나 데이터가 바뀔 때마다 고쳐야 하니까 규칙으로 함
  199. 상위25 = cell["업소수"].quantile(0.75)
  200. 중앙값 = cell["업소수"].quantile(0.5)
  201. 이름들 = []
  202. for c in prof.index:
  203. r = prof.loc[c]
  204. if r["업소수"] >= 상위25:
  205. 밀도 = "고밀"
  206. elif r["업소수"] <= 중앙값:
  207. 밀도 = "저밀"
  208. else:
  209. 밀도 = "중밀"
  210. 배율 = (r[CATS] / 전체평균).sort_values(ascending=False)
  211. if 배율.iloc[0] < 1.3:
  212. 이름 = 밀도 + " 업종혼합"
  213. elif 배율.iloc[1] >= 1.3:
  214. 이름 = f"{밀도} {배율.index[0]}·{배율.index[1]} 특화"
  215. else:
  216. 이름 = f"{밀도} {배율.index[0]} 특화"
  217. 이름들.append(이름)
  218. prof["유형명"] = 이름들
  219. # 같은 이름이 있으면 뒤에 번호 붙이기
  220. 겹침 = prof["유형명"].duplicated(keep=False)
  221. prof.loc[겹침, "유형명"] = prof.loc[겹침, "유형명"] + " " + prof.loc[겹침].index.astype(str)
  222. cell["유형명"] = cell["군집"].map(prof["유형명"])
  223. log("\n[군집(상권 유형) 프로파일]")
  224. log(prof.round(3).to_string())
  225. prof.round(4).to_csv(os.path.join(결과폴더, "table_상권유형_군집프로파일.csv"), encoding="utf-8-sig")
  227. # 동별로 어떤 유형 격자가 몇 개 있는지
  228. dong_type = pd.crosstab(cell["대표동"], cell["유형명"])
  229. log("\n[행정동별 상권 유형 격자 수]")
  230. log(dong_type.to_string())
  232. # 결과 그림
  233. # 2차원으로 줄여서 군집 그림 그리기
  234. pca = PCA(n_components=2, random_state=SEED)
  235. P = pca.fit_transform(X)
  236. plt.figure(figsize=(7, 5))
  237. for c in range(BEST_K):
  238. m = cell["군집"].values == c
  239. plt.scatter(P[m, 0], P[m, 1], s=20, color=COLORS[c % 8], label=prof.loc[c, "유형명"])
  240. plt.title(f"상권 유형 군집 결과 (PCA 2차원, 설명력 {pca.explained_variance_ratio_.sum():.0%})")
  241. plt.legend(fontsize=8)
  242. 그림저장("06_군집_PCA.png")
  244. # 군집별 업종 배율 히트맵 (1보다 크면 평균보다 많음)
  245. lift = prof[CATS] / 전체평균
  246. plt.figure(figsize=(9, 0.6 * BEST_K + 1.5))
  247. plt.imshow(lift.values, cmap="Blues", aspect="auto", vmin=0, vmax=3)
  248. plt.xticks(range(len(CATS)), CATS, rotation=45, ha="right")
  249. plt.yticks(range(BEST_K), prof["유형명"])
  250. for i in range(lift.shape[0]):
  251. for j in range(lift.shape[1]):
  252. 글자색 = "black"
  253. if lift.values[i, j] > 1.8:
  254. 글자색 = "white"
  255. plt.text(j, i, round(lift.values[i, j], 1), ha="center", va="center", fontsize=8, color=글자색)
  256. plt.title("상권 유형별 업종 특화 배율 (전체 평균=1.0)")
  257. plt.colorbar()
  258. 그림저장("07_상권유형_업종특화배율.png")
  260. return cell, prof, BEST_K
파일 확인 정보

SHA-256: f32ac91417bfa4440e8c7ac152477547bb5c0e79f22cf41f6c51460ab3dc7da8