Google Earth Engine: NASADEM 이미지 및 지형정보 생성방법 소개

안녕하세요? 이번 글은 Google Earth Engine에서 NASADEM 이미지를 사용하는 방법을 소개해 보겠습니다.

NASADEM: NASA NASADEM Digital Elevation 30m  |  Earth Engine Data Catalog  |  Google for Developers

NASADEM: NASA NASADEM Digital Elevation 30m

NASADEM은 SRTM 데이터를 재처리한 것으로, ASTER GDEM, ICESat GLAS, PRISM 데이터셋과 같은 보조 데이터를 통합하여 정확도를 향상시킨 데이터입니다. 가장 중요한 처리 개선사항은 향상된 위상 해제(phase unwrapping)를 통한 결측값 감소와 ICESat GLAS 데이터를 제어 지점으로 사용한 것입니다.

 

NASADEM은 총3개 레이어로 구성되어 있습니다: elevation은 고도를, num은 데이터 소스 및 소스 촬영 자료의 수를 나타내는 인덱스를, swb는 업데이트된 SRTM 수체 데이터를 나타냅니다.

# NASADEM
dataset = ee.Image('NASA/NASADEM_HGT/001')

# NASADEM 레이어 목록 출력
layers = dataset.bandNames().getInfo()
for layer in layers:
    print(layer)

elevation
num
swb

대한민국 지역을 중심으로 고도 레이어를 시각화해 보겠습니다.

# 고도 밴드 선택
elevation = dataset.select('elevation')

# 고도가 0 이하인 부분은 투명하게 설정
elevation_masked = elevation.updateMask(elevation.gt(0))

# 고도 시각화 속성
elevation_vis = {
    'min': 0,
    'max': 2000,
}

# 지도 생성
m = geemap.Map(layout={'height':'400px', 'width':'800px'})

# 흰색 배경 이미지
background = ee.Image(1)
m.addLayer(background, {'min': 0, 'max': 1}, 'Background')

# 고도
m.addLayer(elevation_masked, elevation_vis, 'Elevation')

# 대한민국 지역 중심
countries = ee.FeatureCollection("WM/geoLab/geoBoundaries/600/ADM0")
region = countries.filter(ee.Filter.eq('shapeName', 'Korea, South'))
m.centerObject(region, 6)
m

ee.Algorithms.Terrain

ee.Algorithms.Terrain은 지형 DEM에서 경사(slope), 방위(aspect), 그리고 간단한 음영기복(hillshade)을 계산하는 알고리즘입니다. 아래 코드로 지형정보 레이어 목록을 확인하실 수 있습니다.

ee.Algorithms.Terrain  |  Google Earth Engine  |  Google for Developers

# NASADEM 데이터를 불러와서 지형정보 생성
terrain = ee.Algorithms.Terrain(ee.Image('NASA/NASADEM_HGT/001'))

# 지형정보 레이어 목록 출력
layers = terrain.bandNames().getInfo()
for layer in layers:
    print(layer)

elevation
num
swb
slope
aspect
hillshade

ee.Algorithms.Terrain에서 고도(Elevation) 레이어를 시각화해 보겠습니다.

# 수역 마스크 생성
watermask = terrain.select('elevation').gt(0)

# 수역 픽셀을 마스킹하고 관심지역으로 클리핑
predictors = terrain.updateMask(watermask).clip(region)
predictors = predictors.select(['elevation', 'slope', 'aspect'])

# 고도
m = geemap.Map(layout={'height':'400px', 'width':'800px'})
vis_params = {'bands':['elevation'], 'min': 0, 'max': 1800, 'palette': cm.palettes.terrain}
m.addLayer(predictors, vis_params, 'Elevation')
m.add_colorbar(vis_params, label="Elevation (m)", orientation="vertical", layer_name="elevation")
m.centerObject(region, 6)
m

경사 레이어 시각화는 다음과 같습니다.

# 경사
m = geemap.Map(layout={'height':'400px', 'width':'800px'})
vis_params = {'bands':['slope'], 'min': 0, 'max': 25, 'palette': cm.palettes.RdYlGn_r}
m.addLayer(predictors, vis_params, 'Slope')
m.add_colorbar(vis_params, label="Slope", orientation="vertical", layer_name="slope")
m.centerObject(region, 6)
m

경사 방향 시각화는 다음과 같습니다.

# 경사 방향
m = geemap.Map(layout={'height':'400px', 'width':'800px'})
vis_params = {'bands':['aspect'], 'min': 0, 'max': 360, 'palette': cm.palettes.rainbow}
m.addLayer(predictors, vis_params, 'Aspect')
m.add_colorbar(vis_params, label="Aspect", orientation="vertical", layer_name="aspect")
m.centerObject(region, 6)
m

현재 데이터셋의 EPSG 코드는 EPSG:4326, 공간 해상도는 30미터입니다.

# 투영 정보 가져오기
projection = predictors.projection()

# EPSG 코드 및 해상도 가져오기
epsg_code = projection.crs().getInfo()        # 투영의 EPSG 코드
scale = projection.nominalScale().getInfo()   # 해상도(미터 단위)

# 결과 출력
print('Predictors의 EPSG 코드:', epsg_code)
print('Predictors의 공간 해상도:', scale, '미터')

Predictors의 EPSG 코드: EPSG:4326
Predictors의 공간 해상도: 30.922080775909325 미터

이어서 고도, 경사, 경사 방향 레이어에 대한 피어슨 상관계수를 계산해 보겠습니다. 피어슨 상관분석은 요인들 간의 잠재적 다중공선성을 확인하는 데 활용되며, 상관계수 |r|이 0.8 이상인 요인들은 다중공선성 문제를 줄이기 위해 일반적으로 제외됩니다.

# 30m 스케일로 픽셀 5000개 샘플링
sample = predictors.sample(region=region, scale=30, numPixels=5000).getInfo()
data = pd.DataFrame([feat["properties"] for feat in sample["features"]])

# 피어슨 상관계수 계산
corr_matrix = data.corr(method="pearson")

# 상삼각 행렬의 값만 추출하여 중복을 제거한 상관계수 출력
def get_upper_triangle(corr_matrix):
    upper_triangle = corr_matrix.where(
        np.triu(np.ones(corr_matrix.shape), k=1).astype(bool)
    )
    return upper_triangle

upper_triangle_corr = get_upper_triangle(corr_matrix)

# 상삼각 행렬을 1차원으로 변환하고 높은 순으로 정렬
sorted_corr = (
    upper_triangle_corr.unstack().dropna().sort_values(ascending=False)
)

print("피어슨 상관계수 (높은 순):")
print(sorted_corr)

피어슨 상관계수 (높은 순):
slope      elevation    0.555488
           aspect       0.079534
elevation  aspect       0.057541
dtype: float64

피어슨 상관계수를 시각화한 결과는 다음과 같습니다.

# 히트맵 생성
plt.figure(figsize=(5, 4))
plt.title("Pearson Correlation Heatmap")
heatmap = plt.imshow(corr_matrix, cmap='coolwarm', interpolation='none', aspect='auto')
plt.colorbar(label='Correlation Coefficient')
plt.xticks(np.arange(len(corr_matrix.columns)), corr_matrix.columns, rotation=45)
plt.yticks(np.arange(len(corr_matrix.columns)), corr_matrix.columns)

# 상관계수 값 추가
for i in range(len(corr_matrix.columns)):
    for j in range(len(corr_matrix.columns)):
        plt.text(j, i, f'{corr_matrix.iloc[i, j]:.2f}',
                 ha='center', va='center', color='black')

plt.tight_layout()
plt.show()