{ "cells": [ { "cell_type": "markdown", "id": "9afdef8a", "metadata": {}, "source": [ "# 경계정보 생성" ] }, { "cell_type": "code", "execution_count": 18, "id": "014473d8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "🔍 인하대학교 경계 정보 추출 중...\n", "응답 상태 코드: 504\n", "❌ API 오류 발생\n", "🔄 대체 방법을 시도합니다...\n", "📍 Plan B: 좌표 범위 검색 방법...\n", "✓ 대체 방법 성공!\n", "✓ 3개의 요소 발견\n", "🔍 받은 요소 개수: 3\n", " [0] type: way, tags: {'amenity': 'university', 'building': 'yes', 'education': 'university', 'name': '나빌레관', 'name:en': 'Student Club Hall'}\n", " 좌표 개수: 13\n", " 면적: 0.00000008\n", " ✓ 추가됨\n", " [1] type: way, tags: {'amenity': 'university', 'education': 'university', 'name': '인하대학교', 'name:en': 'Inha University', 'name:ja': '仁荷大学校', 'name:zh': '仁荷大学', 'wikidata': 'Q488851'}\n", " 좌표 개수: 37\n", " 면적: 0.00002906\n", " ✓ 추가됨\n", " [2] type: way, tags: {'amenity': 'university', 'education': 'university', 'name': '인하공업전문대학', 'name:en': 'Inha Technical College', 'name:ja': '仁荷工業専門大学', 'name:zh': '仁荷工业专门大学'}\n", " 좌표 개수: 26\n", " 면적: 0.00001058\n", " ✓ 추가됨\n", "📊 총 3개 feature 생성\n", " 가장 큰 polygon 선택 (면적: 0.00002906)\n", "✓ GeoJSON 파일 저장 완료: inha_boundary.geojson\n", "✓ 위성 이미지 (Esri) 추가 완료\n", "✓ 위성 이미지 저장 완료: inha_boundary_satellite.png\n", "✓ 지도 스타일 + 클리핑 추가 완료\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\USER\\AppData\\Local\\Temp\\ipykernel_4744\\1166856351.py:227: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.\n", " poly = gdf_web.geometry.unary_union\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "✓ 지도 이미지 저장 완료: inha_boundary_map.png\n", "✓ 폴리곤 이미지 저장 완료: inha_boundary_polygon.png\n", "✓ Shapefile 저장 완료: inha_boundary.shp\n", "\n", "📍 경계 정보:\n", " - 서쪽: 126.649007, 남쪽: 37.446685\n", " - 동쪽: 126.657847, 북쪽: 37.452137\n", " - 중심점: (126.653435, 37.449753)\n", "\n", "📌 첫 번째 경계의 좌표:\n", " 총 37개의 좌표점\n", " 시작점: (126.6504913, 37.4488536)\n", " 끝점: (126.6504913, 37.4488536)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "c:\\Users\\USER\\anaconda3\\envs\\torch_ev\\lib\\site-packages\\pyogrio\\raw.py:723: RuntimeWarning: Normalized/laundered field name: 'name:en' to 'name_en'\n", " ogr_write(\n", "c:\\Users\\USER\\anaconda3\\envs\\torch_ev\\lib\\site-packages\\pyogrio\\raw.py:723: RuntimeWarning: Normalized/laundered field name: 'name:ja' to 'name_ja'\n", " ogr_write(\n", "c:\\Users\\USER\\anaconda3\\envs\\torch_ev\\lib\\site-packages\\pyogrio\\raw.py:723: RuntimeWarning: Normalized/laundered field name: 'name:zh' to 'name_zh'\n", " ogr_write(\n" ] } ], "source": [ "import numpy as np\n", "import matplotlib.path as mpath\n", "import matplotlib.patches as mpatches\n", "import requests\n", "import json\n", "import geopandas as gpd\n", "from shapely.geometry import shape, Polygon, MultiPolygon\n", "import matplotlib.pyplot as plt\n", "import contextily as ctx\n", "from matplotlib import font_manager\n", "\n", "# 대체 방법: 좌표 범위로 검색\n", "def get_inha_boundary_alternative():\n", " print(\"📍 Plan B: 좌표 범위 검색 방법...\")\n", " overpass_url = \"http://overpass-api.de/api/interpreter\"\n", " \n", " # 인하대학교 주변 영역 (위도, 경도 범위)\n", " overpass_query = f\"\"\"\n", " [out:json][timeout:25];\n", " (\n", " way[\"amenity\"=\"university\"](37.445,126.650,37.456,126.662);\n", " relation[\"amenity\"=\"university\"](37.445,126.650,37.456,126.662);\n", " way[\"name~\"=\"인하\"](37.445,126.650,37.456,126.662);\n", " );\n", " out geom;\n", " \"\"\"\n", " \n", " try:\n", " response = requests.post(overpass_url, data={'data': overpass_query}, timeout=30)\n", " if response.status_code == 200:\n", " data = response.json()\n", " if data.get('elements'):\n", " print(\"✓ 대체 방법 성공!\")\n", " return data\n", " except Exception as e:\n", " print(f\"⚠️ 대체 방법 실패: {e}\")\n", " \n", " print(\"📍 Plan C: 수동 경계 생성...\")\n", " return create_manual_boundary()\n", "\n", "# 수동으로 경계 생성 (실제 인하대 메인 캠퍼스만)\n", "def create_manual_boundary():\n", " print(\"✅ 인하대학교 메인 캠퍼스 경계 생성\")\n", " # 인하대학교 메인 캠퍼스 경계 (정문-5호관-60주년기념관-본관-체육관-후문 순)\n", " coords = [\n", " [126.65168, 37.45100], # 정문\n", " [126.65150, 37.45095],\n", " [126.65140, 37.45020],\n", " [126.65100, 37.44920],\n", " [126.65090, 37.44880],\n", " [126.65120, 37.44860],\n", " [126.65280, 37.44780],\n", " [126.65450, 37.44680],\n", " [126.65580, 37.44680], # 남동쪽 모서리\n", " [126.65750, 37.44720],\n", " [126.65850, 37.44780],\n", " [126.65900, 37.44860],\n", " [126.65920, 37.44980],\n", " [126.65900, 37.45100], # 동쪽 경계\n", " [126.65850, 37.45180],\n", " [126.65750, 37.45240],\n", " [126.65620, 37.45260], # 북동쪽 모서리\n", " [126.65500, 37.45250],\n", " [126.65380, 37.45220],\n", " [126.65250, 37.45180],\n", " [126.65168, 37.45100], # 정문으로 닫기\n", " ]\n", " \n", " return {\n", " 'elements': [{\n", " 'type': 'way',\n", " 'geometry': [{'lon': c[0], 'lat': c[1]} for c in coords],\n", " 'tags': {'name': '인하대학교', 'amenity': 'university'}\n", " }]\n", " }\n", "\n", "# OpenStreetMap Overpass API를 사용하여 인하대학교 경계 추출\n", "def get_inha_boundary():\n", " overpass_url = \"http://overpass-api.de/api/interpreter\"\n", " \n", " # 인하대학교를 검색하는 Overpass 쿼리 (더 정확한 필터링)\n", " overpass_query = \"\"\"\n", " [out:json][timeout:25];\n", " (\n", " way[\"name\"=\"인하대학교\"][\"amenity\"=\"university\"];\n", " relation[\"name\"=\"인하대학교\"][\"amenity\"=\"university\"];\n", " );\n", " out geom;\n", " \"\"\"\n", " \n", " try:\n", " response = requests.get(overpass_url, params={'data': overpass_query}, timeout=30)\n", " print(f\"응답 상태 코드: {response.status_code}\")\n", " \n", " if response.status_code != 200:\n", " print(f\"❌ API 오류 발생\")\n", " print(\"🔄 대체 방법을 시도합니다...\")\n", " return get_inha_boundary_alternative()\n", " \n", " data = response.json()\n", " if not data.get('elements'):\n", " print(\"⚠️ 검색 결과가 없습니다. 대체 방법을 시도합니다...\")\n", " return get_inha_boundary_alternative()\n", " return data\n", " except Exception as e:\n", " print(f\"❌ 요청 실패: {e}\")\n", " print(\"🔄 대체 방법을 시도합니다...\")\n", " return get_inha_boundary_alternative()\n", "\n", "# GeoJSON 형식으로 변환\n", "def convert_to_geojson(osm_data):\n", " features = []\n", " \n", " print(f\"🔍 받은 요소 개수: {len(osm_data.get('elements', []))}\")\n", " \n", " for idx, element in enumerate(osm_data.get('elements', [])):\n", " tags = element.get('tags', {})\n", " print(f\" [{idx}] type: {element['type']}, tags: {tags}\")\n", " \n", " if element['type'] == 'way' and 'geometry' in element:\n", " coords = [(node['lon'], node['lat']) for node in element['geometry']]\n", " print(f\" 좌표 개수: {len(coords)}\")\n", " \n", " if len(coords) > 5: # 최소 5개 이상\n", " try:\n", " polygon = Polygon(coords)\n", " area = polygon.area\n", " print(f\" 면적: {area:.8f}\")\n", " \n", " if polygon.is_valid:\n", " features.append({\n", " 'type': 'Feature',\n", " 'geometry': polygon.__geo_interface__,\n", " 'properties': tags\n", " })\n", " print(f\" ✓ 추가됨\")\n", " except Exception as e:\n", " print(f\" ✗ Polygon 생성 실패: {e}\")\n", " \n", " elif element['type'] == 'relation' and 'members' in element:\n", " all_coords = []\n", " for member in element['members']:\n", " if 'geometry' in member:\n", " coords = [(node['lon'], node['lat']) for node in member['geometry']]\n", " all_coords.extend(coords)\n", " \n", " print(f\" relation 좌표 개수: {len(all_coords)}\")\n", " \n", " if len(all_coords) > 5:\n", " try:\n", " polygon = Polygon(all_coords)\n", " if polygon.is_valid:\n", " features.append({\n", " 'type': 'Feature',\n", " 'geometry': polygon.__geo_interface__,\n", " 'properties': tags\n", " })\n", " print(f\" ✓ relation 추가됨\")\n", " except Exception as e:\n", " print(f\" ✗ Relation polygon 실패: {e}\")\n", " \n", " print(f\"📊 총 {len(features)}개 feature 생성\")\n", " \n", " # 면적이 가장 큰 polygon 선택\n", " if features and len(features) > 1:\n", " features.sort(key=lambda f: Polygon(f['geometry']['coordinates'][0]).area, reverse=True)\n", " print(f\" 가장 큰 polygon 선택 (면적: {Polygon(features[0]['geometry']['coordinates'][0]).area:.8f})\")\n", " features = [features[0]]\n", " \n", " geojson = {\n", " 'type': 'FeatureCollection',\n", " 'features': features\n", " }\n", " \n", " return geojson\n", "\n", "# 경계 시각화 및 저장\n", "def visualize_and_save(geojson_data, output_file='inha_boundary.geojson'):\n", " # GeoJSON 파일로 저장\n", " with open(output_file, 'w', encoding='utf-8') as f:\n", " json.dump(geojson_data, f, ensure_ascii=False, indent=2)\n", " print(f\"✓ GeoJSON 파일 저장 완료: {output_file}\")\n", " \n", " # GeoDataFrame으로 변환하여 시각화\n", " if geojson_data['features']:\n", " gdf = gpd.GeoDataFrame.from_features(geojson_data['features'])\n", " gdf.set_crs(epsg=4326, inplace=True)\n", " \n", " # Web Mercator 투영법으로 변환\n", " gdf_web = gdf.to_crs(epsg=3857)\n", " \n", " # 한글 폰트 설정\n", " try:\n", " font_path = font_manager.findfont(font_manager.FontProperties(family='Malgun Gothic'))\n", " plt.rcParams['font.family'] = 'Malgun Gothic'\n", " except:\n", " plt.rcParams['font.family'] = 'sans-serif'\n", " \n", " # 방법 1: 위성 이미지 (Esri WorldImagery - 최신)\n", " fig1, ax1 = plt.subplots(figsize=(15, 12))\n", " gdf_web.plot(ax=ax1, facecolor='none', edgecolor='red', linewidth=3, alpha=0.9)\n", " \n", " try:\n", " # 최신 위성 이미지 소스 시도\n", " ctx.add_basemap(ax1, source=ctx.providers.Esri.WorldImagery, alpha=0.8)\n", " print(\"✓ 위성 이미지 (Esri) 추가 완료\")\n", " except:\n", " try:\n", " # Google Satellite 대체\n", " ctx.add_basemap(ax1, source='https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}', alpha=0.8)\n", " print(\"✓ 위성 이미지 (Google) 추가 완료\")\n", " except Exception as e:\n", " print(f\"⚠️ 위성 이미지 추가 실패: {e}\")\n", " \n", " ax1.set_title('인하대학교 경계 (위성)', fontsize=18, pad=20, weight='bold')\n", " ax1.axis('off')\n", " plt.tight_layout()\n", " plt.savefig('inha_boundary_satellite.png', dpi=300, bbox_inches='tight', facecolor='white')\n", " print(\"✓ 위성 이미지 저장 완료: inha_boundary_satellite.png\")\n", " plt.close(fig1)\n", " \n", " # 방법 2: 지도 스타일 (OpenStreetMap)\n", " # 방법 2: 지도 스타일 (OpenStreetMap) + 경계 밖 배경 클리핑\n", " fig2, ax2 = plt.subplots(figsize=(15, 12))\n", "\n", " # Web Mercator 기준 폴리곤 준비 (합쳐진 하나의 geometry)\n", " poly = gdf_web.geometry.unary_union\n", "\n", " # MultiPolygon이면 가장 큰 것 하나만 사용 (선택)\n", " if poly.geom_type == \"MultiPolygon\":\n", " poly = max(poly.geoms, key=lambda g: g.area)\n", "\n", " # 폴리곤을 matplotlib PathPatch로 변환 (clip용)\n", " x, y = poly.exterior.coords.xy\n", " coords = np.vstack([x, y]).T\n", "\n", " Path = mpath.Path\n", " path_data = [(Path.MOVETO, coords[0])]\n", " for c in coords[1:]:\n", " path_data.append((Path.LINETO, c))\n", " path_data.append((Path.CLOSEPOLY, coords[0]))\n", "\n", " codes, verts = zip(*path_data)\n", " clip_path = mpath.Path(verts, codes)\n", " clip_patch = mpatches.PathPatch(clip_path, facecolor=\"none\", edgecolor=\"none\")\n", "\n", " ax2.add_patch(clip_patch)\n", "\n", " # 우선 폴리곤을 색칠\n", " gdf_web.plot(ax=ax2, facecolor='lightblue', edgecolor='darkblue', linewidth=2.5, alpha=0.6)\n", "\n", " try:\n", " ctx.add_basemap(ax2, source=ctx.providers.OpenStreetMap.Mapnik, alpha=0.7)\n", " # 방금 추가된 배경 이미지에 클립 적용\n", " if ax2.images:\n", " ax2.images[-1].set_clip_path(clip_patch)\n", " print(\"✓ 지도 스타일 + 클리핑 추가 완료\")\n", " except Exception as e:\n", " print(f\"⚠️ 지도 추가 실패: {e}\")\n", " \n", " ax2.set_title('인하대학교 경계 (지도, boundary 밖 배경 제거)', fontsize=18, pad=20, weight='bold')\n", " ax2.axis('off')\n", " plt.tight_layout()\n", " plt.savefig('inha_boundary_map.png', dpi=300, bbox_inches='tight', facecolor='white')\n", " print(\"✓ 지도 이미지 저장 완료: inha_boundary_map.png\")\n", " plt.close(fig2)\n", "\n", " \n", " # 방법 3: 경계만 (폴리곤만 강조)\n", " fig3, ax3 = plt.subplots(figsize=(12, 10))\n", " gdf.plot(ax=ax3, facecolor='#4CAF50', edgecolor='#1B5E20', linewidth=2.5, alpha=0.4)\n", " ax3.set_title('인하대학교 캠퍼스 영역', fontsize=18, pad=20, weight='bold')\n", " ax3.set_xlabel('경도 (Longitude)', fontsize=12)\n", " ax3.set_ylabel('위도 (Latitude)', fontsize=12)\n", " ax3.grid(True, alpha=0.3, linestyle='--', linewidth=0.5)\n", " ax3.set_facecolor('#f5f5f5')\n", " plt.tight_layout()\n", " plt.savefig('inha_boundary_polygon.png', dpi=300, bbox_inches='tight', facecolor='white')\n", " print(\"✓ 폴리곤 이미지 저장 완료: inha_boundary_polygon.png\")\n", " plt.close(fig3)\n", " \n", " # Shapefile로도 저장 (GIS 소프트웨어에서 사용 가능)\n", " try:\n", " gdf.to_file('inha_boundary.shp')\n", " print(\"✓ Shapefile 저장 완료: inha_boundary.shp\")\n", " except Exception as e:\n", " print(f\"⚠️ Shapefile 저장 실패: {e}\")\n", " \n", " # 경계 정보 출력\n", " print(f\"\\n📍 경계 정보:\")\n", " bounds = gdf.total_bounds\n", " print(f\" - 서쪽: {bounds[0]:.6f}, 남쪽: {bounds[1]:.6f}\")\n", " print(f\" - 동쪽: {bounds[2]:.6f}, 북쪽: {bounds[3]:.6f}\")\n", " \n", " # CRS 변환 후 중심점 계산 (경고 방지)\n", " gdf_projected = gdf.to_crs(epsg=5179) # 한국 평면직각좌표계\n", " centroid_proj = gdf_projected.centroid.iloc[0]\n", " centroid_wgs84 = gpd.GeoSeries([centroid_proj], crs=5179).to_crs(4326).iloc[0]\n", " print(f\" - 중심점: ({centroid_wgs84.x:.6f}, {centroid_wgs84.y:.6f})\")\n", " \n", " return gdf\n", " else:\n", " print(\"❌ 경계 데이터를 찾을 수 없습니다.\")\n", " return None\n", "\n", "# 메인 실행\n", "if __name__ == \"__main__\":\n", " print(\"🔍 인하대학교 경계 정보 추출 중...\")\n", " \n", " # 1. OSM 데이터 가져오기\n", " osm_data = get_inha_boundary()\n", " print(f\"✓ {len(osm_data.get('elements', []))}개의 요소 발견\")\n", " \n", " # 2. GeoJSON으로 변환\n", " geojson = convert_to_geojson(osm_data)\n", " \n", " # 3. 저장 및 시각화\n", " gdf = visualize_and_save(geojson)\n", " \n", " # 4. 좌표 데이터 출력 (필요시)\n", " if geojson['features']:\n", " print(\"\\n📌 첫 번째 경계의 좌표:\")\n", " coords = geojson['features'][0]['geometry']['coordinates']\n", " print(f\" 총 {len(coords[0]) if coords else 0}개의 좌표점\")\n", " if coords and len(coords[0]) > 0:\n", " print(f\" 시작점: {coords[0][0]}\")\n", " print(f\" 끝점: {coords[0][-1]}\")" ] }, { "cell_type": "markdown", "id": "79139ad6", "metadata": {}, "source": [ "# 파일로 불러오기" ] }, { "cell_type": "code", "execution_count": 21, "id": "6f2d7463", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\USER\\AppData\\Local\\Temp\\ipykernel_4744\\2019802498.py:19: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.\n", " poly = gdf_web.geometry.unary_union\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "✓ 저장 완료: inha_boundary_map_from_shp.png\n" ] } ], "source": [ "import geopandas as gpd\n", "import matplotlib.pyplot as plt\n", "import contextily as ctx\n", "import numpy as np\n", "import matplotlib.path as mpath\n", "import matplotlib.patches as mpatches\n", "\n", "def plot_clipped_map_from_shp(\n", " shp_path=\"inha_boundary.shp\",\n", " out_png=\"inha_boundary_map_from_shp.png\"\n", "):\n", " # 1) Shapefile 읽기 (shx/dbf/prj는 같이 자동으로 읽힘)\n", " gdf = gpd.read_file(shp_path)\n", "\n", " # 2) Web Mercator로 변환 (contextily용)\n", " gdf_web = gdf.to_crs(epsg=3857)\n", "\n", " # 3) 폴리곤 하나로 합치기 (MultiPolygon이면 가장 큰 것 선택)\n", " poly = gdf_web.geometry.unary_union\n", " if poly.geom_type == \"MultiPolygon\":\n", " poly = max(poly.geoms, key=lambda g: g.area)\n", "\n", " # 4) clip에 쓸 PathPatch 만들기\n", " x, y = poly.exterior.coords.xy\n", " coords = np.vstack([x, y]).T\n", "\n", " Path = mpath.Path\n", " path_data = [(Path.MOVETO, coords[0])]\n", " for c in coords[1:]:\n", " path_data.append((Path.LINETO, c))\n", " path_data.append((Path.CLOSEPOLY, coords[0]))\n", "\n", " codes, verts = zip(*path_data)\n", " clip_path = mpath.Path(verts, codes)\n", " clip_patch = mpatches.PathPatch(clip_path, facecolor=\"none\", edgecolor=\"none\")\n", "\n", " # 5) 그림 생성\n", " fig, ax = plt.subplots(figsize=(15, 12))\n", "\n", " # 폴리곤 먼저 색칠 (경계 확인용)\n", " gdf_web.plot(ax=ax, facecolor=\"lightblue\", edgecolor=\"darkblue\", alpha=0.6, linewidth=2)\n", "\n", " # clip 패치 추가\n", " ax.add_patch(clip_patch)\n", "\n", " # 6) basemap 추가 + 클리핑\n", " ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, alpha=0.7)\n", " if ax.images:\n", " ax.images[-1].set_clip_path(clip_patch)\n", "\n", " # 7) 캠퍼스 영역만 보이게 축 범위도 폴리곤에 맞게 잘라주기\n", " minx, miny, maxx, maxy = poly.bounds\n", " ax.set_xlim(minx, maxx)\n", " ax.set_ylim(miny, maxy)\n", " ax.axis(\"off\")\n", " plt.tight_layout()\n", " plt.savefig(out_png, dpi=300, bbox_inches=\"tight\", facecolor=\"white\")\n", " plt.close(fig)\n", "\n", " print(f\"✓ 저장 완료: {out_png}\")\n", "\n", "\n", "plot_clipped_map_from_shp()\n" ] }, { "cell_type": "markdown", "id": "9ee6be64", "metadata": {}, "source": [ "# 표본 추출 비교" ] }, { "cell_type": "code", "execution_count": 29, "id": "c68604ed", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\USER\\AppData\\Local\\Temp\\ipykernel_4744\\3417582304.py:15: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.\n", " poly = gdf.geometry.unary_union # 하나의 폴리곤으로 합치기\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "격자 층별 표본 수 할당 (면적비례): [ 4 12 7 12 16 9]\n", "동일면적 층별 표본 수 할당 (동일 표본수): [10 10 10 10 10 10]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\USER\\AppData\\Local\\Temp\\ipykernel_4744\\3417582304.py:205: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.\n", " poly_web = gdf_web.geometry.unary_union\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "✓ 원본 경계 지도 저장: inha_boundary_base.png\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\USER\\AppData\\Local\\Temp\\ipykernel_4744\\3417582304.py:205: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.\n", " poly_web = gdf_web.geometry.unary_union\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "✓ SRS 샘플 지도 저장: inha_sampling_srs.png\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\USER\\AppData\\Local\\Temp\\ipykernel_4744\\3417582304.py:205: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.\n", " poly_web = gdf_web.geometry.unary_union\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "✓ 층화 샘플 지도 저장: inha_sampling_strata.png\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\USER\\AppData\\Local\\Temp\\ipykernel_4744\\3417582304.py:205: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.\n", " poly_web = gdf_web.geometry.unary_union\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "✓ 동일면적 층화 샘플 지도 저장: inha_sampling_strata_equalarea.png\n", "\n", "[SRS 샘플 좌표 (EPSG:5179)]\n", " design geometry\n", "0 SRS POINT (925278.467 1939423.964)\n", "1 SRS POINT (925060.691 1939262.149)\n", "2 SRS POINT (925344.356 1939270.6)\n", "3 SRS POINT (925259.687 1939296.418)\n", "4 SRS POINT (925386.793 1939152.443)\n", "\n", "[격자 층화 샘플 좌표 (EPSG:5179)]\n", " stratum_id geometry\n", "0 0 POINT (924882.652 1939247.19)\n", "1 0 POINT (924924.472 1939181.424)\n", "2 0 POINT (924947.85 1939212.952)\n", "3 0 POINT (924930.148 1939196.918)\n", "4 1 POINT (925147.586 1939174.715)\n", "\n", "[동일면적 층화 샘플 좌표 (EPSG:5179)]\n", " stratum_id geometry\n", "0 0 POINT (924846.854 1939231.61)\n", "1 0 POINT (924830.9 1939503.537)\n", "2 0 POINT (924832.812 1939290.495)\n", "3 0 POINT (924882.953 1939275.119)\n", "4 0 POINT (924776.295 1939433.993)\n" ] } ], "source": [ "import numpy as np\n", "import geopandas as gpd\n", "from shapely.geometry import Point, box\n", "import matplotlib.pyplot as plt\n", "import contextily as ctx\n", "import matplotlib.path as mpath\n", "import matplotlib.patches as mpatches\n", "\n", "# ===========================\n", "# 1. 공통: 경계 불러오기\n", "# ===========================\n", "def load_inha_boundary(shp_path=\"inha_boundary.shp\"):\n", " gdf = gpd.read_file(shp_path) # prj 기준으로 CRS 인식\n", " gdf = gdf.to_crs(epsg=5179) # 한국 평면직각좌표계\n", " poly = gdf.geometry.unary_union # 하나의 폴리곤으로 합치기\n", " return gdf, poly\n", "\n", "\n", "# ===========================\n", "# 2. 폴리곤 내부 랜덤 샘플링 (SRS)\n", "# ===========================\n", "def random_points_in_polygon(poly, n):\n", " minx, miny, maxx, maxy = poly.bounds\n", " points = []\n", " # rejection sampling: bounding box에서 뽑아서, 폴리곤 안만 사용\n", " while len(points) < n:\n", " xs = np.random.uniform(minx, maxx, n * 2)\n", " ys = np.random.uniform(miny, maxy, n * 2)\n", " for x, y in zip(xs, ys):\n", " p = Point(x, y)\n", " if poly.contains(p): # 경계까지 포함하려면 contains 대신 covers 써도 됨\n", " points.append(p)\n", " if len(points) == n:\n", " break\n", " return points\n", "\n", "\n", "# ===========================\n", "# 3-1. 폴리곤을 격자로 나눈 후, 층(구역) 만들기 (기존 방식)\n", "# ===========================\n", "def make_strata(poly, n_rows=2, n_cols=3):\n", " \"\"\"\n", " n_rows x n_cols 격자 생성 후,\n", " 각 cell과 캠퍼스 폴리곤의 교차(intersection)를 strata로 사용\n", " \"\"\"\n", " minx, miny, maxx, maxy = poly.bounds\n", " dx = (maxx - minx) / n_cols\n", " dy = (maxy - miny) / n_rows\n", "\n", " strata_geoms = []\n", " for i in range(n_rows):\n", " for j in range(n_cols):\n", " cell = box(\n", " minx + j * dx,\n", " miny + i * dy,\n", " minx + (j + 1) * dx,\n", " miny + (i + 1) * dy,\n", " )\n", " inter = cell.intersection(poly)\n", " if not inter.is_empty and inter.area > 0:\n", " strata_geoms.append(inter)\n", "\n", " strata_gdf = gpd.GeoDataFrame(\n", " {\"stratum_id\": range(len(strata_geoms))},\n", " geometry=strata_geoms,\n", " crs=\"EPSG:5179\"\n", " )\n", " return strata_gdf\n", "\n", "\n", "# ===========================\n", "# 3-2. 동일면적으로 나누는 층(구역) 만들기\n", "# - x축 방향으로 잘라서 각 밴드의 폴리곤 면적이 거의 같게\n", "# ===========================\n", "# ===========================\n", "# 3-2. 동일면적으로 나누는 층(구역) 만들기 (빠른 근사 버전)\n", "# - 폴리곤 안에 무작위 점 많이 뽑고\n", "# - x좌표 기준으로 6등분(quantile) → 거의 동일면적 스트립\n", "# ===========================\n", "def make_equal_area_strata(poly, n_strata=6, n_mc=10000):\n", " \"\"\"\n", " poly를 x축 방향 스트립 n_strata개로 나누되,\n", " 폴리곤 내부에 균일 샘플링한 점들의 x좌표 분포를 이용해\n", " '근사 동일면적'이 되도록 자른다.\n", " - n_mc: 몬테카를로 포인트 개수 (많을수록 면적 근사 정확, 느려짐)\n", " \"\"\"\n", "\n", " # 1) 폴리곤 내부에 몬테카를로 포인트 생성\n", " mc_points = random_points_in_polygon(poly, n_mc)\n", " xs = np.array([p.x for p in mc_points])\n", "\n", " # 2) x좌표 정렬 후, n_strata등분 quantile 위치 찾기\n", " xs_sorted = np.sort(xs)\n", " cuts = []\n", " for i in range(1, n_strata):\n", " q = i / n_strata\n", " idx = int(q * len(xs_sorted))\n", " cuts.append(xs_sorted[idx])\n", "\n", " # 3) 잘라낼 x 경계 만들기 (좌/우 끝 포함)\n", " minx, miny, maxx, maxy = poly.bounds\n", " x_edges = [minx] + cuts + [maxx]\n", "\n", " # 4) 각 스트립과 폴리곤의 교차로 strata 생성 (intersection은 스트립 개수만큼 딱 6번)\n", " strata_geoms = []\n", " for j in range(n_strata):\n", " cell = box(x_edges[j], miny, x_edges[j + 1], maxy)\n", " inter = poly.intersection(cell)\n", " if not inter.is_empty and inter.area > 0:\n", " strata_geoms.append(inter)\n", "\n", " strata_gdf = gpd.GeoDataFrame(\n", " {\"stratum_id\": range(len(strata_geoms))},\n", " geometry=strata_geoms,\n", " crs=\"EPSG:5179\"\n", " )\n", " return strata_gdf\n", "\n", "\n", "\n", "\n", "# ===========================\n", "# 4-1. 층화 샘플링 (면적 비례 할당, 기존)\n", "# ===========================\n", "def stratified_sampling(poly_strata_gdf, n_total):\n", " areas = poly_strata_gdf.geometry.area.values\n", " props = areas / areas.sum()\n", "\n", " # 면적 비례 초기 할당\n", " n_alloc = np.round(props * n_total).astype(int)\n", "\n", " # 총합 보정 (sum == n_total 되도록)\n", " diff = n_total - n_alloc.sum()\n", " if diff > 0:\n", " # 더 필요한 만큼, 면적이 큰 층부터 1씩 추가\n", " idx_order = np.argsort(-areas)\n", " for k in range(diff):\n", " n_alloc[idx_order[k % len(n_alloc)]] += 1\n", " elif diff < 0:\n", " # 너무 많이 할당된 경우, 면적이 작은 층부터 1씩 감소 (0 밑으로는 안 가게)\n", " idx_order = np.argsort(areas)\n", " k = 0\n", " while diff < 0:\n", " i = idx_order[k % len(n_alloc)]\n", " if n_alloc[i] > 1: # 최소 1개는 유지\n", " n_alloc[i] -= 1\n", " diff += 1\n", " k += 1\n", "\n", " # 각 층에서 샘플링\n", " all_points = []\n", " stratum_ids = []\n", " for geom, sid, n_i in zip(poly_strata_gdf.geometry,\n", " poly_strata_gdf[\"stratum_id\"],\n", " n_alloc):\n", " pts = random_points_in_polygon(geom, n_i)\n", " all_points.extend(pts)\n", " stratum_ids.extend([sid] * len(pts))\n", "\n", " sample_gdf = gpd.GeoDataFrame(\n", " {\"stratum_id\": stratum_ids},\n", " geometry=all_points,\n", " crs=poly_strata_gdf.crs,\n", " )\n", " return sample_gdf, n_alloc\n", "\n", "\n", "# ===========================\n", "# 4-2. 층화 샘플링 (동일 표본수 할당: n_total을 층 개수로 균등 분배)\n", "# ===========================\n", "def stratified_sampling_equal_n(poly_strata_gdf, n_total):\n", " n_strata = len(poly_strata_gdf)\n", " base = n_total // n_strata\n", " rem = n_total % n_strata\n", "\n", " n_alloc = np.array([base] * n_strata)\n", " for i in range(rem):\n", " n_alloc[i] += 1\n", "\n", " all_points = []\n", " stratum_ids = []\n", " for geom, sid, n_i in zip(poly_strata_gdf.geometry,\n", " poly_strata_gdf[\"stratum_id\"],\n", " n_alloc):\n", " pts = random_points_in_polygon(geom, n_i)\n", " all_points.extend(pts)\n", " stratum_ids.extend([sid] * len(pts))\n", "\n", " sample_gdf = gpd.GeoDataFrame(\n", " {\"stratum_id\": stratum_ids},\n", " geometry=all_points,\n", " crs=poly_strata_gdf.crs,\n", " )\n", " return sample_gdf, n_alloc\n", "\n", "\n", "# ===========================\n", "# 5. 클리핑용 패치 생성\n", "# ===========================\n", "def make_clip_patch_and_poly_web(boundary_gdf_5179):\n", " \"\"\"\n", " 경계 gdf(EPGS:5179)에서 Web Mercator(3857) 기준 폴리곤 + clip용 PathPatch 반환\n", " \"\"\"\n", " gdf_web = boundary_gdf_5179.to_crs(epsg=3857)\n", " poly_web = gdf_web.geometry.unary_union\n", "\n", " # MultiPolygon이면 가장 큰 것 사용\n", " if poly_web.geom_type == \"MultiPolygon\":\n", " poly_web = max(poly_web.geoms, key=lambda g: g.area)\n", "\n", " x, y = poly_web.exterior.coords.xy\n", " coords = np.vstack([x, y]).T\n", "\n", " Path = mpath.Path\n", " path_data = [(Path.MOVETO, coords[0])]\n", " for c in coords[1:]:\n", " path_data.append((Path.LINETO, c))\n", " path_data.append((Path.CLOSEPOLY, coords[0]))\n", "\n", " codes, verts = zip(*path_data)\n", " clip_path = mpath.Path(verts, codes)\n", " clip_patch = mpatches.PathPatch(\n", " clip_path,\n", " facecolor=\"none\",\n", " edgecolor=\"none\"\n", " )\n", "\n", " return clip_patch, poly_web, gdf_web\n", "\n", "\n", "# ===========================\n", "# 6. 시각화 1: 원본 경계만 (샘플 없음)\n", "# ===========================\n", "def plot_base(gdf_boundary_5179,\n", " out_png=\"inha_boundary_base.png\"):\n", "\n", " clip_patch, poly_web, boundary_web = make_clip_patch_and_poly_web(gdf_boundary_5179)\n", "\n", " fig, ax = plt.subplots(figsize=(12, 10))\n", "\n", " # clip 패치\n", " ax.add_patch(clip_patch)\n", "\n", " # 경계(선)\n", " boundary_web.boundary.plot(ax=ax, linewidth=2, edgecolor=\"black\")\n", "\n", " # basemap\n", " ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, alpha=0.9)\n", "\n", " # basemap 클리핑\n", " if ax.images:\n", " ax.images[-1].set_clip_path(clip_patch)\n", "\n", " # 축 범위\n", " minx, miny, maxx, maxy = poly_web.bounds\n", " ax.set_xlim(minx, maxx)\n", " ax.set_ylim(miny, maxy)\n", "\n", " # 바깥 투명\n", " ax.set_facecolor(\"none\")\n", " fig.patch.set_alpha(0.0)\n", "\n", " ax.axis(\"off\")\n", " plt.tight_layout()\n", " plt.savefig(out_png, dpi=300, bbox_inches=\"tight\", transparent=True)\n", " plt.close(fig)\n", " print(f\"✓ 원본 경계 지도 저장: {out_png}\")\n", "\n", "\n", "# ===========================\n", "# 7. 시각화 2: SRS 샘플\n", "# ===========================\n", "def plot_srs(gdf_boundary_5179, srs_points_5179,\n", " out_png=\"inha_sampling_srs.png\"):\n", "\n", " clip_patch, poly_web, boundary_web = make_clip_patch_and_poly_web(gdf_boundary_5179)\n", " pts_web = srs_points_5179.to_crs(epsg=3857)\n", "\n", " fig, ax = plt.subplots(figsize=(12, 10))\n", "\n", " ax.add_patch(clip_patch)\n", "\n", " # 경계(선)\n", " boundary_web.boundary.plot(ax=ax, linewidth=2, edgecolor=\"black\")\n", "\n", " # SRS 점\n", " pts_web.plot(ax=ax, markersize=30, color=\"red\", alpha=0.9)\n", "\n", " # basemap\n", " ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, alpha=0.9)\n", "\n", " # basemap 클리핑\n", " if ax.images:\n", " ax.images[-1].set_clip_path(clip_patch)\n", "\n", " # 축 범위\n", " minx, miny, maxx, maxy = poly_web.bounds\n", " ax.set_xlim(minx, maxx)\n", " ax.set_ylim(miny, maxy)\n", "\n", " # 바깥 투명\n", " ax.set_facecolor(\"none\")\n", " fig.patch.set_alpha(0.0)\n", "\n", " ax.axis(\"off\")\n", " plt.tight_layout()\n", " plt.savefig(out_png, dpi=300, bbox_inches=\"tight\", transparent=True)\n", " plt.close(fig)\n", " print(f\"✓ SRS 샘플 지도 저장: {out_png}\")\n", "\n", "\n", "# ===========================\n", "# 8-1. 시각화 3: 층화 샘플 (격자 + 면적비례)\n", "# ===========================\n", "def plot_strata(gdf_boundary_5179, strata_gdf_5179, sample_gdf_5179,\n", " out_png=\"inha_sampling_strata.png\"):\n", "\n", " clip_patch, poly_web, boundary_web = make_clip_patch_and_poly_web(gdf_boundary_5179)\n", " strata_web = strata_gdf_5179.to_crs(epsg=3857)\n", " sample_web = sample_gdf_5179.to_crs(epsg=3857)\n", "\n", " fig, ax = plt.subplots(figsize=(12, 10))\n", "\n", " ax.add_patch(clip_patch)\n", "\n", " # 층(구역) 영역 색칠\n", " strata_web.plot(\n", " ax=ax,\n", " column=\"stratum_id\",\n", " alpha=0.35,\n", " edgecolor=\"black\",\n", " linewidth=1,\n", " legend=False\n", " )\n", "\n", " # 캠퍼스 전체 경계\n", " boundary_web.boundary.plot(ax=ax, linewidth=2, edgecolor=\"black\")\n", "\n", " # 층화 샘플 점 (층별 색 구분)\n", " sample_web.plot(\n", " ax=ax,\n", " column=\"stratum_id\",\n", " markersize=30,\n", " cmap=\"tab10\",\n", " legend=False\n", " )\n", "\n", " # basemap\n", " ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, alpha=0.9)\n", "\n", " # basemap 클리핑\n", " if ax.images:\n", " ax.images[-1].set_clip_path(clip_patch)\n", "\n", " # 축 범위\n", " minx, miny, maxx, maxy = poly_web.bounds\n", " ax.set_xlim(minx, maxx)\n", " ax.set_ylim(miny, maxy)\n", "\n", " # 바깥 투명\n", " ax.set_facecolor(\"none\")\n", " fig.patch.set_alpha(0.0)\n", "\n", " ax.axis(\"off\")\n", " plt.tight_layout()\n", " plt.savefig(out_png, dpi=300, bbox_inches=\"tight\", transparent=True)\n", " plt.close(fig)\n", " print(f\"✓ 층화 샘플 지도 저장: {out_png}\")\n", "\n", "\n", "# ===========================\n", "# 8-2. 시각화 4: 층화 샘플 (동일면적 strata + 동일 표본수)\n", "# ===========================\n", "def plot_strata_equal(gdf_boundary_5179, strata_gdf_5179, sample_gdf_5179,\n", " out_png=\"inha_sampling_strata_equalarea.png\"):\n", "\n", " clip_patch, poly_web, boundary_web = make_clip_patch_and_poly_web(gdf_boundary_5179)\n", " strata_web = strata_gdf_5179.to_crs(epsg=3857)\n", " sample_web = sample_gdf_5179.to_crs(epsg=3857)\n", "\n", " fig, ax = plt.subplots(figsize=(12, 10))\n", "\n", " ax.add_patch(clip_patch)\n", "\n", " # 층(구역) 영역 색칠\n", " strata_web.plot(\n", " ax=ax,\n", " column=\"stratum_id\",\n", " alpha=0.35,\n", " edgecolor=\"black\",\n", " linewidth=1,\n", " legend=False\n", " )\n", "\n", " # 캠퍼스 전체 경계\n", " boundary_web.boundary.plot(ax=ax, linewidth=2, edgecolor=\"black\")\n", "\n", " # 층화 샘플 점 (층별 색 구분)\n", " sample_web.plot(\n", " ax=ax,\n", " column=\"stratum_id\",\n", " markersize=30,\n", " cmap=\"tab10\",\n", " legend=False\n", " )\n", "\n", " # basemap\n", " ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, alpha=0.9)\n", "\n", " # basemap 클리핑\n", " if ax.images:\n", " ax.images[-1].set_clip_path(clip_patch)\n", "\n", " # 축 범위\n", " minx, miny, maxx, maxy = poly_web.bounds\n", " ax.set_xlim(minx, maxx)\n", " ax.set_ylim(miny, maxy)\n", "\n", " # 바깥 투명\n", " ax.set_facecolor(\"none\")\n", " fig.patch.set_alpha(0.0)\n", "\n", "\n", " ax.axis(\"off\")\n", " plt.tight_layout()\n", " plt.savefig(out_png, dpi=300, bbox_inches=\"tight\", transparent=True)\n", " plt.close(fig)\n", " print(f\"✓ 동일면적 층화 샘플 지도 저장: {out_png}\")\n", "\n", "\n", "# ===========================\n", "# 9. 메인 실행\n", "# ===========================\n", "if __name__ == \"__main__\":\n", " # 1) 경계 불러오기\n", " boundary_gdf, boundary_poly = load_inha_boundary(\"inha_boundary.shp\")\n", "\n", " # 2) SRS 샘플링 (예: 40개)\n", " n_srs = 60\n", " srs_points = random_points_in_polygon(boundary_poly, n_srs)\n", " srs_gdf = gpd.GeoDataFrame(\n", " {\"design\": [\"SRS\"] * len(srs_points)},\n", " geometry=srs_points,\n", " crs=\"EPSG:5179\"\n", " )\n", "\n", " # 3-1) 격자 층(구역) 만들고 층화 샘플링 (면적비례)\n", " strata_gdf = make_strata(boundary_poly, n_rows=2, n_cols=3)\n", " strat_sample_gdf, n_alloc = stratified_sampling(strata_gdf, n_total=60)\n", " print(\"격자 층별 표본 수 할당 (면적비례):\", n_alloc)\n", "\n", " # 3-2) 동일면적 strata 만들고 층화 샘플링 (동일 표본수)\n", " equal_strata_gdf = make_equal_area_strata(boundary_poly, n_strata=6)\n", " equal_strat_sample_gdf, n_alloc_equal = stratified_sampling_equal_n(equal_strata_gdf, n_total=60)\n", " print(\"동일면적 층별 표본 수 할당 (동일 표본수):\", n_alloc_equal)\n", "\n", " # 4) 원본 경계 지도 (바깥 투명)\n", " plot_base(boundary_gdf, out_png=\"inha_boundary_base.png\")\n", "\n", " # 5) SRS 지도\n", " plot_srs(boundary_gdf, srs_gdf,\n", " out_png=\"inha_sampling_srs.png\")\n", "\n", " # 6) 층화 지도 (격자 + 면적비례)\n", " plot_strata(boundary_gdf, strata_gdf, strat_sample_gdf,\n", " out_png=\"inha_sampling_strata.png\")\n", "\n", " # 7) 층화 지도 (동일면적 strata + 동일 표본수)\n", " plot_strata_equal(boundary_gdf, equal_strata_gdf, equal_strat_sample_gdf,\n", " out_png=\"inha_sampling_strata_equalarea.png\")\n", "\n", " # 8) 좌표 확인용 출력\n", " print(\"\\n[SRS 샘플 좌표 (EPSG:5179)]\")\n", " print(srs_gdf.head())\n", "\n", " print(\"\\n[격자 층화 샘플 좌표 (EPSG:5179)]\")\n", " print(strat_sample_gdf.head())\n", "\n", " print(\"\\n[동일면적 층화 샘플 좌표 (EPSG:5179)]\")\n", " print(equal_strat_sample_gdf.head())\n" ] } ], "metadata": { "kernelspec": { "display_name": "torch_ev", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 5 }