Rasterio와 Geopandas를 이용한 Point Sampling

QGIS에서 Point Sampling 플러그인을 이용하면 포인트 위치에 해당하는 Raster의 값을 구할 수 있습니다. 예를 들어 DEM 이라면 표고값을 구할 수 있겠죠.

이 부분에 대한 내용은 정보린 님의 블로그에서 자세하게 잘 나와 있습니다.

Point Sampling Tool

그런데, 구글링하다보니 유병혁님의 블로그에서 GeoPandas와 Rasterio를 이용해서 Point Sampling 기능을 구현한 내용이 있어서 따라 만들어봤습니다.

https://foss4g.tistory.com/1752

간만에 Raster 내용 올리면서, Point 데이터 위치에 해당하는 DEM의 표고값을 Point에 받아오는 내용을 공유해 봅니다.

먼저 Rasterio를 이용해서 DEM 데이터를 준비합니다.

import rasterio

# Rasterio로 래스터 데이터 읽기
dataset = rasterio.open(r'data/raster/dem30.tif')

다음 Geopandas를 이용해서 Point 데이터를 준비합니다.

import geopandas as gpd 

# geopandas로 벡터 포인트 읽기 
gdf = gpd.read_file('data/shp/festival_pt5174.shp', encoding='cp949')

두 데이터는 동일한 평면직각좌표계인 상태입니다.

다음은 포인트의 좌표에 대응되는 Raster의 값을 받는 사용자 함수입니다.

# 포인트 샘플링을 수행하는 함수
def sample_elevation(dataset, coords):
    row, col = zip(*[dataset.index(*coord) for coord in coords])
    values = dataset.read(1)[row, col]
    return values

포인트의 좌표를 리스트 형태로 만들어서 위의 샘플링 함수에 적용합니다.

# 포인트의 좌표를 리스트로 추출합니다.
coord_list = [(x,y) for x,y in zip(gdf['geometry'].x , gdf['geometry'].y)]

# 좌표에 대한 샘플링 실행
elevation_samples = sample_elevation(dataset, coord_list) 

elevation_samples의 값은 아래와 같은 형태입니다.

array([ 22, 16, 102, 40, 41, 69, 6, 24, 30, 18, 76, 19, 27, 9, 19, 20, 61, 6, 14, 45, 39, 43, 101, 75, 42, 42], dtype=uint16)

이 결과값 리스트를 포인트 GeoDataframe의 새로운 항목에 붙여줍니다.

gdf['height'] = elevation_samples

gdf 속성을 보면 height 항목에 DEM 표고값이 들어갔습니다.

포인트 Geodataframe에 부여된 DEM의 표고값

지도로 시각화해보면 맞게 적용된 것을 알 수 있습니다.

포인트의 색상이 포인트 샘플링한 DEM의 표고값

꼭 DEM이 아니더라도 Raster에 적용할 수 있겠습니다.