[python GIS] GDAL기반 특정 위치에서 래스터 파일의 값을 추출하기

이번 글은 래스터 파일에서 위치포인트를 갖는 shp파일을 이용해서 래스터 값을 추출하는 코드를 설명한다.

아래 코드는 jupyter notebook에서 실행한다.

 

# 관련 패키지 불러오기

from osgeo import gdal, ogr, osr

 

# 래스터 파일 경로
raster_file_path = "env_data/res_30m_f/seoul_dem_final.tif"

 

# 포인트 데이터셋 (Shapefile) 경로
point_shp_path = "env_data/sinkhole_loc/서울시_싱크홀_내역_e5186.shp"

# 출력 CSV 파일 경로
output_csv_path = " env_data /path_to_output.csv"

 

# 출력폴더 생성
output_floder = "results"
os.makedirs(output_floder, exist_ok=True)

# 래스터 파일 열기
raster_ds = gdal.Open(raster_file_path)

# 포인트 데이터셋 열기
point_ds = ogr.Open(point_shp_path)
point_layer = point_ds.GetLayer()

# 래스터 정보 가져오기
geotransform = raster_ds.GetGeoTransform()
x_origin = geotransform[0]
y_origin = geotransform[3]
pixel_width = geotransform[1]
pixel_height = geotransform[5]

# CSV 파일 열기 (값 저장용)
csv_file = open(f'{output_floder}/{output_csv_path}', "w")
csv_file.write("X,Y,Value\n")

# 포인트 데이터셋의 각 포인트에서 래스터 값 추출
for point_feature in point_layer:
    point_geometry = point_feature.GetGeometryRef()
    x = point_geometry.GetX()
    y = point_geometry.GetY()

    # 좌표를 픽셀 좌표로 변환
    x_pixel = int((x - x_origin) / pixel_width)
    y_pixel = int((y - y_origin) / pixel_height)

    # 래스터 값 읽기
    raster_band = raster_ds.GetRasterBand(1)  # 여기서는 첫 번째 밴드를 사용
    value = raster_band.ReadAsArray(x_pixel, y_pixel, 1, 1)[0][0]

    # CSV 파일에 기록
    csv_file.write(f"{x},{y},{value}\n")

# 파일 닫기
csv_file.close()

print(f"값이 CSV 파일에 저장되었습니다: {output_csv_path}")