R을 활용한 원격탐사 데이터 처리 및 시각화

그리고 인터랙티브 시각화

Author

이상일, 전혜민, 김세창, 김우형

Published

July 27, 2025

1 원격탐사 데이터 다루기

1.1 실습 준비

첫째, RStudio에서 실습을 위한 프로젝트(project)를 실행한다.

둘째, 실습을 실행할 R 스크립트 파일을 생성한다.

셋째, 필수적인 패키지를 불러온다. 만일 패키지가 인스톨 되어 있지 않다면, 오른쪽 아래 윈도우의 Packages 탭이나 install.packages() 함수를 사용하여 먼저 인스톨한다.

install.packages(c("tidyverse", "sf", "terra", "RStoolbox", "tidyterra", "RcolorBrewer"))
library(tidyverse) # R로 데이터사이언스를 하기 위한 패키지
library(sf) # 벡터 데이터를 위한 패키지
library(terra) # 래스터 데이터를 위한 패키지
library(RStoolbox) # 원격탐사 분석용 전문 패키지
library(tidyterra) # ggplot2 기반의 래스터 시각화 패키지
library(RColorBrewer) # 브루어 컬러를 사용하기 위한 패키지

1.2 Landsat 데이터 취득

1.2.1 데이터 내려받기

모든 실습은 Landsat 8/9에서 촬영한 영상을 사용한다. Landsat 영상을 다운받는 절차를 간략히 정리하면 다음과 같다.

  • USGS가 운영하는 EarthExplore 사이트에 접속

  • 로그인(계정이 없다면 우선 계정 만들기)

  • 원하는 도시 검색

  • 해당 도시의 Landsat 8/9 Level-2 데이터 다운로드

해당 웹사이트에서 데이터를 다운로드 받는 방법에 대한 동영상을 프로그램 웹사이트에 올려 두었으니 이를 참고하면 된다. 원하는 도시의 ’영어 이름’을 통해 검색할 수도 있고, ’PATH/ROW’를 통해 검색할 수도 있다.

이번 실습에서는 미리 다운로드된 데이터를 사용한다. 제공된 landsat zip 파일을 프로젝트 폴더에 압축해제한다.

이 중 LC08_L2SP_116034_20250426_20250429_02_T1 폴더에 다운로드된 landsat 데이터가 들어가 있다. 22개 정도의 파일이 들어가 있을텐데, 폴더 이름을 포함하여 개별 파일이 무엇을 의미하는지는 다음의 웹사이트를 참고할 수 있다.

다음의 두 가지가 중요하다.

  • TIF 파일 중 B와 숫자가 결합된 이름으로 끝나는 파일이 밴드에 해당한다. Level-2 데이터는 Level-1 데이터와 달리 1-11의 모든 밴드를 포함하지 않고 1-7과 10만을 가지고 있다.

  • MTL로 끝나는 텍스트 파일에 메타데이터가 포함되어 있다.

1.2.2 데이터 정리하기

밴드(1-7, 10)만 불러들여 하나의 결합(stack) 파일을 만든다. folder_path에는 Landsat 데이터가 포함된 폴더의 경로를 정확히 입력한다. 최종적으로 완성된 결합 파일의 정보(해상도(30m), EPSG, 구성 파일 등)를 확인한다.

결합 파일을 여러 가지의 밴드를 한꺼번에 포함하고 있는 파일로, 추후 분석 및 시각화를 위해서는 이 파일에서 필요한 밴드만을 추출하고 조합할 수 있다.

folder_path <- "landsat/LC08_L2SP_116034_20250426_20250429_02_T1" # 각자 받은 폴더 입력
band_files <- list.files(
    folder_path,
    pattern = "B[0-9]+\\.TIF$", # B[숫자].TIF로 끝나는 파일들
    full.names = TRUE
    ) |> sort()
band_files
[1] "landsat/LC08_L2SP_116034_20250426_20250429_02_T1/LC08_L2SP_116034_20250426_20250429_02_T1_SR_B1.TIF" 
[2] "landsat/LC08_L2SP_116034_20250426_20250429_02_T1/LC08_L2SP_116034_20250426_20250429_02_T1_SR_B2.TIF" 
[3] "landsat/LC08_L2SP_116034_20250426_20250429_02_T1/LC08_L2SP_116034_20250426_20250429_02_T1_SR_B3.TIF" 
[4] "landsat/LC08_L2SP_116034_20250426_20250429_02_T1/LC08_L2SP_116034_20250426_20250429_02_T1_SR_B4.TIF" 
[5] "landsat/LC08_L2SP_116034_20250426_20250429_02_T1/LC08_L2SP_116034_20250426_20250429_02_T1_SR_B5.TIF" 
[6] "landsat/LC08_L2SP_116034_20250426_20250429_02_T1/LC08_L2SP_116034_20250426_20250429_02_T1_SR_B6.TIF" 
[7] "landsat/LC08_L2SP_116034_20250426_20250429_02_T1/LC08_L2SP_116034_20250426_20250429_02_T1_SR_B7.TIF" 
[8] "landsat/LC08_L2SP_116034_20250426_20250429_02_T1/LC08_L2SP_116034_20250426_20250429_02_T1_ST_B10.TIF"
landsat_stack <- rast(band_files)
landsat_stack
class       : SpatRaster 
dimensions  : 7961, 7841, 8  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : 172185, 407415, 4030185, 4269015  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 52N (EPSG:32652) 
sources     : LC08_L2SP_116034_20250426_20250429_02_T1_SR_B1.TIF  
              LC08_L2SP_116034_20250426_20250429_02_T1_SR_B2.TIF  
              LC08_L2SP_116034_20250426_20250429_02_T1_SR_B3.TIF  
              ... and 5 more sources
names       : LC08_~SR_B1, LC08_~SR_B2, LC08_~SR_B3, LC08_~SR_B4, LC08_~SR_B5, LC08_~SR_B6, ... 
names(landsat_stack)
[1] "LC08_L2SP_116034_20250426_20250429_02_T1_SR_B1" 
[2] "LC08_L2SP_116034_20250426_20250429_02_T1_SR_B2" 
[3] "LC08_L2SP_116034_20250426_20250429_02_T1_SR_B3" 
[4] "LC08_L2SP_116034_20250426_20250429_02_T1_SR_B4" 
[5] "LC08_L2SP_116034_20250426_20250429_02_T1_SR_B5" 
[6] "LC08_L2SP_116034_20250426_20250429_02_T1_SR_B6" 
[7] "LC08_L2SP_116034_20250426_20250429_02_T1_SR_B7" 
[8] "LC08_L2SP_116034_20250426_20250429_02_T1_ST_B10"

밴드의 이름을 다음과 같이 수정한다.

names(landsat_stack) <- c(
    "B1_Coastal", "B2_Blue", "B3_Green", "B4_Red", "B5_NIR", 
    "B6_SWIR1", "B7_SWIR2", "B10_TIR"
)

terra 패키지의 plot() 함수를 이용하여 데이터를 플로팅해본다.

plot(landsat_stack)

이제 두 가지의 핵심적인 과정이 남아 있다.

  • 우리나라의 기본 CRS(EPSG:5179)로 전환한다.

  • 영상을 해당 도시에 맞추어 잘라낸다(cropping).

여기서는 서울을 사례로 설명한다. 방금 전에 프로젝트 폴더에 정상적으로 zip 파일을 압축해제 했다면, landsat 폴더 안에 SIGUN이라는 폴더가 포함되어 있을 것이다. 여기서 서울 폴리곤만을 추출한 뒤 그 모양으로 영상을 잘라낼 것이다.

우선 SIGUN 폴더의 SIGUN.shp을 불러와 서울만을 추출한다.

SIGUN <- st_read("landsat/SIGUN/SIGUN.shp", options = "ENCODING=CP949")
seoul <- SIGUN |>  
    filter(
        SG1_NM == "서울특별시"
        )
landsat_seoul <- landsat_stack |> 
    project("EPSG:5179") |> 
    crop(seoul)

이 데이터를 저장한다.

writeRaster(landsat_seoul, filename = "landsat/landsat_seoul.tif", overwrite = TRUE)

terra 패키지의 rast() 함수를 이용하여 데이터를 불러온다. 다음부터는 앞의 과정은 생략하고 여기서부터 시작하면 된다. 한번 더 plot() 함수를 통해 시각화한다.

landsat_seoul <- rast("landsat/landsat_seoul.tif")
plot(landsat_seoul)

앞에서 실행한 잘라내기는 상하좌우 끝점을 바탕으로 한 최소포괄사각형에 기반한 것임을 알아챌 수 있을 것이다. 불규칙한 바운더리로 잘라내는 것은 불가능하며, 단지 바운더리 바깥에 위치한 픽셀을 NA로 처리할 수 있을 뿐이다. 이것을 매스킹(masking)이라고 한다. 잘라내는 것이 아니라 감추는 것이다. terra 패키지의 mask() 함수를 이용한다.

landsat_seoul_mask <- landsat_seoul |> 
    mask(seoul)
plot(landsat_seoul_mask)

2 원격탐사 데이터 시각화하기

2.1 영상 매핑(image mapping)

영상 매핑은 밴드에 포함된 정보를 시각 변수를 활용하여 지도 형태로 표현하는 과정을 의미한다. 영상 매핑에서 가장 중요한 것은 밴드 조합(band composition)인데, 3개 밴드의 RGB 조합을 통해 새로운 정보를 생성하고 시각화하는 과정을 의미한다. 여기서는 아래의 이미지에 있는 첫번째, Natural Color에 대해서만 실습한다.

출처: https://www.usgs.gov/media/images/common-landsat-band-combinations

R(B4_Red) + G(B3_Green) + B(B2_Blue) 조합의 영상을 생성한다. 현실 세계와 매우 비슷한 자연색상을 표현한다. terra 패키지의 plotRGB() 함수를 활용한다.

plotRGB(
    landsat_seoul_mask,
    r = 4, g = 3, b = 2, stretch = "lin"
)

stretch

plotRGB() 의 stretch 옵션은 대비 확장의 방법을 결정한다. 사용가능한 옵션에는 "lin"과 "hist"가 있다.

  • stretch = "lin" : 최소-최대(선형) 대비 확장

    • 전체 픽셀 값의 하위 2%와 상위 2%에 해당하는 극단값은 무시하고 선형 스트레치 수행

  • stretch = "hist" : 히스토그램 균등화

    • 사용자가 설정한 DN에 대략 동일한 수의 픽셀이 할당되게 하는 방법

혹은 ggplot2와 tidyterra 패키지를 활용하여 조금 더 아름다운 지도를 생성할 수 있다. geom_spatraster_rgb() 함수가 tidyterra 패키지의 함수이다.

ggplot() +
    geom_spatraster_rgb(
        data = landsat_seoul_mask, maxcell = 2000000,
        r = 4, g = 3, b = 2, stretch = "lin"
    ) +
    labs(
        title = "Seoul Color-Infrared Image",
        subtitle = "Landsat 8 OLI Surface Reflectance, Acquired: April 26, 2025"
    ) +
    theme_void()

Note

래스터 데이터에 대한 플롯 함수들은 보통 플로팅의 속도를 높이기 위해 재샘플링(resampling)을 통해 데이터의 공간해상도를 떨어뜨린다. maxcell 아규먼트에 영상의 총픽셀수(여기서는 1,248,585개) 이상의 값을 넣으면 이것이 방지되어 원 공간해상도로 플로팅할 수 있게 된다.

2.2 영상 변환(image transformation)

2.2.1 식생 지수

영상 대수(image algebra)를 이용하여 NDVI(normalized difference vegetation index)를 산출한다. 수식은 다음과 같다. 식생은 적색광(Red)을 흡수하고, 근적외선(NIR)을 반사하기 때문에, NDVI가 높을수록 일반적으로 건강한 녹색 식생의 존재를 의미한다.

\[ NDVI=\frac{\rho_{NIR}-\rho_{Red}}{\rho_{NIR}+\rho_{Red}} \]

수식에 맞추어 NDVI 영상을 생성한다.

numerator <- landsat_seoul_mask[["B5_NIR"]] - landsat_seoul_mask[["B4_Red"]]
denominator <- landsat_seoul_mask[["B5_NIR"]] + landsat_seoul_mask[["B4_Red"]]
NDVI_2025 <- numerator / denominator

ggplot2와 tidyterra를 활용하여 지도를 제작한다.

ggplot() +
  geom_spatraster(
    data = NDVI_2025, maxcell = 2000000
  ) + 
  scale_fill_gradientn(
    colors = brewer.pal(11, "RdYlGn"), # RColorBrewer 패키지 사용 
    na.value = NA
    ) +
  labs(
    title = "NDVI(normalized difference vegetation index), 2025.04.26.", 
    fill = "Values"
  )

2.3 변화 탐지(change detection)

변화 탐지는 동일 지역을 서로 다른 시점에 촬영한 영상을 비교 분석하여, 시간에 따른 지표면의 변화를 식별하고 정량화화는 과정이다.

변화 탐지 알고리즘은 두 부류로 나뉜다.

  • 변화 유무 탐지(binary change detection)

  • 클래스 변화 탐지(“from-to” change detection)

실습에서는 변화 유무 탐지 기법 중 영상차 변화 탐지(image differencing change detection) 기법만 다룬다.

목적은 식생 피복 상의 변화를 탐지하는 것이다. 실질적으로 최근의 NDVI 영상과 과거의 NDVI 영상을 비교하여 식생 피복 변화가 현저한 지점을 찾아 나타내는 것이다. 최근의 NDVI는 앞에서 얻었고, 문제는 과거의 NDVI를 생성하는 것이다. 이를 위해 과거(최소한 30년 정도 전) 비슷한 계절에 촬영된 Landsat Level-2 데이터를 얻어야 한다. 이를 위해 EarthExplorer에서 탐색하고 데이터를 얻어낸다. Data Sets 탭에서 Landsat 4-5를 선택하고, 비교하려는 월과 일이 유사한 데이터를 선택하면 된다.

실습에서는 landsat 폴더 안에 미리 제공된 landsat_seoul_1990 파일을 사용한다. 이 파일은 앞선 landsat_seoul 파일과 동일하게 이름 수정, 투영법 변환, cropping을 1990년의 데이터에 대해 적용한 것이다.

데이터를 불러와 masking한다.

landsat_seoul_1990 <- rast("landsat/landsat_seoul_1990.tif")
landsat_seoul_mask_1990 <- landsat_seoul_1990 |> 
  mask(seoul)

두 시점의 변화를 분석하기 위해 한 가지 매우 중요한 작업을 해야 한다. 두 연도의 데이터를 비교해보자.

landsat_seoul_mask_2025 <- landsat_seoul_mask
landsat_seoul_mask_2025
class       : SpatRaster 
dimensions  : 1011, 1235, 8  (nrow, ncol, nlyr)
resolution  : 29.98595, 29.98595  (x, y)
extent      : 935035, 972067.7, 1936664, 1966980  (xmin, xmax, ymin, ymax)
coord. ref. : KGD2002 / Unified CS (EPSG:5179) 
source(s)   : memory
varname     : landsat_seoul 
names       : B1_Coastal,   B2_Blue,  B3_Green,    B4_Red,    B5_NIR,  B6_SWIR1, ... 
min values  :   2479.531,  3730.502,  7606.901,  7485.591,  7500.738,  7602.218, ... 
max values  :  28204.867, 30010.512, 34866.984, 39143.594, 43328.496, 51793.164, ... 
landsat_seoul_mask_1990
class       : SpatRaster 
dimensions  : 1011, 1235, 7  (nrow, ncol, nlyr)
resolution  : 29.98589, 29.98589  (x, y)
extent      : 935033.9, 972066.5, 1936667, 1966983  (xmin, xmax, ymin, ymax)
coord. ref. : KGD2002 / Unified CS (EPSG:5179) 
source(s)   : memory
varname     : landsat_seoul_1990 
names       :   B1_Blue,  B2_Green,    B3_Red,    B4_NIR,  B5_SWIR1,  B7_SWIR2, ... 
min values  :  7930.364,  8045.801,  8121.254,  8072.596,  7392.665,  7275.608, ... 
max values  : 60326.312, 30200.496, 29073.572, 32792.176, 61055.176, 44568.754, ...

extent와 resolution이 미세하게 차이가 난다는 점을 알 수 있다. 이것이 에러를 발생시키기 때문에 1990년 데이터를 2025년 데이터에 맞추어 조정한다. terra 패키지의 resample() 함수를 이용한다. method 아규먼트는 보간(interpolation) 방식을 결정한다. 연속형 변수에서는 bilinear 옵션을 추천하며, 해당 옵션은 3x3 cell window를 활용하여 bilinear interportation을 수행한다. (범주형 변수는 near 옵션 추천)

landsat_seoul_mask_1990 <- resample(
  landsat_seoul_mask_1990, # resampling될 영상
  landsat_seoul_mask_2025, # resampling 기준 영상
  method = "bilinear" # interpolation 방식(3*3 cell window)
)

1990년의 NDVI를 계산한다. 밴드가 다르게 지정되어야 함에 유의한다.

numerator <- landsat_seoul_mask_1990[["B4_NIR"]] - landsat_seoul_mask_1990[["B3_Red"]]
denominator <- landsat_seoul_mask_1990[["B4_NIR"]] + landsat_seoul_mask_1990[["B3_Red"]]

NDVI_1990 <- numerator / denominator

빠르게 지도를 그려본다.

plot(NDVI_1990)

둘 사이의 영상 일치도를 높이기 위해 ’영상대영상 대조 매칭(image to image contrast matching)’을 실시한다. RSToolbox 패키지의 histMatch() 함수를 활용한다.

NDVI_1990 <- histMatch(NDVI_1990, NDVI_2025)

최종적인 1990년의 NDVI 지도를 2025년의 지도와 동일한 방식으로 제작한다.

ggplot() +
  geom_spatraster(
    data = NDVI_1990, maxcell = 2000000
  ) + 
  scale_fill_gradientn(
    colors = brewer.pal(11, "RdYlGn"),
    na.value = NA
    ) +
  labs(
    title = "NDVI(normalized difference vegetation index), 1990.04.26.", 
    fill = "Values"
  )

차분을 이용해 변화를 탐지한다. NDVI 값이 높아졌다는 것은 해당 지역의 식생이 증가한 것이라고 해석할 수 있다. 국립중앙박물관의 용산 이전(2005), 올림픽 공원 조성(2002), 한강 수변 공원 확대 등을 파악할 수 있다.

NDVI_diff <- NDVI_2025 - NDVI_1990
plot(NDVI_diff)

표준편차를 계산한다. 이 값은 개별 관측값이 평균값으로부터 평균적으로 얼마만큼 떨어져 있는가를 나타낸다.

global(NDVI_diff, fun = sd, na.rm = TRUE) |> round(digits = 1)
        sd
B5_NIR 0.1

0.1을 임계값으로 하여 범주를 구분한다.

NDVI_diff_cat <- classify(
  NDVI_diff, 
  matrix(c(-Inf, -0.1, 1, # from, to, value
           -0.1, 0.1, 2, 
           0.1, Inf, 3), ncol = 3, byrow = TRUE))
NDVI_diff_cat <- as.factor(NDVI_diff_cat)

최종적인 변화 탐지 지도를 제작한다.

ggplot() +
  geom_spatraster(
    data = NDVI_diff_cat, maxcell = 2000000
  ) +
  scale_fill_manual(
    labels = c("Loss", "No Change", "Gain"),
    values = c("red", "gray80", "darkgreen"),
    na.translate = FALSE,
    na.value = NA
  ) +
  labs(
    title = "NDVI Change Detection",
    fill = "Categories"
  ) +
  theme_void()

3 인터랙티브 시각화

3.1 인터랙티브 지도 제작 - leaflet 패키지

leaflet은 웹 상의 반응형 지도 제작에 특화된 JavaScript 라이브러리다. 이 라이브러리를 R에서 쓸 수 있게 도와주는 래퍼 패키지가 leaflet 이다.

먼저 필요한 패키지를 설치하고 불러온다.

install.packages("leaflet")

원격탐사 데이터를 가지고 leaflet으로 지도를 그리기 위해서는 addRasterImage() 함수를 사용한다.

3.1.1 수치형 지도 그리기

앞서 만든 1990년과 2025년의 NDVI 변화량 지도를 leaflet으로 시각화한다.

pal_numeric <- colorNumeric(palette = "viridis",
                            values(NDVI_diff),
                            na.color = "transparent"
                            )

leaflet() |> 
  addTiles() |> 
  addRasterImage(NDVI_diff, colors = pal_numeric, opacity = 0.8) |> 
  addLegend(pal = pal_numeric,
            values = values(NDVI_diff),
            title = "NDVI Change"
            )

3.1.2 범주형 지도 그리기

앞서 만든 변화 탐지 지도를 leaflet으로 시각화한다. 우선 카테고리를 leaflet에서 사용할 수 있도록 약간 변형해야 한다.

levels(NDVI_diff_cat) <- data.frame(
  ID = 1:3,
  category = c("Loss", "No Change", "Gain")
)

pal_category <- colorFactor(
  palette = c("red", "white", "darkgreen"),
  levels = c(1, 2, 3),
  na.color = "transparent"
)

addRasterImage를 활용하여 시각화한다. labFormat 부분은 범례에 각 카테고리가 “1, 2, 3”이 아닌 “Loss, No Change, Gain”으로 나타나게 하는 역할을 한다.

leaflet() |>
  addTiles() |>
  addRasterImage(NDVI_diff_cat, colors = pal_category, opacity = 0.8) |>
  addLegend(
    pal = pal_category,
    values = values(NDVI_diff_cat),
    title = "NDVI Change",
    labFormat = function(type, cuts, p) {
      c("Loss", "No Change", "Gain")[as.numeric(cuts)]
    }
  )

3.2 인터랙티브 그래프 제작 - plotly 패키지

plotly 패키지는 데이터 시각화용 JavaScript 라이브러리로, 반응형 시각화 도구로 최근 널리 각광받고 있다. 이 라이브러리는 다양한 오픈소스 프로그래밍 언어에서 사용이 가능하며, R에서 이를 사용할 수 있도록 만든 래퍼(wrapper) 패키지가 plotly이다. ’plotly’는 원래 별도의 시각화용 문법을 가지고 있지만, ggplotly()라는 함수를 사용하면 ggplot2로 만든 그래프를 손쉽게 plotly 그래프로 바꿀 수 있다.

우선, 사전학습 자료에서 소개된 iris 데이터를 가지고 정적인 그래프를 그려보자.

가장 먼저 필요한 패키지를 로드한다.

종에 따른 꽃잎의 길이와 너비의 관계를 보여주는 ggplot 그래프를 그리고, 이를 변수로 저장한다.

p <- ggplot(data = iris, aes(x = Petal.Length, y = Petal.Width, color = Species)) +
        geom_point() +
        labs(
            title = "Relationship Between Petal Length and Width", # 그래프 제목
            x = "Petal Length (cm)", # x축 제목
            y = "Petal Width (cm)", # y축 제목
            color = "Species", # 범례 제목
        )

p

이 그래프를 반응형으로 만들고자 한다면, 단순히 ggplot 그래프를 ggplotly() 함수의 인자로 넣으면 된다.

줌(zoom), 팬(pan), 박스 선택(box select), 라소 선택(Lasso select) 등과 같은 상호작용을 할 수 있다. 상호작용 이후 원래 그래프로 돌아오려면 더블클릭을 하면 된다. 또한 그래프 상의 데이터 포인트 위에 마우스를 올리면 정보가 표시되며, 범례의 특정 종을 클릭하면 해당 종에 대한 데이터 포인트를 보이지 않게 할 수 있다.