Earth Engine & Geemap: 토지피복 분류(supervised classification)

안녕하세요? 이번 글은 Google Earth Engine과 geemap을 활용한 토지피복 분류를 실습해 보도록 하겠습니다. 태안해안국립공원 지역의 Landsat 8호 위성영상을 대상으로 감독 분류(supervised classification)을 다뤄보는 내용입니다. 해당 내용은 geemap 공식 튜토리얼 내용을 일부 수정, 보완한 것입니다.

32 supervised classification - geemap

일단, ee와 geemap 모듈을 가져오고, Earth Engine을 초기화하는 작업을 수행합니다.

import ee
import geemap

# Earth Engine 초기화
ee.Initialize()

geemap 라이브러리에서 Map 클래스를 사용하여 빈 지도를 변수 Map에 할당합니다.

Map = geemap.Map()
Map

WDPA: World Database on Protected Areas (polygons)에서 태안해안국립공원 경계를 선택고 중앙점을 계산합니다. 중앙점의 위경도 좌표는 Landsat 8호 위성영상을 검색하는데 사용하고자 합니다.

WDPA: World Database on Protected Areas (polygons)  |  Earth Engine Data Catalog  |  Google for Developers

# 태안해안국립공원 영역을 추출하는 함수 정의
def get_boundary():
    roi = ee.FeatureCollection("WCMC/WDPA/current/polygons") \
        .filter(ee.Filter.eq("NAME", "Taeanhaean")) \
        .geometry()
    return roi

# 태안해안국립공원 경계 추출
region = get_boundary()

# 중앙점 계산
centroid = region.centroid()

# 중앙점의 위도와 경도 값 얻기
lon_lat = centroid.getInfo()['coordinates']
longitude = lon_lat[0]
latitude = lon_lat[1]

# 위도와 경도 값을 소수점 넷째 자리까지 출력
print("중앙점의 위도: {:.4f}".format(latitude))
print("중앙점의 경도: {:.4f}".format(longitude))

중앙점의 위도: 36.6125
중앙점의 경도: 126.2256

USGS Landsat 8 Collection 2 Tier 1 TOA Reflectance에서 Landsat 8호 위성영상을 검색합니다. 태안해안국립공원 지역 주변의 이미지만 검색하고 2022년 1월 1일부터 2022년 12월 31일까지 이미지 중 선택합니다. 또한 구름 피복 기준으로 영상을 정렬하여 구름 덮인 정도가 가장 적은 이미지로부터 1부터 7까지의 밴드만 선택합니다.

USGS Landsat 8 Collection 2 Tier 1 TOA Reflectance  |  Earth Engine Data Catalog  |  Google for Developers

point = ee.Geometry.Point([126.2256, 36.6125])

image = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_TOA')
    .filterBounds(point)
    .filterDate('2022-01-01', '2022-12-31')
    .sort('CLOUD_COVER')
    .first()
    .select('B[1-7]')
)

vis_params = {'min': 0, 'max': 0.4, 'bands': ['B5', 'B4', 'B3']}

# Landsat-8 데이터 추가
Map.centerObject(point, 9)
Map.addLayer(image, vis_params, "Landsat-8")
Map

선택한 이미지의 해당 날짜와 구름 피복 정보를 확인해 봅니다.

ee.Date(image.get('system:time_start')).format('YYYY-MM-dd').getInfo()

'2022-10-27'

image.get('CLOUD_COVER').getInfo()

0.02

이번에는 ESA WorldCover 10m v200 이미지 컬렉션을 추가하고, Landsat 8호 위성영상과 동일한 크기로 자르겠습니다.

ESA WorldCover 10m v200  |  Earth Engine Data Catalog  |  Google for Developers

# ESA WorldCover 10m v200
esa_worldcover = ee.ImageCollection('ESA/WorldCover/v200').first()
esa_worldcover_kr = esa_worldcover.clip(image.geometry())

Map.addLayer(esa_worldcover_kr, {}, 'ESA WorldCover')
Map

ESA WorldCover 10m v200를 참조하여 학습 데이터를 생성해 보겠습니다. 이미지 영역에서 30m급 해상도를 기준으로 5,000개 픽셀을 샘플로 추출합니다. 0난수 생성 시드(seed) 파라미터는 동일한 시드를 사용하면 재현 가능한 결과를 얻을 수 있습니다. 지오메트리 정보(geometries)를 가져올지 여부를 True로 설정했으므로, 샘플의 각 지오메트리 정보도 포함될 것입니다.

# 학습 데이터 생성
points = esa_worldcover_kr.sample(
    **{
        'region': image.geometry(),
        'scale': 30,
        'numPixels': 5000,
        'seed': 0,
        'geometries': True,
    }
)

Map.addLayer(points, {}, 'training', False)

print(points.size().getInfo())

4103

print(points.first().getInfo())

{'type': 'Feature', 'geometry': {'type': 'Point', 'coordinates': [125.95303141182198, 36.149676107086066]}, 'id': '0', 'properties': {'Map': 80}}

토지피복을 예측하는 분류기를 만드는 과정입니다. 예측에 사용할 특성(feature) 밴드들을 지정합니다. 각 밴드는 특정 스펙트럼 대역을 나타냅니다.  학습 데이터의 각 지점은 토지피복 클래스 레이블을 'Map' 속성에 저장하고 있을 것입니다. 각 지점을 이미지 위에 오버레이하여 훈련 데이터를 샘플링합니다. 이때 특성 값과 토지피복 값이 모두 포함되도록 설정합니다. 여기서는 CART(Classification and Regression Trees) 분류기를 훈련했습니다. 이렇게 학습된 분류기 trained를 사용하여 새로운 데이터의 토지피복을 예측할 수 있습니다. 

# 예측에 사용하는 밴드
bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7']


# 토지피복 레이블 저장 테이블의 속성
label = 'Map'

# 이미지 위에 점들을 오버레이하여 훈련 데이터 얻기
training = image.select(bands).sampleRegions(
    **{'collection': points, 'properties': [label], 'scale': 30}
)

# 기본 매개변수를 사용하여 CART 분류기 훈련
trained = ee.Classifier.smileCart().train(training, label, bands)

print(training.first().getInfo())

{'type': 'Feature', 'geometry': None, 'id': '0_0', 'properties': {'B1': 0.1613909751176834, 'B2': 0.13543030619621277, 'B3': 0.08888012915849686, 'B4': 0.05841150879859924, 'B5': 0.03379363194108009, 'B6': 0.011541628278791904, 'B7': 0.006362282671034336, 'Map': 80}}

학습에 사용된 동일한 밴드들을 선택하여 이미지를 분류합니다. 이전에 학습된 분류기 trained를 사용하여 이미지를 분류하고, 분류 결과를 result 변수에 할당합니다.

# 훈련에 사용된 동일한 밴드로 이미지 분류
result = image.select(bands).classify(trained)

# 랜덤한 색상으로 클러스터 표시
Map.addLayer(result.randomVisualizer(), {}, 'classified')
Map

현재 추출 결과는 각 분류 클래스에 랜덤한 색상을 적용하여 시각화한 것입니다. 이번에는 토지피복 클래스와 색상을 지정해 보겠습니다.

class_values = esa_worldcover_kr.get('Map_class_values').getInfo()
class_values

[10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100]

class_palette = esa_worldcover_kr.get('Map_class_palette').getInfo()
class_palette

['006400',
 'ffbb22',
 'ffff4c',
 'f096ff',
 'fa0000',
 'b4b4b4',
 'f0f0f0',
 '0064c8',
 '0096a0',
 '00cf75',
 'fae6a0']

landcover = result.set('classification_class_values', class_values)
landcover = landcover.set('classification_class_palette', class_palette)

Map.addLayer(landcover, {}, 'Land cover')
Map

분류 결과 레이어의 투명도를 조정하고, 미리 정해진 ESA_WorldCover 범례를 지도에 추가해 보겠습니다.

print('Change layer opacity:')
cluster_layer = Map.layers[-1]
cluster_layer.interact(opacity=(0, 1, 0.1))

Map.add_legend(builtin_legend='ESA_WorldCover')
Map

geemap 라이브러리를 통해 Earth Engine에 생성된 이미지를 로컬 디렉토리에 내보낼 수 있습니다. 주석 처리된 코드 부분은 이미지를 Google Drive에 저장하고 싶은 경우에 사용할 수 있습니다.

import os

out_dir = 'D:/GEODATA'
out_file = os.path.join(out_dir, 'landcover.tif')

geemap.ee_export_image(landcover, filename=out_file, scale=100)

"""
geemap.ee_export_image_to_drive(
    landcover, description='landcover', folder='export', scale=30
)
"""

1912_geemap_Earth Engine을 활용한 토지피복 분류.ipynb

0.01MB