Dokumen ini berisi proses pengolahan citra satelit Landsat untuk menghitung Normalized Difference Vegetation Index (NDVI). NDVI adalah indeks yang digunakan untuk mengidentifikasi vegetasi berdasarkan perbedaan reflektansi cahaya Near-Infrared (NIR) dan Red.

NDVI dihitung menggunakan rumus:

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

Dimana:

Near-Infrared (NIR) = Band 5
Red = Band 4

library(terra)
library(tidyverse)

rootDir <- "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/"
lst_files <- list.files(rootDir, pattern = "_B[0-9]+\\.TIF$", full.names = TRUE)
print(lst_files)

##  [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B1.TIF" 
##  [2] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B10.TIF"
##  [3] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B11.TIF"
##  [4] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B2.TIF" 
##  [5] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B3.TIF" 
##  [6] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B4.TIF" 
##  [7] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B5.TIF" 
##  [8] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B6.TIF" 
##  [9] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B7.TIF" 
## [10] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B8.TIF" 
## [11] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B9.TIF"

Kode ini digunakan untuk mencari file raster dalam format .TIF yang berada dalam folder.

Fungsi list.files() dengan pola “_B[0-9]+\.TIF$” memastikan hanya file yang memiliki band numerik (1 hingga 11) yang dipilih.

Perintah print(lst_files) mencetak daftar file yang ditemukan untuk memastikan data telah dimuat dengan benar.

1 Periksa Cakupan Extent

for (file in lst_files) {
  r <- rast(file)
  print(file)
  print(ext(r))  # Menampilkan extent tiap file
}

## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B1.TIF"
## SpatExtent : 218685, 446415, -276315, -44385 (xmin, xmax, ymin, ymax)
## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B10.TIF"
## SpatExtent : 218685, 446415, -276315, -44385 (xmin, xmax, ymin, ymax)
## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B11.TIF"
## SpatExtent : 218685, 446415, -276315, -44385 (xmin, xmax, ymin, ymax)
## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B2.TIF"
## SpatExtent : 218685, 446415, -276315, -44385 (xmin, xmax, ymin, ymax)
## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B3.TIF"
## SpatExtent : 218685, 446415, -276315, -44385 (xmin, xmax, ymin, ymax)
## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B4.TIF"
## SpatExtent : 218685, 446415, -276315, -44385 (xmin, xmax, ymin, ymax)
## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B5.TIF"
## SpatExtent : 218685, 446415, -276315, -44385 (xmin, xmax, ymin, ymax)
## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B6.TIF"
## SpatExtent : 218685, 446415, -276315, -44385 (xmin, xmax, ymin, ymax)
## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B7.TIF"
## SpatExtent : 218685, 446415, -276315, -44385 (xmin, xmax, ymin, ymax)
## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B8.TIF"
## SpatExtent : 218692.5, 446407.5, -276307.5, -44392.5 (xmin, xmax, ymin, ymax)
## [1] "D:/2. Pengembangan diri/1 Exercise Bagus/Remote Sensing in R/Citra Jambi/LC09_L1TP_125061_20240703_20240703_02_T1_B9.TIF"
## SpatExtent : 218685, 446415, -276315, -44385 (xmin, xmax, ymin, ymax)

Dalam langkah ini, setiap file raster dibaca menggunakan rast(file), kemudian cakupan spasial (extent) dari masing-masing band diperiksa menggunakan ext(r). Semua band harus memiliki extent yang sama agar bisa digabungkan dalam satu stack raster. Namun, jika ada perbedaan, seperti pada Band 8, hal ini bisa menyebabkan error saat proses stacking. ¹.

Resampling Band

ref_raster <- rast(lst_files[1])

# Resample semua band agar sesuai dengan referensi
aligned_rasters <- lapply(lst_files, function(f) {
  r <- rast(f)
  resample(r, ref_raster, method = "bilinear")
})

## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================

Untuk mengatasi perbedaan extent, resampling dilakukan dengan memilih Band 1 sebagai referensi (ref_raster <- rast(lst_files[1])).

Fungsi resample() digunakan untuk menyesuaikan setiap band dengan resolusi dan extent referensi.

Metode “bilinear” digunakan karena data raster ini bersifat kontinu. ².

2 Stack Raster

lst_images <- rast(aligned_rasters)

Setelah semua band disesuaikan, lalu digabung menjadi satu stack raster menggunakan rast(aligned_rasters). Objek lst_images sekarang berisi seluruh band yang telah sejajar.

Menampilkan Citra Raster

print(lst_images)  # Mengecek hasil akhir

## class       : SpatRaster 
## dimensions  : 7731, 7591, 11  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 218685, 446415, -276315, -44385  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 48N (EPSG:32648) 
## sources     : spat_56b83535432b_22200_o8up7yuDdJ6XHBX.tif  
##               spat_56b83fe77509_22200_zxzxx9jyyx6UnnI.tif  
##               spat_56b836106a61_22200_80Q6udnsxVen2Mq.tif  
##               ... and 8 more sources
## varnames    : LC09_L1TP_125061_20240703_20240703_02_T1_B1 
##               LC09_L1TP_125061_20240703_20240703_02_T1_B1 
##               LC09_L1TP_125061_20240703_20240703_02_T1_B1 
##               ...
## names       : LC09_~T1_B1, LC09_~1_B10, LC09_~1_B11, LC09_~T1_B2, LC09_~T1_B3, LC09_~T1_B4, ... 
## min values  :        9195,       18512,       18183,        8217,        7050,        6165, ... 
## max values  :       47775,       26099,       24272,       50373,       51746,       54803, ...

plot(lst_images)   # Menampilkan citra hasil stack

Setelah stack raster dibuat, hasilnya diperiksa dengan print(lst_images), yang akan menampilkan informasi tentang jumlah band, resolusi, dan cakupan spasial. Kemudian, plot raster menggunakan plot(lst_images) untuk memastikan bahwa data telah dimuat dengan benar.

Citra RGB

names(lst_images)

##  [1] "LC09_L1TP_125061_20240703_20240703_02_T1_B1" 
##  [2] "LC09_L1TP_125061_20240703_20240703_02_T1_B10"
##  [3] "LC09_L1TP_125061_20240703_20240703_02_T1_B11"
##  [4] "LC09_L1TP_125061_20240703_20240703_02_T1_B2" 
##  [5] "LC09_L1TP_125061_20240703_20240703_02_T1_B3" 
##  [6] "LC09_L1TP_125061_20240703_20240703_02_T1_B4" 
##  [7] "LC09_L1TP_125061_20240703_20240703_02_T1_B5" 
##  [8] "LC09_L1TP_125061_20240703_20240703_02_T1_B6" 
##  [9] "LC09_L1TP_125061_20240703_20240703_02_T1_B7" 
## [10] "LC09_L1TP_125061_20240703_20240703_02_T1_B8" 
## [11] "LC09_L1TP_125061_20240703_20240703_02_T1_B9"

names(lst_images) <- paste0("Band_", 1:nlyr(lst_images))
lst_rgb <- c(lst_images[["Band_7"]], lst_images[["Band_5"]], lst_images[["Band_3"]])
plotRGB(lst_rgb, main="RGB", stretch="lin")

Setiap band dalam stack diberi nama baru menggunakan paste0(“Band_”, 1:nlyr(lst_images)), sehingga memudahkan akses dalam analisis selanjutnya.

Perintah plotRGB() digunakan untuk menampilkan citra RGB dengan penyesuaian kontras linear (stretch=“lin”), yang dapat meningkatkan visualisasi.

Simpan file stack raster

writeRaster(lst_images, str_c(rootDir,"stacked_raster_Jambi.TIF"), overwrite=TRUE)

## |---------|---------|---------|---------|=========================================

3 Hitung NDVI

landsat_stack_csf <- lst_images
landsat_ndvi <- (landsat_stack_csf[[5]] - landsat_stack_csf[[4]]) / (landsat_stack_csf[[5]] + landsat_stack_csf[[4]])
plot(landsat_ndvi)

Hasil penghitungan NDVI diplot menggunakan plot(landsat_ndvi), yang akan menampilkan distribusi indeks vegetasi. ³.

Nilai Statistik NDVI yang dihasilkan

# Mengambil semua nilai NDVI tanpa NA
ndvi_values <- na.omit(values(landsat_ndvi))

# minimum dan maksimum
min_ndvi <- min(ndvi_values)
max_ndvi <- max(ndvi_values)

# kuartil
ndvi_quantiles <- quantile(ndvi_values, probs = c(0.25, 0.5, 0.75), na.rm = FALSE)

# Cetak nilai statistik dasar
print(paste("Minimum NDVI:", min_ndvi))

## [1] "Minimum NDVI: -0.201712902625209"

print(paste("Maksimum NDVI:", max_ndvi))

## [1] "Maksimum NDVI: 0.172377084729724"

print("Kuartil NDVI:ndvi_quantiles")

## [1] "Kuartil NDVI:ndvi_quantiles"

# Menampilkan range histogram untuk melihat apakah mencakup nilai minimum
hist_breaks <- hist(ndvi_values, plot = FALSE)$breaks
print(paste("Range bin histogram: ", min(hist_breaks), "sampai", max(hist_breaks)))

## [1] "Range bin histogram:  -0.21 sampai 0.18"

4 References

Kaitlyn. (n.d.). Calculating NDVI. RPubs. Diakses pada 16 Februari 2025, dari https://rpubs.com/kmp24/280687

Waswa, R. (n.d.). Stacking Landsat 9 Images. RPubs. Diakses pada 16 Februari 2025, dari http://rpubs.com/roywaswa/stacking_lst_9

Direktorat Statistik Kesejahteraan Rakyat, BPS, saptahas@bps.go.id

Dari hasil yang ditampilkan, semua band memiliki extent yang sama kecuali Band 8, yang memiliki perbedaan kecil di batas xmin, xmax, ymin, dan ymax. Ini yang menyebabkan error ketika mencoba membuat stack raster↩︎
Jika data berupa kategori (misalnya kelas penggunaan lahan), metode “near” lebih sesuai↩︎
NDVI bernilai mendekati 1 menunjukkan vegetasi sehat, sedangkan nilai mendekati -1 menunjukkan daerah tanpa vegetasi↩︎

Stacking Raster & Calculating NDVI

Kota Jambi (Akuisisi Citra Juni-Desember 2024)

Sapta Hastho Ponco

2025-02-20

1 Periksa Cakupan Extent

2 Stack Raster

3 Hitung NDVI

4 References