Analisis Perubahan Luasan Tutupan dan Penggunaan Lahan dengan R dan MapBiomas Indonesia

Tutorial R, terbit pertama di RPubs (2026-04-14).

Kode

MapBiomas Indonesia Koleksi 3 · Reproyeksi UTM Zone 50N · Visualisasi Seri Waktu

Penulis

Afiliasi

Defani Arman Alfitriansyah

Fakultas Kehutanan dan Lingkungan, Universitas Kuningan

Diterbitkan

14 April 2026

1 Pendahuluan¶

Dokumen ini menyajikan prosedur teknis analisis data Tutupan dan Penggunaan Lahan (TUPLAH) berbasis citra satelit secara multitemporal menggunakan bahasa pemrograman R. Data yang digunakan merupakan produk klasifikasi tutupan lahan dari platform MapBiomas Indonesia Koleksi 3, yang mencakup 14 kelas tutupan dan penggunaan lahan untuk seluruh wilayah Indonesia [mapbiomas2024].

Analisis yang dilakukan meliputi empat tahapan utama:

1. Pemuatan data raster multitahun ke dalam satu objek SpatRaster multilayer.

2. Reproyeksi ke sistem koordinat metrik UTM Zone 50N (EPSG: 32650) untuk perhitungan luas yang akurat secara geometrik.

3. Pembuatan peta seri waktu 1990–2024 dengan skema warna resmi MapBiomas.

4. Kuantifikasi dan visualisasi perubahan luas setiap kelas tutupan lahan pada delapan titik waktu pengamatan.

Data input yang digunakan:

Tabel 1: Spesifikasi data input TUPLAH

2 Prasyarat dan Pemuatan Paket¶

Sebelum menjalankan skrip, pastikan keempat paket berikut telah terpasang pada lingkungan R (versi minimum ≥ 4.2). Instalasi cukup dilakukan satu kali per mesin.

# Instalasi paket — jalankan satu kali
install.packages(c("terra", "ggplot2", "dplyr", "scales"))

Setiap sesi analisis dimulai dengan memanggil seluruh paket yang dibutuhkan:

library(terra)   # Operasi data raster: import, proyeksi, statistik
library(ggplot2) # Pembuatan peta dan grafik ilmiah
library(dplyr)   # Manipulasi data frame dengan sintaks pipe
library(scales)  # Format angka ribuan pada label sumbu grafik

Tabel berikut merangkum versi minimum dan peran masing-masing paket dalam alur kerja ini.

Tabel 2: Paket R yang digunakan beserta versi minimum dan fungsinya

3 Fungsi Kustom rasterdf()¶

Paket ggplot2 tidak dapat memproses objek SpatRaster secara langsung. Fungsi geom_raster() membutuhkan input berupa data.frame yang memuat kolom koordinat spasial x, y, nilai atribut, dan nama layer. Fungsi kustom rasterdf() mengotomatiskan proses konversi tersebut.

Fungsi ini menerima dua argumen:

• x : objek SpatRaster yang akan dikonversi.

• aggregate : faktor pengurangan resolusi spasial sebelum konversi. Nilai aggregate = 3 berarti setiap 3 × 3 piksel direpresentasikan oleh satu piksel dalam data.frame, yang penting untuk mencegah kelebihan memori pada raster berskala regional atau nasional.

rasterdf <- function(x, aggregate = 1) {
  resampleFactor <- aggregate
  inputRaster    <- x
  inCols         <- ncol(inputRaster)
  inRows         <- nrow(inputRaster)

  # Buat raster baru dengan resolusi setelah agregasi
  resampledRaster <- rast(
    ncol = (inCols / resampleFactor),
    nrow = (inRows / resampleFactor),
    crs  = crs(inputRaster)
  )

  # Sesuaikan extent dengan raster asli
  ext(resampledRaster) <- ext(inputRaster)

  # Resample dengan nearest-neighbor — wajib untuk data kategorik
  y <- resample(inputRaster, resampledRaster, method = "near")

  # Ekstrak koordinat dan nilai piksel ke dalam data frame
  coords <- xyFromCell(y, seq_len(ncell(y)))
  dat    <- stack(values(y, dataframe = TRUE))
  names(dat) <- c("value", "variable")
  dat <- cbind(coords, dat)
  dat
}

PentingPenggunaan method = "near" pada Data Kategorik

Metode near (nearest-neighbor) wajib digunakan untuk data raster kategorik seperti TUPLAH. Metode interpolasi lain (bilinear, cubic) hanya sesuai untuk data kontinu karena akan menghasilkan nilai piksel berupa interpolasi antara dua kode kelas — misalnya nilai 4.3 di antara kelas Mangrove (pixel_id = 5) dan Formasi Hutan (pixel_id = 3) — yang tidak memiliki makna semantik dalam skema klasifikasi MapBiomas.

4 Legenda Klasifikasi MapBiomas Indonesia Koleksi 3¶

MapBiomas Indonesia Koleksi 3 menggunakan sistem klasifikasi dua tingkat yang mengacu pada kerangka FAO Land Cover Classification System [DiGregorio2000]. Setiap kelas tutupan dan penggunaan lahan direpresentasikan oleh nilai piksel unik (pixel_id) pada data raster, serta memiliki kode warna heksadesimal resmi yang konsisten di seluruh produk peta MapBiomas [mapbiomas2024].

mapbiomas_legend <- data.frame(
  pixel_id = c(3, 5, 76, 13, 40, 35, 9, 21, 30, 24, 25, 31, 33, 27),
  class_id = c(
    "Formasi Hutan",            "Mangrove",
    "Hutan Rawa Gambut",        "Tumbuhan Non-Hutan Lainnya",
    "Sawah",                    "Sawit",
    "Kebun Kayu",               "Pertanian Lainnya",
    "Lubang Tambang",           "Permukiman",
    "Non-Vegetasi Lainnya",     "Tambak",
    "Sungai/Danau/Laut",        "Citra Tertutup Awan"
  ),
  hex_color = c(
    "#1f8d49", "#04381d", "#2f7360", "#d89f5c",
    "#c71585", "#9065d0", "#7a5900", "#ffefc3",
    "#9c0027", "#d4271e", "#db4d4f", "#091077",
    "#2532e4", "#ffffff"
  ),
  stringsAsFactors = FALSE
)

# Vektor bernama untuk scale_fill_manual() dan scale_color_manual()
mb_colors <- setNames(mapbiomas_legend$hex_color,
                      as.character(mapbiomas_legend$pixel_id))
mb_labels <- setNames(mapbiomas_legend$class_id,
                      as.character(mapbiomas_legend$pixel_id))

Tabel berikut merangkum seluruh 14 kelas beserta kode warna resminya.

Tabel 3: Kelas tutupan dan penggunaan lahan MapBiomas Indonesia Koleksi 3. Sumber: ATBD MapBiomas Indonesia [mapbiomas2024].

5 Pemuatan Data Raster TUPLAH¶

Delapan file GeoTIFF yang merepresentasikan tutupan dan penggunaan lahan pada titik waktu 1990, 1995, 2000, 2005, 2010, 2015, 2020, dan 2024 dimuat sekaligus ke dalam satu objek SpatRaster berlapis-ganda (multilayer) menggunakan satu pemanggilan fungsi rast().

# Definisikan path file raster — sesuaikan dengan direktori penyimpanan data
tuplah_files <- c(
  "D:/MapBiomas/MapBiomas_LULC_1990.tif",
  "D:/MapBiomas/MapBiomas_LULC_1995.tif",
  "D:/MapBiomas/MapBiomas_LULC_2000.tif",
  "D:/MapBiomas/MapBiomas_LULC_2005.tif",
  "D:/MapBiomas/MapBiomas_LULC_2010.tif",
  "D:/MapBiomas/MapBiomas_LULC_2015.tif",
  "D:/MapBiomas/MapBiomas_LULC_2020.tif",
  "D:/MapBiomas/MapBiomas_LULC_2024.tif"
)

tahun_list <- c("1990", "1995", "2000", "2005", "2010", "2015", "2020", "2024")

# Muat seluruh layer ke dalam satu objek SpatRaster multilayer
tuplah_stk <- rast(tuplah_files)

# Berikan label bermakna pada setiap layer
names(tuplah_stk) <- tahun_list

# Periksa ringkasan objek
print(tuplah_stk)

CatatanFungsi names() pada Objek SpatRaster

Fungsi names() digunakan untuk memberikan label bermakna pada setiap layer. Label ini akan muncul sebagai judul panel pada grafik facet_wrap() dan sebagai identifikasi layer pada output freq(). Sebelum analisis, verifikasi CRS seluruh file dapat dilakukan dengan:

sapply(tuplah_files, function(f) crs(rast(f), describe = TRUE)$name)

6 Reproyeksi ke UTM Zone 50N (EPSG: 32650)¶

Produk raster MapBiomas Indonesia umumnya disimpan dalam sistem koordinat geografis WGS 84 (EPSG: 4326) dengan satuan derajat. Sistem ini tidak cocok untuk perhitungan luas berbasis piksel karena panjang satu derajat bervariasi menurut garis lintang, sehingga menghasilkan estimasi luas yang tidak akurat secara geometrik.

Proyeksi UTM Zone 50N (EPSG: 32650) dipilih karena mencakup wilayah pada rentang bujur 114°E–120°E. Dalam sistem UTM, satuan koordinat adalah meter, sehingga luas piksel dapat dihitung secara langsung sebagai perkalian resolusi spasial dalam satuan meter persegi.

# Tampilkan CRS sebelum reproyeksi
cat("CRS asli:", crs(tuplah_stk, describe = TRUE)$name, "\n")

# Reproyeksi ke UTM Zone 50N
# method = "near" wajib digunakan untuk data raster kategorik
tuplah_stk <- project(
  tuplah_stk,
  y      = "EPSG:32650",
  method = "near"
)

# Kembalikan nama layer (fungsi project() mereset nama layer)
names(tuplah_stk) <- tahun_list

# Konfirmasi CRS dan metadata setelah reproyeksi
cat("CRS sesudah:", crs(tuplah_stk, describe = TRUE)$name, "\n")
print(tuplah_stk)

PeringatanPemilihan Zona UTM yang Tepat

Zona UTM yang digunakan harus sesuai dengan lokasi geografis wilayah studi. Penggunaan zona yang tidak sesuai menghasilkan distorsi geometrik dan kesalahan perhitungan luas.

Tabel 4: Referensi kode EPSG untuk zona UTM di wilayah Indonesia

7 Verifikasi Nilai Piksel terhadap Legenda¶

Sebelum melanjutkan ke tahap visualisasi dan analisis, perlu dipastikan bahwa seluruh nilai piksel unik dalam raster memiliki padanan dalam tabel legenda mapbiomas_legend. Nilai yang tidak dikenali akan direpresentasikan sebagai NA pada peta dan grafik, yang berpotensi menyesatkan interpretasi hasil.

# Kumpulkan semua nilai unik dari seluruh layer
all_vals <- unique(na.omit(as.vector(values(tuplah_stk))))
cat("Nilai unik:", sort(all_vals), "\n")

# Identifikasi nilai yang tidak ada di legenda MapBiomas
unmatched <- setdiff(all_vals, mapbiomas_legend$pixel_id)

if (length(unmatched) > 0) {
  warning("Nilai tidak dikenal: ", paste(unmatched, collapse = ", "))
} else {
  cat("Semua nilai piksel bersesuaian dengan legenda MapBiomas.\n")
}

CatatanPenanganan Nilai NoData

Fungsi setdiff(A, B) mengembalikan elemen-elemen yang ada di vektor A tetapi tidak ada di vektor B. Jika hasilnya adalah integer(0) (vektor kosong), berarti seluruh nilai piksel telah dikenali. Nilai seperti 0 atau 255 umumnya merupakan indikasi piksel NoData yang perlu ditangani sebelum analisis dilanjutkan:

tuplah_stk[tuplah_stk == 0] <- NA

8 Konversi Raster ke Data Frame¶

Objek SpatRaster dikonversi ke data.frame menggunakan fungsi rasterdf() yang didefinisikan pada Bagian 3. Argumen aggregate = 3 berarti setiap 3 × 3 piksel direpresentasikan oleh satu piksel dalam data.frame. Pengurangan resolusi ini penting untuk mencegah penggunaan memori yang berlebihan pada raster dengan cakupan wilayah yang luas.

# Konversi ke data frame dengan pengurangan resolusi 3x
# Untuk cakupan wilayah kecil : aggregate = 1
# Untuk cakupan regional/nasional : aggregate = 3 atau lebih
tuplah_df <- rasterdf(tuplah_stk, aggregate = 3)

# Tetapkan urutan panel pada facet_wrap() melalui faktor berurutan
tuplah_df$variable <- factor(tuplah_df$variable, levels = tahun_list)

# Hapus piksel NA (area di luar batas atau tertutup awan)
tuplah_df <- tuplah_df[!is.na(tuplah_df$value), ]

cat("Dimensi data frame:", nrow(tuplah_df), "baris\n")

data.frame hasil konversi tuplah_df memiliki struktur empat kolom sebagai berikut:

Tabel 5: Struktur data.frame tuplah_df hasil konversi fungsi rasterdf()

9 Pembuatan Peta Seri Waktu 1990–2024¶

Peta tutupan lahan multitemporal divisualisasikan menggunakan facet_wrap() yang membagi tampilan menjadi delapan panel — masing-masing mewakili satu titik waktu pengamatan. Fungsi geom_raster() merender piksel sebagai grid spasial kontinu, sedangkan scale_fill_manual() menetapkan warna resmi MapBiomas untuk setiap kelas.

# Ekstrak kelas yang benar-benar hadir dalam data
# (menghindari entri legenda kosong pada grafik)
kelas_ada  <- as.character(sort(unique(tuplah_df$value)))
warna_plot <- mb_colors[kelas_ada]
label_plot <- mb_labels[kelas_ada]
warna_plot <- warna_plot[!is.na(warna_plot)]
label_plot <- label_plot[!is.na(label_plot)]

p_peta <- ggplot(data = tuplah_df) +
  geom_raster(aes(x = x, y = y, fill = as.character(value))) +
  scale_fill_manual(
    name       = "Tutupan/Penggunaan\nLahan",
    values     = warna_plot,
    labels     = label_plot,
    na.translate = FALSE
  ) +
  facet_wrap(facets = vars(variable), ncol = 3) +
  coord_sf(expand = FALSE) +
  labs(
    title    = "Tutupan dan Penggunaan Lahan (TUPLAH) 1990–2024",
    subtitle = "Klasifikasi MapBiomas Indonesia Koleksi 3",
    x        = "Easting (m)",
    y        = "Northing (m)",
    caption  = "Sumber: MapBiomas Indonesia Koleksi 3 | ATBD Desember 2024"
  ) +
  theme_bw(base_size = 11) +
  theme(
    strip.text.x  = element_text(size = 11, face = "bold"),
    legend.position  = "right",
    legend.key.size  = unit(0.5, "cm"),
    legend.text      = element_text(size = 8),
    plot.title       = element_text(face = "bold", size = 13)
  )

ggsave("TUPLAH_Peta_2019_2024.png",
       plot   = p_peta,
       width  = 16,
       height = 10,
       dpi    = 300,
       bg     = "white")

p_peta

CatatanKonversi Tipe Data pada fill

Argumen fill = as.character(value) mengkonversi nilai piksel dari numerik ke karakter agar ggplot2 memperlakukannya sebagai variabel kategorik. Tanpa konversi ini, scale_fill_manual() tidak akan bekerja sebagaimana mestinya karena ggplot2 akan mengasumsikan data bersifat kontinu. Argumen na.translate = FALSE mencegah entri NA muncul sebagai kategori tersendiri dalam legenda.

10 Perhitungan Luas Kelas dengan freq()¶

Kuantifikasi luas per kelas tutupan lahan dilakukan menggunakan fungsi freq() dari paket terra, yang menghitung jumlah piksel untuk setiap nilai unik pada setiap layer dalam objek SpatRaster. Hasilnya adalah data.frame dengan kolom layer (nama tahun), value (kode kelas piksel), dan count (jumlah piksel).

# Hitung jumlah piksel per kelas per tahun pada raster UTM 50N
freq_raw <- freq(tuplah_stk, usenames = TRUE)
glimpse(freq_raw)

# Ambil resolusi piksel dalam meter (valid karena proyeksi UTM — satuan meter)
res_m   <- res(tuplah_stk)          # c(lebar_x, tinggi_y) dalam meter
luas_m2 <- res_m[1] * res_m[2]     # Luas satu piksel dalam m²

cat("Resolusi piksel:", res_m[1], "m ×", res_m[2], "m\n")
cat("Luas satu piksel:", luas_m2, "m²\n")

# Konversi jumlah piksel ke hektar dan tambahkan label kelas
tuplah_chg <- freq_raw %>%
  mutate(
    luas_ha  = count * luas_m2 / 10000,   # 1 ha = 10.000 m²
    luas_km2 = count * luas_m2 / 1e6,     # 1 km² = 1.000.000 m²
    kelas    = factor(value,
                      levels = mapbiomas_legend$pixel_id,
                      labels = mapbiomas_legend$class_id),
    tahun    = as.numeric(layer)
  ) %>%
  filter(!is.na(kelas))                    # Hapus nilai di luar legenda MapBiomas

Formula konversi luas yang digunakan adalah sebagai berikut:

di mana $r_x$ dan $r_y$ adalah resolusi piksel dalam meter pada arah $x$ dan $y$. Penggunaan fungsi res() memastikan kode tetap valid meskipun resolusi raster berbeda dari 30 meter.

CatatanPenetapan Level Faktor secara Eksplisit

Kolom kelas dibuat sebagai factor dengan level yang diurutkan berdasarkan mapbiomas_legend$pixel_id dan label dari mapbiomas_legend$class_id. Penetapan level eksplisit ini memastikan urutan tampilan kelas pada legenda grafik konsisten dan dapat direproduksi, tidak bergantung pada urutan kemunculan kelas dalam data.

11 Visualisasi Line Chart Perubahan Luas¶

Tren perubahan luas seluruh kelas tutupan lahan divisualisasikan dalam satu grafik garis menggunakan geom_line() dan geom_point(). Setiap kelas diwakili oleh satu garis berwarna, konsisten dengan warna pada peta (Bagian 9), agar pembaca dapat mengaitkan pola spasial dengan tren temporal secara intuitif.

p_line <- ggplot(data = tuplah_chg) +
  geom_line(
    aes(x = tahun, y = luas_ha, color = kelas),
    linewidth = 1.0
  ) +
  geom_point(
    aes(x = tahun, y = luas_ha, color = kelas),
    size = 2.5
  ) +
  scale_color_manual(
    name   = "Kelas Tutupan/\nPenggunaan Lahan",
    values = setNames(mapbiomas_legend$hex_color,
                      mapbiomas_legend$class_id)
  ) +
  scale_y_continuous(labels = label_comma(big.mark = ".")) +
  scale_x_continuous(breaks = as.numeric(tahun_list)) +
  labs(
    title    = "Perubahan Luas Kelas Tutupan dan Penggunaan Lahan 1990–2024",
    subtitle = "Berdasarkan klasifikasi MapBiomas Indonesia Koleksi 3",
    x        = "Tahun",
    y        = "Luas (Hektar)",
    caption  = "Sumber: MapBiomas Indonesia Koleksi 3 | ATBD Desember 2024"
  ) +
  theme_classic(base_size = 11) +
  theme(
    legend.text      = element_text(size = 8),
    legend.key.size  = unit(0.5, "cm"),
    plot.title       = element_text(face = "bold", size = 12)
  )

ggsave("TUPLAH_Perubahan_Luas_Line.png",
       plot   = p_line,
       width  = 12,
       height = 7,
       dpi    = 300,
       bg     = "white")

p_line

Catatan

12 Konstruksi Vektor Warna Bernama¶

Konstruksi setNames(mapbiomas_legend$hex_color, mapbiomas_legend$class_id) menghasilkan vektor bernama di mana nama elemennya adalah nama kelas. Hal ini diperlukan karena kolom kelas pada tuplah_chg bertipe factor berlabel nama kelas, dan ggplot2 mencocokkan warna berdasarkan nama level faktor tersebut.

13 Visualisasi Faceted Line Chart per Kelas¶

Ketika semua kelas ditampilkan dalam satu panel, kelas-kelas dengan luas kecil — seperti Mangrove, Tambak, atau Lubang Tambang — sulit terbaca karena tertindih secara visual oleh kelas dominan. Untuk mengatasi keterbatasan ini, digunakan faceted line chart di mana setiap kelas ditampilkan dalam panel tersendiri dengan argumen scales = "free_y", sehingga tren relatif kelas minor tetap dapat diamati secara jelas.

p_facet_line <- ggplot(data = tuplah_chg) +
  geom_line(aes(x = tahun, y = luas_ha), linewidth = 0.8) +
  geom_point(aes(x = tahun, y = luas_ha), size = 2) +
  facet_wrap(
    facets = vars(kelas),
    scales = "free_y",   # Skala sumbu-y berbeda per panel
    ncol   = 4
  ) +
  scale_x_continuous(breaks = as.numeric(tahun_list)) +
  scale_y_continuous(labels = label_comma(big.mark = ".")) +
  labs(
    title   = "Perubahan Luas per Kelas Tutupan dan Penggunaan Lahan 1990–2024",
    x       = "Tahun",
    y       = "Luas (Hektar)",
    caption = "Skala sumbu-y bebas (free_y) per kelas"
  ) +
  theme_bw(base_size = 10) +
  theme(
    strip.text.x  = element_text(size = 8, face = "bold"),
    axis.text.x   = element_text(angle = 45, hjust = 1, size = 7),
    plot.title    = element_text(face = "bold", size = 12)
  )

ggsave("TUPLAH_Perubahan_Luas_Facet.png",
       plot   = p_facet_line,
       width  = 16,
       height = 10,
       dpi    = 300,
       bg     = "white")

p_facet_line

Perbedaan antara kedua jenis grafik dirangkum dalam tabel berikut:

Tabel 6: Perbandingan antara line chart dan faceted line chart

14 Ringkasan Output¶

Setelah menjalankan skrip secara lengkap, tiga file PNG akan tersimpan di direktori kerja R aktif. Direktori aktif dapat diperiksa dengan getwd() dan diubah menggunakan setwd().

Tabel 7: Daftar file output yang dihasilkan skrip analisis

14.1 Tips Pemecahan Masalah¶

TipMemori tidak mencukupi saat konversi data frame

Naikkan nilai argumen aggregate pada rasterdf() (misalnya dari 3 ke 5 atau 10). Ini mengurangi jumlah baris data.frame tanpa mengubah raster asli.

tuplah_df <- rasterdf(tuplah_stk, aggregate = 5)

TipReproyeksi lambat pada data besar

Lakukan reproyeksi di luar R menggunakan gdalwarp (GDAL) atau QGIS, simpan sebagai GeoTIFF berproyeksi UTM, kemudian muat langsung dengan rast().

TipNama layer hilang setelah project()

Tambahkan selalu names(tuplah_stk) <- tahun_list segera setelah project() karena fungsi ini mereset atribut nama layer.

TipWarna tidak muncul atau salah di grafik

Pastikan kolom kelas bertipe factor dan nama level faktornya identik dengan nama dalam vektor warna yang dibuat menggunakan setNames().

TipNilai piksel tidak dikenali

Nilai seperti 0 atau 255 umumnya merupakan piksel NoData. Tangani sebelum analisis dilanjutkan:

tuplah_stk[tuplah_stk == 0] <- NA

15 Informasi Sesi R¶

sessionInfo()

## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_Indonesia.utf8  LC_CTYPE=English_Indonesia.utf8   
## [3] LC_MONETARY=English_Indonesia.utf8 LC_NUMERIC=C                      
## [5] LC_TIME=English_Indonesia.utf8    
## 
## time zone: Asia/Jakarta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] htmlwidgets_1.6.4 compiler_4.5.2    fastmap_1.2.0     cli_3.6.5        
##  [5] tools_4.5.2       htmltools_0.5.9   otel_0.2.0        rstudioapi_0.17.1
##  [9] yaml_2.3.12       rmarkdown_2.29    knitr_1.51        jsonlite_2.0.0   
## [13] xfun_0.56         digest_0.6.39     rlang_1.1.7       evaluate_1.0.5

Referensi¶

MapBiomas Indonesia. (2024). Algorithm Theoretical Basis Document (Dokumen Landasan Teoritis Algoritma), Koleksi 3, Versi 1. Jakarta: Auriga Nusantara. Desember 2024.

Di Gregorio, A. & Jansen, L.J.M. (2000). Land Cover Classification System (LCCS): Classification Concepts and User Manual for Software Version 1.0. Food and Agriculture Organization of the United Nations (FAO), Roma.

Hijmans, R.J. (2024). terra: Spatial Data Analysis. R package version 1.7. https://CRAN.R-project.org/package=terra

Souza, C.M. Jr. et al. (2020). Reconstructing Three Decades of Land Use and Land Cover Changes in Brazilian Biomes with Landsat Archive and Earth Engine. Remote Sensing, 12(17), 2735. Souza et al. (2020)

Walker, J.D. (2024). Geographic Data Science with R: Visualizing and Analyzing Environmental Change. Edisi Pertama. CRC Press / Taylor & Francis Group.

Wickham, H. et al. (2023). ggplot2: Elegant Graphics for Data Analysis. Edisi Ketiga. Springer. https://ggplot2-book.org

Wickham, H. et al. (2023). dplyr: A Grammar of Data Manipulation. R package version 1.1. https://CRAN.R-project.org/package=dplyr