Analisis Tutupan dan Penggunaan Lahan (TUPLAH) dengan R dan MapBiomas Indonesia

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:

Pemuatan data raster multitahun ke dalam satu objek SpatRaster multilayer.
Reproyeksi ke sistem koordinat metrik UTM Zone 50N (EPSG: 32650) untuk perhitungan luas yang akurat secara geometrik.
Pembuatan peta seri waktu 1990–2024 dengan skema warna resmi MapBiomas.
Kuantifikasi dan visualisasi perubahan luas setiap kelas tutupan lahan pada delapan titik waktu pengamatan.

Data input yang digunakan:

Tabel 1: Spesifikasi data input TUPLAH

Parameter	Nilai
File raster	`MapBiomas_LULC_1990.tif` — `MapBiomas_LULC_2024.tif`
Direktori	`D:/MapBiomas/`
Format	GeoTIFF (`.tif`)
Sistem koordinat asal	WGS 84 / Geographic (EPSG: 4326)
Sistem koordinat target	WGS 84 / UTM Zone 50N (EPSG: 32650)
Resolusi spasial	30 m × 30 m
Skema klasifikasi	MapBiomas Indonesia Koleksi 3 [@mapbiomas2024]

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

Paket	Versi Min.	Fungsi dalam Analisis
`terra`	≥ 1.7	`rast()`, `project()`, `freq()`, `res()`, `values()`, `resample()`
`ggplot2`	≥ 3.4	`ggplot()`, `geom_raster()`, `geom_line()`, `facet_wrap()`, `ggsave()`
`dplyr`	≥ 1.1	`mutate()`, `filter()`, operator pipe `%>%`
`scales`	≥ 1.2	`label_comma()` — format ribuan pada sumbu numerik

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
}

Penggunaan 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].

`pixel_id`	Nama Kelas	Warna (Hex)
3	Formasi Hutan	`#1f8d49`
5	Mangrove	`#04381d`
76	Hutan Rawa Gambut	`#2f7360`
13	Tumbuhan Non-Hutan Lainnya	`#d89f5c`
40	Sawah	`#c71585`
35	Sawit	`#9065d0`
9	Kebun Kayu	`#7a5900`
21	Pertanian Lainnya	`#ffefc3`
30	Lubang Tambang	`#9c0027`
24	Permukiman	`#d4271e`
25	Non-Vegetasi Lainnya	`#db4d4f`
31	Tambak	`#091077`
33	Sungai/Danau/Laut	`#2532e4`
27	Citra Tertutup Awan	`#ffffff`

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)

Fungsi 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)

Pemilihan 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

Zona UTM	EPSG	Cakupan Bujur	Wilayah Representatif
UTM Zone 47N	32647	96°–102°E	Sumatera bagian barat
UTM Zone 48N	32648	102°–108°E	Sumatera timur, Kalimantan barat
UTM Zone 49S	32749	108°–114°E	Kalimantan tengah, Jawa
UTM Zone 50N	32650	114°–120°E	Kalimantan timur
UTM Zone 53S	32753	144°–150°E	Papua timur

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")
}

Penanganan 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()

Kolom	Tipe Data	Deskripsi
`x`	numerik	Koordinat easting (meter, UTM 50N)
`y`	numerik	Koordinat northing (meter, UTM 50N)
`value`	numerik	Nilai piksel = `pixel_id` kelas MapBiomas
`variable`	faktor	Nama layer = tahun (1990, 1995, …, 2024)

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

Konversi 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:

\[ \text{luas\_ha} = n_{\text{piksel}} \times \frac{r_x \cdot r_y}{10{,}000} \]

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.

Penetapan 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

Aspek	Line Chart	Faceted Line Chart
Tujuan	Perbandingan luas absolut antar kelas	Tren temporal tiap kelas secara individual
Jumlah panel	Satu panel untuk semua kelas	Satu panel per kelas
Skala sumbu-y	Seragam — perbandingan absolut valid	Bebas (`free_y`) — tren relatif tampak jelas
Kelas minor	Tertindih oleh kelas dominan	Setara dengan kelas lain dalam panel
Argumen kunci	`scale_color_manual()`	`facet_wrap(scales = "free_y")`

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

File Output	Deskripsi
`TUPLAH_Peta_2019_2024.png`	Peta seri waktu 8 panel (1990–2024) dengan warna resmi MapBiomas, resolusi 300 dpi
`TUPLAH_Perubahan_Luas_Line.png`	Line chart semua kelas dalam satu panel, sumbu-y seragam, resolusi 300 dpi
`TUPLAH_Perubahan_Luas_Facet.png`	Faceted line chart satu panel per kelas, sumbu-y bebas (`free_y`), resolusi 300 dpi

14.1 Tips Pemecahan Masalah

Memori 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)

Reproyeksi lambat pada data besar

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

Nama layer hilang setelah project()

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

Warna 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().

Nilai 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. https://doi.org/10.3390/rs12172735

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

--- title: "Analisis Tutupan dan Penggunaan Lahan (TUPLAH) dengan R dan MapBiomas Indonesia" subtitle: "MapBiomas Indonesia Koleksi 3 · Reproyeksi UTM Zone 50N · Visualisasi Seri Waktu" author: - name: "Defani Arman Alfitriansyah" affiliation: "Fakultas Kehutanan dan Lingkungan, Universitas Kuningan" email: "" date: today date-format: "D MMMM YYYY" lang: id format: html: theme: cosmo toc: true toc-depth: 3 toc-title: "Daftar Isi" toc-location: left number-sections: true number-depth: 3 code-fold: false code-tools: true code-copy: true code-line-numbers: true highlight-style: github df-print: kable fig-align: center fig-cap-location: bottom tbl-cap-location: top embed-resources: true smooth-scroll: true page-layout: full execute: eval: false echo: true warning: false message: false cache: false knitr: opts_chunk: comment: "##" collapse: false --- ```{=html} <style> body { font-family: "Source Sans Pro", "Segoe UI", sans-serif; } h1, h2, h3 { font-weight: 600; } .callout { border-radius: 6px; } code { font-size: 0.88em; } .sourceCode { border-left: 3px solid #4a7fbf; } </style> ``` ------------------------------------------------------------------------ ## Pendahuluan {#sec-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:** | Parameter | Nilai | |:---|:---| | File raster | `MapBiomas_LULC_1990.tif` — `MapBiomas_LULC_2024.tif` | | Direktori | `D:/MapBiomas/` | | Format | GeoTIFF (`.tif`) | | Sistem koordinat asal | WGS 84 / Geographic (EPSG: 4326) | | Sistem koordinat target | WGS 84 / UTM Zone 50N (EPSG: 32650) | | Resolusi spasial | 30 m × 30 m | | Skema klasifikasi | MapBiomas Indonesia Koleksi 3 [@mapbiomas2024] | : Spesifikasi data input TUPLAH {#tbl-input} ------------------------------------------------------------------------ ## Prasyarat dan Pemuatan Paket {#sec-paket} Sebelum menjalankan skrip, pastikan keempat paket berikut telah terpasang pada lingkungan R (versi minimum ≥ 4.2). Instalasi cukup dilakukan satu kali per mesin. ```{r} #| label: install-packages #| eval: false # Instalasi paket — jalankan satu kali install.packages(c("terra", "ggplot2", "dplyr", "scales")) ``` Setiap sesi analisis dimulai dengan memanggil seluruh paket yang dibutuhkan: ```{r} #| label: load-packages 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. | Paket | Versi Min. | Fungsi dalam Analisis | |:-----------------------|:----------------------:|:-----------------------| | `terra` | ≥ 1.7 | `rast()`, `project()`, `freq()`, `res()`, `values()`, `resample()` | | `ggplot2` | ≥ 3.4 | `ggplot()`, `geom_raster()`, `geom_line()`, `facet_wrap()`, `ggsave()` | | `dplyr` | ≥ 1.1 | `mutate()`, `filter()`, operator pipe `%>%` | | `scales` | ≥ 1.2 | `label_comma()` — format ribuan pada sumbu numerik | : Paket R yang digunakan beserta versi minimum dan fungsinya {#tbl-packages} ------------------------------------------------------------------------ ## Fungsi Kustom `rasterdf()` {#sec-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. ```{r} #| label: func-rasterdf 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 } ``` ::: callout-important ## Penggunaan `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. ::: ------------------------------------------------------------------------ ## Legenda Klasifikasi MapBiomas Indonesia Koleksi 3 {#sec-legenda} 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]. ```{r} #| label: define-legend 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. | `pixel_id` | Nama Kelas | Warna (Hex) | |:----------:|:---------------------------|:-----------:| | 3 | Formasi Hutan | `#1f8d49` | | 5 | Mangrove | `#04381d` | | 76 | Hutan Rawa Gambut | `#2f7360` | | 13 | Tumbuhan Non-Hutan Lainnya | `#d89f5c` | | 40 | Sawah | `#c71585` | | 35 | Sawit | `#9065d0` | | 9 | Kebun Kayu | `#7a5900` | | 21 | Pertanian Lainnya | `#ffefc3` | | 30 | Lubang Tambang | `#9c0027` | | 24 | Permukiman | `#d4271e` | | 25 | Non-Vegetasi Lainnya | `#db4d4f` | | 31 | Tambak | `#091077` | | 33 | Sungai/Danau/Laut | `#2532e4` | | 27 | Citra Tertutup Awan | `#ffffff` | : Kelas tutupan dan penggunaan lahan MapBiomas Indonesia Koleksi 3. Sumber: ATBD MapBiomas Indonesia [@mapbiomas2024]. {#tbl-legenda} ------------------------------------------------------------------------ ## Pemuatan Data Raster TUPLAH {#sec-load} 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()`. ```{r} #| label: load-raster # 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) ``` ::: callout-note ## Fungsi `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: ``` r sapply(tuplah_files, function(f) crs(rast(f), describe = TRUE)$name) ``` ::: ------------------------------------------------------------------------ ## Reproyeksi ke UTM Zone 50N (EPSG: 32650) {#sec-reproyeksi} 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. ```{r} #| label: reproject # 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) ``` ::: callout-warning ## Pemilihan 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. | Zona UTM | EPSG | Cakupan Bujur | Wilayah Representatif | |:-------------|:-----:|:-------------:|:---------------------------------| | UTM Zone 47N | 32647 | 96°–102°E | Sumatera bagian barat | | UTM Zone 48N | 32648 | 102°–108°E | Sumatera timur, Kalimantan barat | | UTM Zone 49S | 32749 | 108°–114°E | Kalimantan tengah, Jawa | | UTM Zone 50N | 32650 | 114°–120°E | Kalimantan timur | | UTM Zone 53S | 32753 | 144°–150°E | Papua timur | : Referensi kode EPSG untuk zona UTM di wilayah Indonesia {#tbl-epsg} ::: ------------------------------------------------------------------------ ## Verifikasi Nilai Piksel terhadap Legenda {#sec-verifikasi} 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. ```{r} #| label: verify-pixels # 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") } ``` ::: callout-note ## Penanganan 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: ``` r tuplah_stk[tuplah_stk == 0] <- NA ``` ::: ------------------------------------------------------------------------ ## Konversi Raster ke Data Frame {#sec-konversi} Objek `SpatRaster` dikonversi ke `data.frame` menggunakan fungsi `rasterdf()` yang didefinisikan pada @sec-rasterdf. 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. ```{r} #| label: raster-to-df # 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: | Kolom | Tipe Data | Deskripsi | |:-----------|:---------:|:------------------------------------------| | `x` | numerik | Koordinat easting (meter, UTM 50N) | | `y` | numerik | Koordinat northing (meter, UTM 50N) | | `value` | numerik | Nilai piksel = `pixel_id` kelas MapBiomas | | `variable` | faktor | Nama layer = tahun (1990, 1995, …, 2024) | : Struktur `data.frame` `tuplah_df` hasil konversi fungsi `rasterdf()` {#tbl-df} ------------------------------------------------------------------------ ## Pembuatan Peta Seri Waktu 1990–2024 {#sec-peta} 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. ```{r} #| label: fig-peta #| fig-cap: "Peta seri waktu tutupan dan penggunaan lahan 1990–2024 berdasarkan klasifikasi MapBiomas Indonesia Koleksi 3." # 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 ``` ::: callout-note ## ![](images/TUPLAH_Peta_1990_2024.png)Konversi 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. ::: ------------------------------------------------------------------------ ## Perhitungan Luas Kelas dengan `freq()` {#sec-luas} 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). ```{r} #| label: calc-area # 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: $$ \text{luas\_ha} = n_{\text{piksel}} \times \frac{r_x \cdot r_y}{10{,}000} $$ 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. ::: callout-note ## Penetapan 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. ::: ------------------------------------------------------------------------ ## Visualisasi Line Chart Perubahan Luas {#sec-linechart} 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 (@sec-peta), agar pembaca dapat mengaitkan pola spasial dengan tren temporal secara intuitif. ```{r} #| label: fig-linechart #| fig-cap: "Line chart perubahan luas kelas tutupan dan penggunaan lahan 1990–2024. Warna garis konsisten dengan skema warna peta MapBiomas." 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 ``` ::: callout-note ![](images/TUPLAH_Perubahan_Luas_Line.png) ## 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. ::: ------------------------------------------------------------------------ ## Visualisasi Faceted Line Chart per Kelas {#sec-facetchart} 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. ```{r} #| label: fig-facetchart #| fig-cap: "Faceted line chart perubahan luas per kelas tutupan dan penggunaan lahan 1990–2024. Skala sumbu-y bebas (free_y) per kelas." 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 ``` ![](images/TUPLAH_Perubahan_Luas_Facet.png) Perbedaan antara kedua jenis grafik dirangkum dalam tabel berikut: | Aspek | Line Chart | Faceted Line Chart | |:-----------------------|:-----------------------|:-----------------------| | Tujuan | Perbandingan luas absolut antar kelas | Tren temporal tiap kelas secara individual | | Jumlah panel | Satu panel untuk semua kelas | Satu panel per kelas | | Skala sumbu-y | Seragam — perbandingan absolut valid | Bebas (`free_y`) — tren relatif tampak jelas | | Kelas minor | Tertindih oleh kelas dominan | Setara dengan kelas lain dalam panel | | Argumen kunci | `scale_color_manual()` | `facet_wrap(scales = "free_y")` | : Perbandingan antara *line chart* dan *faceted line chart* {#tbl-chart-comparison} ------------------------------------------------------------------------ ## Ringkasan Output {#sec-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()`. | File Output | Deskripsi | |:-----------------------------------|:-----------------------------------| | `TUPLAH_Peta_2019_2024.png` | Peta seri waktu 8 panel (1990–2024) dengan warna resmi MapBiomas, resolusi 300 dpi | | `TUPLAH_Perubahan_Luas_Line.png` | Line chart semua kelas dalam satu panel, sumbu-y seragam, resolusi 300 dpi | | `TUPLAH_Perubahan_Luas_Facet.png` | Faceted line chart satu panel per kelas, sumbu-y bebas (`free_y`), resolusi 300 dpi | : Daftar file output yang dihasilkan skrip analisis {#tbl-output} ### Tips Pemecahan Masalah {#sec-tips} ::: {.callout-tip collapse="true"} ## Memori 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. ``` r tuplah_df <- rasterdf(tuplah_stk, aggregate = 5) ``` ::: ::: {.callout-tip collapse="true"} ## Reproyeksi lambat pada data besar Lakukan reproyeksi di luar R menggunakan `gdalwarp` (GDAL) atau QGIS, simpan sebagai GeoTIFF berproyeksi UTM, kemudian muat langsung dengan `rast()`. ::: ::: {.callout-tip collapse="true"} ## Nama layer hilang setelah `project()` Tambahkan selalu `names(tuplah_stk) <- tahun_list` segera setelah `project()` karena fungsi ini mereset atribut nama layer. ::: ::: {.callout-tip collapse="true"} ## Warna 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()`. ::: ::: {.callout-tip collapse="true"} ## Nilai piksel tidak dikenali Nilai seperti `0` atau `255` umumnya merupakan piksel NoData. Tangani sebelum analisis dilanjutkan: ``` r tuplah_stk[tuplah_stk == 0] <- NA ``` ::: ------------------------------------------------------------------------ ## Informasi Sesi R {#sec-session} ```{r} #| label: session-info #| eval: true #| echo: true sessionInfo() ``` ------------------------------------------------------------------------ ## Referensi {#sec-referensi .unnumbered} ::: {#refs} :::  **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. <https://doi.org/10.3390/rs12172735> **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>