GIS

Python: 폴리곤 내에서 포인트 무작위 표집

유병혁 2022. 10. 24. 17:21

이 글은 국립공원공단에서 <AI 기반 야생동물 예측 서비스> 프로젝트를 진행한 2022년 데이터 청년 캠퍼스(빅리더)팀의 소스코드를 재구성했습니다.

안녕하세요? 이번 글은 Python을 통해 폴리곤 내에서 원하는 개수만큼 포인트를 무작위 표집(random sampling)하는 과정을 정리해 보겠습니다. 다양한 방법이 존재할 텐데요, 이 글에서 소개하는 방법은 GIS StackExchange에서 논의된 두가지 방법을 다룹니다.

 

Randomly sample from geopandas DataFrame in Python

I am reading a shapefile as geopandas DataFrame and them using pandas subset method to select a region. geodata = gpd.read_file(bayshp) geodata.dtypes geodata.head(10) OBJECTID FIPSSTCO ...

gis.stackexchange.com

먼저, 필요한 모듈을 호출해 봅니다.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
import shapely.geometry

폴리곤으로 사용할 우리나라 행정구역 데이터를 열어보겠습니다.

gdf = gpd.GeoDataFrame.from_file("DATA/ADM_KOR.gpkg")
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
gdf.plot(ax=ax, color='lightgray', edgecolor='white');

행정구역은 다수의 객체로 구성되는데요, 전체 도형을 합집합(union)으로 표현하는 폴리곤을 생성해 보겠습니다.

# 합집합(union) 도형
polygon = gdf['geometry'].unary_union
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
gpd.GeoSeries(polygon).plot(ax=ax, color='lightgray', edgecolor='white');

폴리곤의 영역은 다음과 같이 정의할 수 있습니다. 이것은 행정구역을 아우르는 경계상자(envolope) 영역입니다.

min_x, min_y, max_x, max_y = polygon.bounds # 폴리곤 영역

자, 그럼 포인트 무작위 표집을 만들어 보겠습니다. 여기서 size는 포인트 개수, overestimate는 실제보다 포인트 개수를 과대 산정하기 위한 파라미터, ratio는 폴리곤과 경계상자 면적의 비율을 나타냅니다. 만약 폴리곤 면적과 경계상자 면적이 일치한다면 ratio는 1이 되겠죠?!

size = 100 # 포인트 개수
overestimate = 2 # 과대 산정
ratio = polygon.area / polygon.envelope.area # 면적 비율 = 폴리곤 면적 / 경계상자 면적

다음과 같이 size, overestimate, ratio가 반영된 무작위 표집을 진행합니다.

samples = np.random.uniform((min_x, min_y), (max_x, max_y), (int(size / ratio * overestimate), 2))
multipoint = shapely.geometry.MultiPoint(samples)

결과는 다음과 같습니다.

base = gpd.GeoSeries(polygon).plot(color='lightgray', edgecolor='white')
gpd.GeoSeries(multipoint).plot(ax=base, marker='o', color='red', markersize=5);

폴리곤과 중첩되는 포인트만을 선택해 줍니다.

multipoint = multipoint.intersection(polygon)
base = gpd.GeoSeries(polygon).plot(color='lightgray', edgecolor='white')
gpd.GeoSeries(multipoint).plot(ax=base, marker='o', color='red', markersize=5);

포인트 개수는 100개로 정의되어 있습니다. 즉, 전체 개수 중 100개를 임의 선택해서 추출해 줍니다.

samples = np.array(multipoint)
df = pd.DataFrame(samples[np.random.choice(len(samples), size)])
df.columns = ['lon', 'lat']
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lon, df.lat))
gdf_pt.head()

물론 반복문을 통해 포인트 생성시마다 폴리곤에 포함되는지 여부를 체크하는 아래와 같은 함수를 정의할 수도 있습니다.

from shapely.geometry import Point

# 단일 폴리곤에서 무작위 표본 추출
def random_points_in_polygon(number, polygon):
    points = []
    min_x, min_y, max_x, max_y = polygon.bounds
    i= 0
    while i < number:
        point = Point(np.random.uniform(min_x, max_x), np.random.uniform(min_y, max_y))
        if polygon.contains(point):
            points.append(point)
            i += 1
    return points

문제는 속도인데요, 반복문을 쓰는 만큼 100개 포인트를 추출하는데 더딘 처리속도를 체감할 수 있습니다.

polygon = gdf['geometry'].unary_union
samples = random_points_in_polygon(size, polygon)
multipoint = shapely.geometry.MultiPoint(samples)
base = gpd.GeoSeries(polygon).plot(color='lightgray', edgecolor='white')
gpd.GeoSeries(multipoint).plot(ax=base, marker='o', color='red', markersize=5);

앞서 소개한 과정을 하나의 함수로 정의합니다. 여기까지 소스코드는 GIS Stack Exchange에 소개되어 있습니다. 아래 소스코드에서 저는 while 문을 수정하고자 하는데요, overestimate가 포인트 개수를 조정하는데 그럼에도 불구하고 목표하는 개수에 도달하지 못하는 경우 다시 반복문으로 돌아가기 때문입니다.

# 지오시리즈에서 무작위 표본 추출
def sample_geoseries(geoseries, size, overestimate=2):
    polygon = geoseries.unary_union
    min_x, min_y, max_x, max_y = polygon.bounds
    ratio = polygon.area / polygon.envelope.area
    samples = np.random.uniform((min_x, min_y), (max_x, max_y), (int(size / ratio * overestimate), 2))
    multipoint = shapely.geometry.MultiPoint(samples)
    multipoint = multipoint.intersection(polygon)
    samples = np.array(multipoint)
    while samples.shape[0] < size:
        samples = np.concatenate([samples, random_points_in_polygon(size, polygon)])
    return samples[np.random.choice(len(samples), size)]

 

points = sample_geoseries(gdf['geometry'], 100)
df = pd.DataFrame(points)
df.columns = ['lon', 'lat']
gdf_pt = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lon, df.lat))
gdf_pt.head()

아래와 같이 함수를 재정리했고 결과는 다음과 같습니다.

def random_points_in_gdf(gdf, size, overestimate=2):
    polygon = gdf['geometry'].unary_union
    min_x, min_y, max_x, max_y = polygon.bounds
    ratio = polygon.area / polygon.envelope.area
    samples = np.random.uniform((min_x, min_y), (max_x, max_y), (int(size / ratio * overestimate), 2))
    multipoint = shapely.geometry.MultiPoint(samples)
    multipoint = multipoint.intersection(polygon)
    samples = np.array(multipoint)
    points = samples[np.random.choice(len(samples), size)]
    df = pd.DataFrame(points, columns=['lon', 'lat'])
    return gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lon, df.lat))
base = gpd.GeoSeries(polygon).plot(color='lightgray', edgecolor='white')
gdf_pt = random_points_in_gdf(gdf, 100)
gdf_pt.plot(ax=base, marker='o', color='red', markersize=5);

해당 소스코드는 아래 링크에서 확인하실 수 있습니다.

 

GitHub - osgeokr/sdm-tools

Contribute to osgeokr/sdm-tools development by creating an account on GitHub.

github.com