Pendahuluan

Analisis ini bertujuan untuk mencari desa di Kalimantan Barat yang memiliki karakteristik serupa (desa kontrol) dengan desa yang menerima program/intervensi (desa intervensi), desa kontrol digunakan sebagai pembanding terhadap desa intervensi. Ini memungkinkan kita menjawab pertanyaan penting: Apa yang terjadi di desa intervensi jika program tidak dilakukan?

Dengan menggunakan pendekatan Dynamic Time Warping (DTW) memungkinkan pencocokan pola perubahan vegetasi secara historis, sekaligus mempertimbangkan variabel tambahan, seperti:

Dengan demikian, desa kontrol yang dipilih tidak hanya mirip dari segi vegetasi, tetapi juga dari konteks ekologi dan pengelolaan wilayah.

Dataset yang Digunakan

  • NDVI (2001–2015): Indikator kepadatan dan kesehatan vegetasi, menggambarkan perubahan lingkungan desa dari waktu ke waktu.
  • Status Kawasan Konservasi (SK.733/Menhut-II/2014): Menandai apakah desa berada dalam area konservasi, yang mempengaruhi regulasi dan akses sumber daya.
  • Keberadaan Mangrove: Sebagai salah satu elemen penting dalam ekosistem pesisir.

Kenapa Menggunakan NDVI?

NDVI (Normalized Difference Vegetation Index) dipilih karena:

  • Mudah diperoleh secara konsisten dari data satelit jangka panjang.
  • Merepresentasikan dinamika vegetasi, yang mencerminkan dampak aktivitas manusia dan perubahan lingkungan.
  • Sensitif terhadap degradasi atau pemulihan lahan, seperti deforestasi, pertanian, dan restorasi.
  • Kuat sebagai indikator proxy untuk memahami kondisi ekologis desa, terutama sebelum dan sesudah intervensi.

Dengan NDVI, kita dapat menilai bagaimana pola vegetasi berubah secara historis, dan menggunakannya untuk membandingkan desa yang memiliki “sejarah ekologis” serupa.

Proses Analisis

  1. Eksplorasi & Visualisasi NDVI
    • Pemeriksaan data hilang
    • Statistik deskriptif (rata-rata, median, dll.)
    • Histogram sebaran NDVI secara keseluruhan
    • Boxplot variasi NDVI antar tahun
  2. Penerapan DTW
    • Mencocokkan pola NDVI desa intervensi dengan seluruh desa lainnya
    • Mengukur kemiripan berdasarkan bentuk tren, bukan hanya nilai mutlak
  3. Evaluasi Hasil Pencocokan
    • Menggunakan metrik: DTW Distance dan RMSE (Root Mean Square Error)
    • Memastikan desa kontrol benar-benar memiliki pola vegetasi yang mirip sebelum intervensi
  4. Visualisasi Spasial
    • Menampilkan hasil pencocokan desa intervensi-kontrol dalam bentuk interaktif
    • Menggunakan leaflet dan visNetwork

Mengapa DTW?

Untuk membandingkan dua desa, kita tidak cukup hanya melihat angka rata-rata NDVI. Kita perlu memahami bagaimana pola perubahan NDVI-nya dari tahun ke tahun.

Dynamic Time Warping (DTW) adalah metode yang sangat berguna dalam hal ini karena:

  • Bisa mencocokkan pola NDVI yang naik-turun secara mirip, walaupun terjadi di waktu yang sedikit bergeser
  • Tidak terpengaruh perbedaan panjang data atau noise kecil
  • Membantu mengidentifikasi desa yang benar-benar mirip secara tren historis, bukan hanya sesaat

Dengan DTW, kita bisa mendapatkan desa kontrol yang adil untuk membandingkan dampak program di desa intervensi, sehingga evaluasi menjadi lebih akurat dan terpercaya.

Library yang Digunakan

Analisis ini menggunakan beberapa pustaka R yang memiliki fungsi berbeda-beda, di antaranya:
1. dplyr – Manipulasi data seperti filter, select, mutate, join, dan summarise.
2. ggplot2 – Visualisasi data (grafik histogram, boxplot, scatter plot, dll.).
3. naniar – Analisis dan visualisasi data yang hilang (missing values).
4. VIM – Visualisasi dan imputasi data hilang (contohnya grafik aggr).
5. tidyr – Merapikan dan mentransformasi data (fungsi pivot_longer, pivot_wider).
6. readr – Membaca dan menulis data CSV/TXT dengan cepat dan aman.
7. dtw – Dynamic Time Warping untuk pencarian kesamaan pola pada data time-series.
8. purrr – Pemrograman fungsional untuk operasi list (misalnya map, map2).
9. stringr – Manipulasi string (replace, detect, substring, regex).
10. sf – Analisis dan visualisasi data spasial (shapefile, GeoJSON).
11. Metrics – Evaluasi model/statistik error (RMSE, MAE, dll.).
12. visNetwork – Visualisasi jaringan interaktif (nodes-edges) berbasis JavaScript.
13. leaflet dan leaflet.extras – Peta interaktif berbasis leaflet.js untuk R.

library(dplyr)            # Manipulasi data (filter, select, mutate, join, summarise)
## Warning: package 'dplyr' was built under R version 4.5.1
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)          # Visualisasi data (grafik histogram, boxplot, scatter plot, dll.)
## Warning: package 'ggplot2' was built under R version 4.5.1
library(naniar)           # Analisis dan visualisasi data yang hilang (missing values)
library(VIM)              # Visualisasi dan imputasi data hilang (misalnya grafik aggr)
## Loading required package: colorspace
## Loading required package: grid
## VIM is ready to use.
## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues
## 
## Attaching package: 'VIM'
## The following object is masked from 'package:datasets':
## 
##     sleep
library(tidyr)            # Merapikan/transformasi data (pivot_longer, pivot_wider)
## Warning: package 'tidyr' was built under R version 4.5.1
library(readr)            # Membaca/menulis data CSV/TXT dengan cepat dan aman
library(dtw)              # Dynamic Time Warping (mencari kesamaan pola time-series)
## Loading required package: proxy
## 
## Attaching package: 'proxy'
## The following objects are masked from 'package:stats':
## 
##     as.dist, dist
## The following object is masked from 'package:base':
## 
##     as.matrix
## Loaded dtw v1.23-1. See ?dtw for help, citation("dtw") for use in publication.
library(purrr)            # Pemrograman fungsional (map, map2, peta list menjadi data frame)
library(stringr)          # Manipulasi string (replace, detect, substring, regex)
library(sf)               # Analisis dan visualisasi data spasial (shapefile, GeoJSON)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(Metrics)          # Evaluasi model/statistik error (RMSE, MAE, dll.)
library(visNetwork)       # Visualisasi jaringan interaktif (nodes-edges) berbasis JavaScript
library(leaflet)          # Peta interaktif berbasis leaflet.js untuk R
library(leaflet.extras)   # Tambahkan ini di awal kalau belum
## Warning: package 'leaflet.extras' was built under R version 4.5.1

Bagian A: Eksplorasi & Visualisasi NDVI

Bagian ini bertujuan untuk memahami karakteristik awal data NDVI yang digunakan sebagai indikator kepadatan dan kesehatan vegetasi desa.

Langkah-langkah eksplorasi data mencakup:

langkah ini dilakukan untuk memastikan bahwa data NDVI memiliki kualitas yang baik sebelum diterapkan dalam proses pencocokan desa menggunakan metode Dynamic Time Warping (DTW).

Langkah awal pada analisis ini adalah memeriksa kelengkapan data pada variabel NDVI.
Kode berikut menghitung jumlah nilai yang hilang (NA) pada kolom mean_ndvi:

# 1. Impor & Persiapan Data NDVI
ndvi <- read.csv("input/bahan_ndvi_2001-2015.csv") %>%
  mutate(id_desa = paste(kecamatan, desa, sep = "_"))

# Cek isi data
head(ndvi)
##   tahun_data kecamatan desa mean_ndvi intervensi tahun_intervensi       id_desa
## 1       2001  sebangki agak 0.7036637      tidak                0 sebangki_agak
## 2       2002  sebangki agak 0.7404532      tidak                0 sebangki_agak
## 3       2003  sebangki agak 0.7448543      tidak                0 sebangki_agak
## 4       2004  sebangki agak 0.7334519      tidak                0 sebangki_agak
## 5       2005  sebangki agak 0.7410669      tidak                0 sebangki_agak
## 6       2006  sebangki agak 0.7081221      tidak                0 sebangki_agak
# Pastikan mean_ndvi numerik
ndvi <- ndvi %>%
  mutate(mean_ndvi = as.numeric(as.character(mean_ndvi)))

# 2. Statistik Deskriptif
# Cek jumlah missing value
ndvi %>%
  summarise(na_ndvi = sum(is.na(mean_ndvi)))
##   na_ndvi
## 1       0

Hasil yang diperoleh adalah 0, yang berarti tidak ada data hilang pada variabel mean_ndvi. Hal ini menunjukkan bahwa seluruh desa dan periode waktu memiliki nilai NDVI yang lengkap, sehingga data siap digunakan untuk analisis lebih lanjut.

# Statistik deskriptif NDVI
tabel_ndvi <- ndvi %>%
  summarise(
    rata_ndvi   = mean(mean_ndvi, na.rm = TRUE),
    sd_ndvi     = sd(mean_ndvi, na.rm = TRUE),
    min_ndvi    = min(mean_ndvi, na.rm = TRUE),
    max_ndvi    = max(mean_ndvi, na.rm = TRUE),
    median_ndvi = median(mean_ndvi, na.rm = TRUE)
  )

print(tabel_ndvi)
##   rata_ndvi    sd_ndvi  min_ndvi  max_ndvi median_ndvi
## 1 0.6798429 0.05745434 0.3500694 0.8159224   0.6853771

Interpretasi ini memberikan gambaran umum mengenai kondisi tutupan vegetasi pada lokasi dan periode yang dianalisis, serta dapat membantu dalam proses pemilihan desa kontrol maupun evaluasi dampak intervensi.

# 3. Visualisasi NDVI
# Histogram NDVI
ggplot(ndvi, aes(x = mean_ndvi)) +
  geom_histogram(bins = 30, fill = "lightgreen", color = "black") +
  ggtitle("Distribusi NDVI") +
  xlab("NDVI") +
  theme_minimal()

Histogram di atas menggambarkan distribusi nilai rata-rata NDVI (Normalized Difference Vegetation Index) pada seluruh unit observasi dalam dataset. Mayoritas nilai NDVI berada pada rentang 0.6 hingga 0.8, yang menunjukkan tingkat kehijauan atau tutupan vegetasi yang relatif tinggi. Pola distribusi terlihat cukup simetris dan tidak mencerminkan adanya outlier ekstrem.

Hal ini sejalan dengan nilai rata-rata NDVI sebesar 0.68 dan median 0.685 yang sebelumnya telah dihitung, mengindikasikan bahwa sebagian besar area memiliki kondisi vegetasi yang baik dan relatif seragam. Sebagian kecil data memiliki NDVI di bawah 0.4, yang dapat mengindikasikan area terbuka, terganggu, atau memiliki tutupan vegetasi rendah.

Distribusi ini memberikan gambaran umum yang penting untuk memahami karakteristik ekologis wilayah studi, serta dapat digunakan sebagai dasar awal dalam membandingkan perubahan tutupan vegetasi antar waktu atau antar kelompok desa (misalnya antara desa intervensi dan desa kontrol).

# Boxplot NDVI per Tahun
ggplot(ndvi, aes(x = as.factor(tahun_data), y = mean_ndvi)) +
  geom_boxplot(fill = "darkseagreen", outlier.color = "gray") +
  ggtitle("Distribusi NDVI per Tahun") +
  xlab("Tahun") +
  ylab("NDVI") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90))

Boxplot di atas menunjukkan sebaran nilai NDVI rata-rata untuk setiap tahun pengamatan. Setiap kotak merepresentasikan distribusi kuartil NDVI pada tahun tertentu, dengan garis tengah menunjukkan nilai median.

Secara umum, pola NDVI dari tahun ke tahun terlihat cukup stabil, dengan nilai median yang konsisten berada di kisaran 0.65 hingga 0.75. Hal ini menunjukkan bahwa tutupan vegetasi cenderung tetap dan tidak mengalami fluktuasi besar selama periode pengamatan. Rentang antar-kuartil (IQR) yang relatif sempit di sebagian besar tahun juga mengindikasikan bahwa variasi tutupan vegetasi antar lokasi dalam tahun yang sama cukup kecil.

Namun, terdapat beberapa outlier di beberapa tahun, yang ditunjukkan oleh titik-titik di luar rentang whisker. Outlier ini kemungkinan berasal dari lokasi atau waktu dengan perubahan vegetasi ekstrem, seperti pembukaan lahan, kebakaran, atau konversi tutupan.

Visualisasi ini penting untuk memahami dinamika vegetasi tahunan secara umum, serta dapat membantu mengidentifikasi tahun-tahun yang perlu ditelusuri lebih lanjut terkait perubahan tutupan lahan, terutama jika akan dilakukan evaluasi dampak program berbasis waktu.

Bagian B: Data Matching Desa Kontrol dengan DTW

Pada bagian ini dilakukan proses pencarian desa kontrol yang memiliki karakteristik serupa dengan desa intervensi.
Metode Dynamic Time Warping (DTW) digunakan untuk mencocokkan pola deret waktu NDVI antar desa, sehingga dapat ditemukan desa dengan dinamika vegetasi yang paling mirip.

Selain kesamaan pola NDVI, proses pencocokan juga mempertimbangkan beberapa kriteria tambahan berikut:

Hasil dari tahap ini berupa daftar pasangan desa intervensi dan tiga desa kontrol terbaik untuk masing-masing, yang dipilih berdasarkan nilai DTW terkecil.

# DYNAMIC TIME WARPING MATCHING UNTUK DESA KONTROL
# 1. NDVI (time-series) + Mangrove + Status Kawasan Konservasi (status_kk)

# 2. Load Data
ndvi <- read_csv("input/bahan_ndvi_2001-2015.csv") %>%
  mutate(id_desa = str_to_lower(
    str_replace_all(paste(kecamatan, desa, sep = "_"), "\\s+", "_")
  ))
## Rows: 9960 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): kecamatan, desa, intervensi
## dbl (3): tahun_data, mean_ndvi, tahun_intervensi
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
mangrove_kk <- read_csv("input/mangrove_kk.csv")  # Info mangrove dan status_kk
## Rows: 664 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): id_desa, mangrove, status_kk
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
desa_coords <- st_read("input/desa_dalam_kawasan_konservasi.shp") %>%
  st_transform(4326)
## Reading layer `desa_dalam_kawasan_konservasi' from data source 
##   `D:\MEL\proyek_desa_kontrol\input\desa_dalam_kawasan_konservasi.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 664 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.6773 ymin: -3.067918 xmax: 114.2054 ymax: 2.081301
## Geodetic CRS:  WGS 84
fungsi_kawasan <- st_read("input/kk_kalbar.shp") %>%
  st_transform(4326) %>%
  select(kk_in)
## Reading layer `kk_kalbar' from data source 
##   `D:\MEL\proyek_desa_kontrol\input\kk_kalbar.shp' using driver `ESRI Shapefile'
## Simple feature collection with 861 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.6506 ymin: -2.995644 xmax: 114.2059 ymax: 2.082183
## Geodetic CRS:  WGS 84
head(ndvi)
## # A tibble: 6 × 7
##   tahun_data kecamatan desa  mean_ndvi intervensi tahun_intervensi id_desa      
##        <dbl> <chr>     <chr>     <dbl> <chr>                 <dbl> <chr>        
## 1       2001 sebangki  agak      0.704 tidak                     0 sebangki_agak
## 2       2002 sebangki  agak      0.740 tidak                     0 sebangki_agak
## 3       2003 sebangki  agak      0.745 tidak                     0 sebangki_agak
## 4       2004 sebangki  agak      0.733 tidak                     0 sebangki_agak
## 5       2005 sebangki  agak      0.741 tidak                     0 sebangki_agak
## 6       2006 sebangki  agak      0.708 tidak                     0 sebangki_agak
head(mangrove_kk)
## # A tibble: 6 × 3
##   id_desa                     mangrove status_kk
##   <chr>                       <chr>    <chr>    
## 1 hulu_sungai_beginci_darat   tidak    HL       
## 2 seluas_bengkawan            tidak    CA       
## 3 air_besar_bentiang_madomang tidak    CA       
## 4 kepulauan_karimata_betok    ya       CA       
## 5 suti_semarang_cempaka_putih tidak    CA       
## 6 kubu_dabong                 ya       HL
head(desa_coords)
## Simple feature collection with 6 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.6773 ymin: -1.727216 xmax: 111.2622 ymax: 1.257432
## Geodetic CRS:  WGS 84
##      kabupaten          kecamatan              desa intervensi tahun_inte
## 1     ketapang        hulu_sungai     beginci_darat         ya       2022
## 2   bengkayang             seluas         bengkawan         ya       2017
## 3       landak          air_besar bentiang_madomang         ya       2019
## 4 kayong_utara kepulauan_karimata             betok         ya       2022
## 5   bengkayang      suti_semarang     cempaka_putih         ya       2021
## 6    kubu_raya               kubu            dabong         ya       2020
##      luas_ha mangrove status_kk                     id_desa status_asl
## 1 102360.698    tidak        HL   hulu_sungai_beginci_darat         HL
## 2  17082.620    tidak        CA            seluas_bengkawan         CA
## 3  11523.536    tidak        CA air_besar_bentiang_madomang         CA
## 4  12306.810       ya        CA    kepulauan_karimata_betok         CA
## 5   3055.406    tidak        CA suti_semarang_cempaka_putih         CA
## 6   9587.332       ya        HL                 kubu_dabong         HL
##                         geometry
## 1 MULTIPOLYGON (((111.1533 -1...
## 2 MULTIPOLYGON (((109.9089 1....
## 3 MULTIPOLYGON (((109.9769 0....
## 4 MULTIPOLYGON (((108.939 -1....
## 5 MULTIPOLYGON (((109.7976 0....
## 6 MULTIPOLYGON (((109.3084 -0...
head(fungsi_kawasan)
## Simple feature collection with 6 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 110.7622 ymin: -2.86578 xmax: 110.8759 ymax: -2.672358
## Geodetic CRS:  WGS 84
##   kk_in                       geometry
## 1 tidak MULTIPOLYGON (((110.8052 -2...
## 2 tidak MULTIPOLYGON (((110.8699 -2...
## 3 tidak MULTIPOLYGON (((110.8753 -2...
## 4 tidak MULTIPOLYGON (((110.875 -2....
## 5 tidak MULTIPOLYGON (((110.8748 -2...
## 6 tidak MULTIPOLYGON (((110.8747 -2...
# 3. Spatial Join: Tambahkan info kk_in ke masing-masing desa
desa_kawasan_join <- st_join(desa_coords, fungsi_kawasan, join = st_intersects) %>%
  st_drop_geometry() %>%
  select(id_desa, kk_in) %>%
  group_by(id_desa) %>%
  summarise(
    kk_in = ifelse(any(kk_in == "ya", na.rm = TRUE), "ya", "tidak"),
    .groups = "drop"
  )

# 4. Gabungkan semua info ke data_gabungan
data_gabungan <- ndvi %>%
  mutate(id_desa = str_to_lower(str_replace_all(id_desa, "\\s+", "_"))) %>%
  left_join(
    mangrove_kk %>%
      mutate(id_desa = str_to_lower(str_replace_all(id_desa, "\\s+", "_"))),
    by = "id_desa"
  ) %>%
  left_join(desa_kawasan_join, by = "id_desa")

# Tambahkan kolom kabupaten dari shapefile ke data_gabungan
desa_kabupaten <- desa_coords %>%
  st_drop_geometry() %>%
  select(id_desa, kabupaten)

data_gabungan <- data_gabungan %>%
  left_join(desa_kabupaten, by = "id_desa")

# 5. Siapkan Desa Intervensi & Kontrol
desa_intervensi <- data_gabungan %>%
  filter(tolower(intervensi) == "ya") %>%
  select(id_desa, kabupaten, tahun_data, mean_ndvi, mangrove, status_kk, kk_in)

# 6. Kabupaten & KK intervensi untuk filter kontrol
kab_intervensi <- unique(desa_intervensi$kabupaten)
kk_intervensi  <- unique(desa_intervensi$kk_in[desa_intervensi$kk_in == "ya"])

# 7. Desa kontrol eligible
desa_kontrol <- data_gabungan %>%
  filter(
    tolower(intervensi) != "ya",
    !kabupaten %in% kab_intervensi,
    !(kk_in %in% kk_intervensi)
  ) %>%
  select(id_desa, kabupaten, tahun_data, mean_ndvi, mangrove, status_kk)

# 8. Fungsi Matching DTW Multivariat
cari_kontrol_terdekat <- function(desa_id, kabupaten_desa, mangrove_val, status_kk_val, ts_desa_intervensi) {
  kandidat_kontrol <- desa_kontrol %>%
    filter(
      mangrove == mangrove_val,
      status_kk == status_kk_val,
      !is.na(mean_ndvi)
    )
  
  ts_kontrol_list <- kandidat_kontrol %>%
    arrange(tahun_data) %>%
    group_by(id_desa) %>%
    summarise(ts_kontrol = list(mean_ndvi), .groups = "drop")
  
  if (nrow(ts_kontrol_list) == 0) return(NULL)
  
  hasil <- ts_kontrol_list %>%
    mutate(nilai_dtw = map_dbl(ts_kontrol, ~ dtw(ts_desa_intervensi, .x)$distance)) %>%
    arrange(nilai_dtw) %>%
    slice_head(n = 3) %>%
    mutate(id_intervensi = desa_id)
  
  return(hasil)
}

# 9. Jalankan Matching
hasil_semua <- list()

for (desa_id in unique(desa_intervensi$id_desa)) {
  desa_data     <- desa_intervensi %>% filter(id_desa == desa_id)
  kabupaten_desa <- unique(desa_data$kabupaten)
  mangrove_val   <- unique(desa_data$mangrove)
  status_kk_val  <- unique(desa_data$status_kk)
  
  ts_desa <- desa_data %>%
    arrange(tahun_data) %>%
    pull(mean_ndvi)
  
  hasil <- cari_kontrol_terdekat(desa_id, kabupaten_desa, mangrove_val, status_kk_val, ts_desa)
  
  if (!is.null(hasil)) hasil_semua[[desa_id]] <- hasil
}

dtw_matching_result <- bind_rows(hasil_semua)

# 10. Tambahkan Metadata Nama Desa & Kabupaten
final_result <- dtw_matching_result %>%
  left_join(
    data_gabungan %>%
      select(id_desa, desa_kontrol = desa, kabupaten_kontrol = kabupaten) %>%
      distinct(id_desa, .keep_all = TRUE),
    by = "id_desa"
  ) %>%
  left_join(
    data_gabungan %>%
      select(id_desa, desa_intervensi = desa, kabupaten_intervensi = kabupaten) %>%
      distinct(id_desa, .keep_all = TRUE),
    by = c("id_intervensi" = "id_desa")
  )

# 11. Tambahkan Mangrove dan Status_KK
info_mangrove_kk <- data_gabungan %>%
  select(id_desa, mangrove, status_kk) %>%
  distinct(id_desa, .keep_all = TRUE)

final_result <- final_result %>%
  left_join(
    info_mangrove_kk %>%
      rename(mangrove_kontrol = mangrove, status_kk_kontrol = status_kk),
    by = "id_desa"
  ) %>%
  left_join(
    info_mangrove_kk %>%
      rename(
        id_intervensi = id_desa,
        mangrove_intervensi = mangrove,
        status_kk_intervensi = status_kk
      ),
    by = "id_intervensi"
  )

# 12. Ranking Hasil Matching
final_result <- final_result %>%
  group_by(id_intervensi) %>%
  arrange(nilai_dtw, .by_group = TRUE) %>%
  mutate(rank_kontrol = row_number()) %>%
  ungroup()

# 13. Berapa banyak desa intervensi unik?
n_desa_intervensi <- final_result %>% distinct(id_intervensi) %>% nrow()

# 14. Berapa banyak total baris hasil matching?
n_matching_total <- nrow(final_result)

cat("Desa intervensi:", n_desa_intervensi,
    "\nTotal baris matching:", n_matching_total,
    "\nSeharusnya:", n_desa_intervensi * 3, "\n")
## Desa intervensi: 27 
## Total baris matching: 81 
## Seharusnya: 81

Hasil Sementara Matching Desa Intervensi dan Kontrol Berdasarkan nilai DTW

Analisis awal terhadap hasil pencocokan desa intervensi dengan desa kontrol menunjukkan bahwa:

  • Proses pencocokan menggunakan metode Dynamic Time Warping (DTW) telah menghasilkan 81 baris data pencocokan, yang merepresentasikan hubungan antara desa intervensi dan desa kontrol potensial berdasarkan kemiripan pola waktu NDVI.
  • Berdasarkan jumlah baris tersebut, diasumsikan terdapat 27 desa intervensi, dengan masing-masing desa dipasangkan dengan 3 desa kontrol. Namun, jumlah desa intervensi yang diproses masih perlu diverifikasi ulang untuk memastikan kesesuaiannya.
  • Untuk setiap desa intervensi, desa kontrol diperingkat berdasarkan nilai DTW, di mana:
    • Ranking 1 menunjukkan desa kontrol dengan pola NDVI yang paling mirip dengan desa intervensi, yaitu dengan nilai DTW terkecil.
    • Ranking 2 dan 3 mewakili desa kontrol berikutnya dengan nilai DTW yang masih relatif kecil, namun tidak semirip ranking 1.

Validasi Data Desa Intervensi Sesudah Proses DTW

Validasi memastikan bahwa seluruh desa intervensi telah berhasil dipasangkan dengan desa kontrol setelah proses pencocokan berbasis Dynamic Time Warping (DTW). Berikut ini cara sederhana untuk melakukan validasi dengan cek total data sebelum dan sesudah proses Dynamic Time Warping (DTW):

# 15. Cek total desa intervensi sebelum proses DTW
desa_intervensi_awal <- data_gabungan %>%
  filter(tolower(intervensi) == "ya") %>%
  distinct(id_desa)

nrow(desa_intervensi_awal)
## [1] 29

Hasil ini menunjukkan terdapat 29 desa intervensi unik yang teridentifikasi (total Desa yang menerima program/intervensi oleh YPI, 29 diluar Desa Togan Baru dan Desa Ganjang karena program baru berjalan di 2025).

# 16. Cek Desa Intervensi yang Tidak Terproses dalam Matching
id_intervensi_awal  <- desa_intervensi_awal %>% pull(id_desa)
id_intervensi_final <- final_result %>% distinct(id_intervensi) %>% pull(id_intervensi)

id_intervensi_tak_terproses <- setdiff(id_intervensi_awal, id_intervensi_final)

cat("Jumlah desa intervensi tidak terproses:", length(id_intervensi_tak_terproses), "\n")
## Jumlah desa intervensi tidak terproses: 2

Pengecekan kembali terhadap ID desa intervensi yang berhasil diproses. Hasilnya menunjukkan bahwa:

  • Jumlah desa intervensi yang tidak terproses dalam pencocokan DTW: 2 desa
  • Dengan demikian, jumlah desa intervensi yang berhasil dipasangkan dengan desa kontrol: 27 desa
# 17. Simpan Hasil
final_result_clean <- final_result %>%
  mutate(across(where(is.list), ~ sapply(., function(x) if (length(x) == 1) x else NA)))

str(final_result_clean)
## tibble [81 × 13] (S3: tbl_df/tbl/data.frame)
##  $ id_desa             : chr [1:81] "singkawang_timur_kelurahan_bagak_sahwa" "singkawang_timur_kelurahan_nyarumkop" "singkawang_timur_kelurahan_sanggau_kulor" "singkawang_timur_kelurahan_sanggau_kulor" ...
##  $ ts_kontrol          : logi [1:81] NA NA NA NA NA NA ...
##  $ nilai_dtw           : num [1:81] 0.346 0.372 0.386 0.333 0.376 ...
##  $ id_intervensi       : chr [1:81] "air_besar_bentiang_madomang" "air_besar_bentiang_madomang" "air_besar_bentiang_madomang" "air_besar_dange_aji" ...
##  $ desa_kontrol        : chr [1:81] "kelurahan_bagak_sahwa" "kelurahan_nyarumkop" "kelurahan_sanggau_kulor" "kelurahan_sanggau_kulor" ...
##  $ kabupaten_kontrol   : chr [1:81] "kota_singkawang" "kota_singkawang" "kota_singkawang" "kota_singkawang" ...
##  $ desa_intervensi     : chr [1:81] "bentiang_madomang" "bentiang_madomang" "bentiang_madomang" "dange_aji" ...
##  $ kabupaten_intervensi: chr [1:81] "landak" "landak" "landak" "landak" ...
##  $ mangrove_kontrol    : chr [1:81] "tidak" "tidak" "tidak" "tidak" ...
##  $ status_kk_kontrol   : chr [1:81] "CA" "CA" "CA" "CA" ...
##  $ mangrove_intervensi : chr [1:81] "tidak" "tidak" "tidak" "tidak" ...
##  $ status_kk_intervensi: chr [1:81] "CA" "CA" "CA" "CA" ...
##  $ rank_kontrol        : int [1:81] 1 2 3 1 2 3 1 2 3 1 ...
write_csv(final_result_clean, "output/dtw_matching_mangrove_kk_ndvi.csv")

https://docs.google.com/spreadsheets/d/1WE_48VQrsaESXMA31qIFeGky6WwqN1Nbp9Qcg35UcJk/edit?gid=2047071437#gid=2047071437

Bagian C: Analisis & Visualisasi Tambahan

Bagian ini berfokus pada analisis lanjutan terhadap hasil pencocokan desa intervensi dan desa kontrol.
Distribusi nilai DTW antara pasangan desa dianalisis untuk melihat sebaran kesamaan pola NDVI yang dihasilkan oleh proses matching.

Beberapa pendekatan visualisasi digunakan:
- Histogram untuk melihat distribusi nilai DTW secara keseluruhan.
- Boxplot untuk memeriksa variasi nilai DTW berdasarkan urutan ranking kontrol (1, 2, dan 3).
- Scatter plot per desa intervensi untuk mengevaluasi variasi nilai antar pasangan.

Selain nilai DTW, dihitung juga nilai Root Mean Square Error (RMSE) pada periode tahun yang sama untuk menilai tingkat kesamaan NDVI antar desa secara kuantitatif.
Hasil analisis divisualisasikan menggunakan jaringan (visNetwork) dan peta interaktif (leaflet) untuk menunjukkan keterkaitan spasial antara desa intervensi dan desa kontrol.

# Distribusi nilai DTW
# Lihat sebaran nilai DTW antara desa intervensi dan kontrol yang dipilih.

# 1. Gabungkan semua nilai ke satu kolom
hasil_dtw_long_2 <- final_result %>%
  pivot_longer(
    cols      = starts_with("nilai"),
    names_to  = "ranking",
    values_to = "nilai_dtw"
  )

# 2. Plot distribusi nilai DTW
ggplot(hasil_dtw_long_2, aes(x = nilai_dtw, fill = factor(ranking))) +
  geom_histogram(bins = 30, position = "identity", alpha = 0.6, color = "white") +
  scale_fill_brewer(palette = "Set1") +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(
    title = "Distribusi Nilai Dynamic Time Warping (DTW)",
    x     = "Nilai DTW",
    y     = "Jumlah Desa Kontrol"
  )

# 3. Statistik Ringkasan nilai DTW
summary(hasil_dtw_long_2$nilai_dtw)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1829  0.2563  0.3637  0.5307  0.7738  1.5486

Distribusi Nilai Dynamic Time Warping (DTW)

Histogram di atas menunjukkan sebaran nilai Dynamic Time Warping (DTW) antara desa intervensi dan desa kontrol yang telah dipilih. Nilai DTW merepresentasikan tingkat kemiripan pola temporal antar desa: semakin rendah nilai DTW, semakin mirip pola perubahannya.

Ringkasan Statistik Nilai DTW

  • Minimum: 0.18
  • Kuartil 1 (Q1): 0.26
  • Median: 0.36
  • Rata-rata: 0.53
  • Kuartil 3 (Q3): 0.77
  • Maksimum: 1.55

Distribusi

  • Histogram menunjukkan bahwa sebagian besar desa kontrol memiliki nilai DTW antara 0.2 hingga 0.5, yang menunjukkan kemiripan pola temporal yang cukup tinggi dengan desa intervensi.
  • Terjadi penurunan frekuensi yang signifikan setelah nilai DTW melewati 0.5, dan semakin jarang setelah 0.8, menunjukkan bahwa hanya sebagian kecil desa yang memiliki pola yang sangat berbeda.
  • Distribusi bersifat right-skewed (condong ke kanan), ini berarti sebagian besar data (frekuensi tinggi) berada di nilai kecil (kiri grafik), sedikit data memiliki nilai besar. (ekor ke kanan)
# 4. Visualisasi Ranking vs nilai DTW
ggplot(final_result, aes(x = rank_kontrol, y = nilai_dtw)) +
  geom_boxplot(aes(group = rank_kontrol, fill = as.factor(rank_kontrol))) +
  theme_minimal() +
  labs(
    title = "Nilai DTW Berdasarkan Ranking",
    x     = "Ranking",
    y     = "Nilai DTW"
  )

Boxplot di atas menunjukkan distribusi nilai Dynamic Time Warping (DTW) untuk masing-masing ranking desa kontrol:

  • Ranking 1 memiliki nilai DTW paling rendah dan sebaran yang sempit.
    Ini menunjukkan bahwa desa-desa dalam kelompok ini memiliki pola temporal yang paling mirip dan konsisten dengan desa intervensi.

  • Ranking 2 menunjukkan sedikit peningkatan nilai DTW dan sebaran yang lebih lebar dibanding ranking 1,
    menandakan variasi kemiripan yang mulai muncul antar desa kontrol terhadap desa intervensi.

  • Ranking 3 memiliki nilai DTW tertinggi dengan sebaran paling lebar,
    mengindikasikan bahwa desa-desa dalam kelompok ini memiliki pola temporal yang berbeda-beda dan cenderung tidak konsisten terhadap desa intervensi.

Semakin tinggi ranking (dari 1 ke 3), nilai DTW cenderung meningkat dan sebarannya semakin lebar.
Hal ini menunjukkan bahwa desa dengan ranking lebih rendah (ranking 3) kurang cocok dijadikan desa kontrol, karena pola indikator lingkungannya tidak cukup mirip dengan desa intervensi.

# 5. Cek Distribusi per Intervensi
final_result %>%
  ggplot(aes(x = id_intervensi, y = nilai_dtw, color = as.factor(rank_kontrol))) +
  geom_point() +
  coord_flip() +
  theme_minimal() +
  labs(
    title = "Nilai DTW per Desa Intervensi",
    x     = "Desa Intervensi",
    y     = "Nilai DTW"
  )

Sebaran nilai Dynamic Time Warping (DTW) pada grafik di atas menunjukkan kedekatan pola temporal antara masing-masing desa intervensi dan desa kontrol, berdasarkan tiga peringkat teratas (Rank 1, 2, dan 3). Titik-titik pada grafik mewakili nilai DTW dari pasangan desa intervensi dan kontrol pada masing-masing peringkat.

  • Rank 1 menunjukkan rentang nilai DTW yang umumnya lebih sempit, berkisar antara 0.3 hingga 1.2. Ini menandakan bahwa desa kontrol pada peringkat pertama cenderung memiliki pola temporal yang paling mirip dan konsisten dengan desa intervensi.

  • Rank 2 memiliki rentang nilai DTW yang lebih lebar, yakni antara 0.4 hingga 1.5. Hal ini mengindikasikan bahwa desa kontrol peringkat kedua mulai menunjukkan variasi pola temporal yang lebih besar terhadap desa intervensi.

  • Rank 3 menunjukkan sebaran nilai DTW yang paling luas, bahkan mencapai nilai di atas 1.5 untuk beberapa kasus. Ini menunjukkan bahwa desa kontrol pada peringkat ketiga memiliki pola temporal yang paling berbeda dari desa intervensi.

Sebaran nilai DTW yang semakin lebar pada peringkat yang lebih rendah mencerminkan bahwa kesesuaian pola temporal antar desa semakin berkurang. Sebaliknya, sebaran yang sempit menunjukkan konsistensi yang lebih kuat dalam kemiripan pola.

Dengan demikian, pemilihan desa kontrol yang valid sebaiknya mempertimbangkan nilai DTW yang rendah serta kestabilan pola temporalnya, agar interpretasi dampak intervensi lebih akurat dan dapat dipertanggungjawabkan.

Desa Intervensi yang tidak memiliki Desa Kontrol

Dalam tahap pencocokan desa intervensi dengan desa kontrol, ditemukan dua desa yang tidak berhasil mendapatkan pasangan kontrol yang sesuai. Hal ini diidentifikasi melalui fungsi anti_join, yang membandingkan daftar awal desa intervensi dengan hasil akhir pencocokan (berdasarkan kemiripan pola temporal indikator lingkungan). Kedua desa tersebut adalah:

# 6. Analisis Kenapa 2 Desa Gagal
desa_gagal <- anti_join(
  desa_intervensi_awal,
  final_result %>% distinct(id_intervensi),
  by = c("id_desa" = "id_intervensi")
)

data_gabungan %>%
  filter(id_desa %in% desa_gagal$id_desa) %>%
  select(id_desa, kabupaten, mangrove, status_kk, kk_in) %>%
  distinct()
## # A tibble: 2 × 5
##   id_desa                   kabupaten    mangrove status_kk kk_in
##   <chr>                     <chr>        <chr>    <chr>     <chr>
## 1 kepulauan_karimata_betok  kayong_utara ya       CA        ya   
## 2 kepulauan_karimata_padang kayong_utara ya       CA        ya
  • Kepulauan Karimata - Betok
  • Kepulauan Karimata - Padang

Keduanya berada di Kabupaten Kayong Utara, memiliki ekosistem mangrove, dan termasuk di dalam kawasan konservasi Cagar Alam (CA). Karakteristik ini menjadikan kedua desa tersebut unik secara spasial dan ekologis, sehingga sulit menemukan desa kontrol lain dengan kondisi yang benar-benar sebanding.

Setelah ditelaah lebih lanjut terhadap data mentah (raw data), ditemukan beberapa poin penting berikut:

  • Cagar Alam Laut Kepulauan Karimata merupakan kawasan konservasi berstatus Cagar Alam (CA) di Kalimantan Barat yang memiliki ekosistem mangrove yang masih relatif utuh.
  • Kawasan CA lainnya di Kalimantan Barat sebagian besar berada di daerah pegunungan, sehingga secara ekologis tidak memungkinkan memiliki mangrove.
  • Ekosistem mangrove di luar kawasan CA banyak ditemukan di kawasan Hutan Lindung (HL) pesisir.

Namun, ditemukan lima desa di Kecamatan Kendawangan yang secara administratif juga berada di dalam kawasan Cagar Alam dan memiliki ekosistem mangrove. Akan tetapi, hasil analisis deret waktu (time-series) NDVI menunjukkan bahwa:

  • Tingkat kehijauan dan tutupan vegetasi pada lima desa di Kendawangan ini jauh lebih rendah dibandingkan dengan desa-desa di Kepulauan Karimata.
  • Pola NDVI menunjukkan karakteristik yang lebih terbuka, yang mengindikasikan kemungkinan adanya tekanan manusia lebih tinggi, degradasi ekosistem, atau konversi lahan.

Oleh karena itu, Meskipun secara spasial dan administratif kelima desa Kendawangan tampak serupa, pola ekologis dan kondisi vegetasi aktualnya tidak sesuai untuk dijadikan desa kontrol yang valid secara metodologis dalam analisis perbandingan dengan desa-desa di Kepulauan Karimata.

Menghitung RMSE

Untuk menilai seberapa mirip perubahan vegetasi (NDVI) antara desa intervensi dan desa kontrol, digunakan metode Root Mean Square Error (RMSE).
Perbandingan dilakukan berdasarkan rata-rata NDVI tahunan selama periode tetap: 2001–2015.

# 7. Menghitung RMSE
# Tahun tetap (range 2001–2015)
tahun_fix <- 2001:2015

# Fungsi menghitung RMSE hanya pada tahun yang sinkron
hitung_rmse_sinkron <- function(id_intv, id_ctrl) {
  ts_i <- desa_intervensi %>%
    filter(id_desa == id_intv, tahun_data %in% tahun_fix)
  
  ts_k <- desa_kontrol %>%
    filter(id_desa == id_ctrl, tahun_data %in% tahun_fix)
  
  # Gabung berdasarkan tahun yang tersedia
  ts_gabung <- inner_join(
    ts_i, ts_k,
    by = "tahun_data",
    suffix = c("_i", "_k")
  )
  
  # Minimal 3 tahun overlap agar sah untuk RMSE
  if (nrow(ts_gabung) >= 3) {
    rmse(ts_gabung$mean_ndvi_i, ts_gabung$mean_ndvi_k)
  } else {
    NA_real_
  }
}

# Tambahkan hasil RMSE ke final_result
final_result <- final_result %>%
  mutate(rmse_ndvi = map2_dbl(id_intervensi, id_desa, hitung_rmse_sinkron))

# Statistik ringkasan RMSE
summary(final_result$rmse_ndvi)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01425 0.02355 0.02977 0.04008 0.05885 0.08967

Distribusi nilai RMSE menunjukkan bahwa:

  • Sekitar 50% pasangan desa memiliki nilai RMSE di bawah 0.030, menunjukkan kemiripan pola NDVI yang cukup tinggi.
  • 25% pasangan desa memiliki nilai RMSE lebih dari 0.059, menandakan ketidaksesuaian yang lebih besar.
  • Nilai rata-rata (0.040) cukup moderat dan tetap berada dalam kisaran yang menunjukkan adanya kesamaan pola pada sebagian besar pasangan.
# 8. Plot Distribusi RMSE
library(ggplot2)

ggplot(final_result, aes(x = rmse_ndvi)) +
  geom_histogram(bins = 20, fill = "steelblue", color = "white", alpha = 0.7) +
  labs(
    title = "Distribusi RMSE NDVI Antar Pasangan (Sinkron Tahun)",
    x     = "RMSE",
    y     = "Jumlah Pasangan"
  ) +
  theme_minimal()

Visualisasi Matching Hasil dari Model DTW (visNetwork)

Bagian ini menyajikan hasil pencocokan desa intervensi dan desa kontrol dalam bentuk jaringan interaktif menggunakan paket visNetwork.
Setiap desa intervensi digambarkan sebagai node (simpul) berwarna tomato (merah), sedangkan desa kontrol dibedakan berdasarkan ranking pencocokan:
- Kontrol Rank 1 (biru tua),
- Kontrol Rank 2 (emas),
- Kontrol Rank 3 (hijau).

Visualisasi ini memudahkan analisis interaktif, seperti:
- Menyoroti hubungan desa tertentu,
- Memfilter kelompok berdasarkan tipe (intervensi vs kontrol),
- Melihat struktur hubungan secara menyeluruh.
Hasilnya memberikan gambaran yang jelas mengenai keterkaitan desa intervensi dengan kandidat desa kontrol yang paling sesuai.

# 1. Siapkan nodes dan edges
# Buat nodes
nodes <- final_result %>%
  transmute(id = id_intervensi, label = desa_intervensi, group = "Intervensi") %>%
  bind_rows(
    final_result %>%
      mutate(group = case_when(
        rank_kontrol == 1 ~ "Kontrol_rank_1",
        rank_kontrol == 2 ~ "Kontrol_rank_2",
        rank_kontrol == 3 ~ "Kontrol_rank_3"
      )) %>%
      transmute(id = id_desa, label = desa_kontrol, group)
  ) %>%
  distinct(id, .keep_all = TRUE)

# 2. Buat edges seperti biasa
edges <- final_result %>%
  transmute(from = id_intervensi,
            to = id_desa,
            width = 2 - scales::rescale(nilai_dtw, to = c(0.2, 1.5)))

# 3. Buat jaringan dengan visNetwork
visNetwork(nodes, edges, height = "650px", width = "100%") %>%
  visNodes(shape = "dot", size = 10) %>%
  visEdges(arrows = "to") %>%
  visGroups(groupname = "Intervensi", color = "tomato") %>%
  visGroups(groupname = "Kontrol_rank_1", color = "steelblue") %>%
  visGroups(groupname = "Kontrol_rank_2", color = "gold") %>%
  visGroups(groupname = "Kontrol_rank_3", color = "forestgreen") %>%
  visOptions(highlightNearest = TRUE, nodesIdSelection = TRUE) %>%
  visLegend(
    addNodes = list(
      list(label = "Intervensi", shape = "dot", color = "tomato"),
      list(label = "Rank 1", shape = "dot", color = "steelblue"),
      list(label = "Rank 2", shape = "dot", color = "gold"),
      list(label = "Rank 3", shape = "dot", color = "forestgreen")
    ),
    useGroups = FALSE,
    position = "right"
  ) %>%
  visLayout(randomSeed = 123)

Visualisasi Matching Hasil dari Model DTW (Leaflet)

Bagian ini memanfaatkan paket leaflet untuk menampilkan hasil pencocokan desa intervensi dan desa kontrol secara spasial pada peta interaktif.
Visualisasi meliputi:

  • Batas desa dalam kawasan konservasi (garis poligon tipis).
  • Kawasan konservasi dengan pewarnaan berbeda berdasarkan fungsinya.
  • Titik desa intervensi yang ditandai dengan warna merah.
  • Titik desa kontrol yang dibedakan berdasarkan ranking pencocokan:
    • Kontrol Rank 1 (biru tua),
    • Kontrol Rank 2 (biru muda),
    • Kontrol Rank 3 (biru sangat muda).
  • Garis penghubung antara desa intervensi dan kontrol yang menunjukkan hubungan hasil pencocokan.

Visualisasi ini memudahkan analisis sebaran spasial hubungan antara desa intervensi dan desa kontrol yang dipilih, sehingga mempermudah pemahaman konteks geografis pencocokan serta mendukung validasi visual hasil model Dynamic Time Warping (DTW).

# 1 Baca dan Siapkan Data
desa_coords <- st_read("input/desa_dalam_kawasan_konservasi.shp") %>%
  st_transform(4326) %>%
  mutate(
    centroid = st_centroid(geometry),
    lon      = st_coordinates(centroid)[, 1],
    lat      = st_coordinates(centroid)[, 2]
  )
## Reading layer `desa_dalam_kawasan_konservasi' from data source 
##   `D:\MEL\proyek_desa_kontrol\input\desa_dalam_kawasan_konservasi.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 664 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.6773 ymin: -3.067918 xmax: 114.2054 ymax: 2.081301
## Geodetic CRS:  WGS 84
fungsi_kawasan <- st_read("input/kk_kalbar.shp") %>%
  rename(fungsi = Fungsi)
## Reading layer `kk_kalbar' from data source 
##   `D:\MEL\proyek_desa_kontrol\input\kk_kalbar.shp' using driver `ESRI Shapefile'
## Simple feature collection with 861 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 108.6506 ymin: -2.995644 xmax: 114.2059 ymax: 2.082183
## Geodetic CRS:  WGS 84
# 2. Color Palette untuk Fungsi Kawasan
pal_fungsi <- colorFactor("Set2", domain = fungsi_kawasan$fungsi)

# 3. Siapkan Data Koordinat Desa yang Terlibat
desa_terlibat <- final_result %>%
  select(id_intervensi, id_desa, rank_kontrol) %>%
  rename(id_kontrol = id_desa) %>%
  pivot_longer(
    cols      = c(id_intervensi, id_kontrol),
    names_to  = "tipe",
    values_to = "id_desa"
  ) %>%
  distinct(id_desa)

coords_matching <- desa_coords %>%
  filter(id_desa %in% desa_terlibat$id_desa) %>%
  select(id_desa, geometry, lon, lat)

# 4. Titik Desa Intervensi dan Kontrol
intervensi_pts <- coords_matching %>%
  filter(id_desa %in% final_result$id_intervensi)

kontrol_pts <- coords_matching %>%
  inner_join(final_result %>% select(id_desa, rank_kontrol), by = "id_desa")

# 5. Garis Relasi Intervensi ↔ Kontrol
buat_garis <- function(rank) {
  final_result %>%
    filter(rank_kontrol == rank) %>%
    left_join(coords_matching, by = c("id_intervensi" = "id_desa")) %>%
    rename(lon_i = lon, lat_i = lat) %>%
    left_join(coords_matching, by = c("id_desa" = "id_desa")) %>%  # id_desa = id_kontrol
    rename(lon_k = lon, lat_k = lat) %>%
    filter(!is.na(lon_i) & !is.na(lat_i) & !is.na(lon_k) & !is.na(lat_k)) %>%
    mutate(
      line = pmap(
        list(lon_i, lat_i, lon_k, lat_k),
        ~ st_linestring(matrix(c(..1, ..3, ..2, ..4), ncol = 2))
      )
    ) %>%
    mutate(line = st_sfc(line, crs = 4326)) %>%
    st_sf(geometry = .$line)
}

garis_matching_1 <- buat_garis(1)
garis_matching_2 <- buat_garis(2)
garis_matching_3 <- buat_garis(3)
# 6. Visualisasi Leaflet dengan Full Screen
leaflet() %>%
  addTiles() %>%
  
  # Batas Desa
  addPolygons(
    data       = desa_coords,
    fillColor  = "transparent",
    color      = "black",
    weight     = 0.4,
    label      = ~id_desa,
    group      = "Desa Dalam KK"
  ) %>%
  
  # Kawasan Konservasi
  addPolygons(
    data        = fungsi_kawasan,
    fillColor   = ~pal_fungsi(fungsi),
    fillOpacity = 0.4,
    color       = "gray",
    weight      = 1,
    label       = ~fungsi,
    group       = "Status Kawasan"
  ) %>%
  
  # Garis hubungan intervensi ↔ kontrol
  addPolylines(data = garis_matching_1, color = "black",   weight = 2.0, opacity = 0.8, group = "Kontrol Rank 1") %>%
  addPolylines(data = garis_matching_2, color = "orange",  weight = 1.5, opacity = 0.6, group = "Kontrol Rank 2") %>%
  addPolylines(data = garis_matching_3, color = "purple",  weight = 1.5, opacity = 0.5, group = "Kontrol Rank 3") %>%
  
  # Titik Desa Intervensi
  addCircleMarkers(
    data        = intervensi_pts,
    lng         = ~lon,
    lat         = ~lat,
    color       = "red",
    radius      = 6,
    stroke      = TRUE,
    fillOpacity = 1,
    label       = ~id_desa,
    group       = "Desa Intervensi"
  ) %>%
  
  # Titik Desa Kontrol per rank
  {
    for (rank in 1:3) {
      pts   <- kontrol_pts %>% filter(rank_kontrol == rank)
      warna <- c("blue", "dodgerblue", "lightblue")[rank]
      grp   <- paste0("Kontrol Rank ", rank)
      . <- addCircleMarkers(
        .,
        data        = pts,
        lng         = ~lon,
        lat         = ~lat,
        color       = warna,
        radius      = 5,
        stroke      = TRUE,
        fillOpacity = 0.9,
        label       = ~id_desa,
        group       = grp
      )
    }
    .
  } %>%
  
  # Judul di pojok bawah
  addControl(
    html = "
      <div style='background-color: white; padding: 6px; font-size: 14px; border-radius: 6px;'>
        <b>Sebaran Desa Kontrol Terhadap Desa Intervensi (Matching)</b><br>
        Berdasarkan metode Dynamic Time Warping<br>
        (Kesamaan Status KK, Keberadaan ekosistem Mangrove dan NDVI 2001–2015)
      </div>",
    position = "bottomleft"
  ) %>%
  
  # Layer Control
  addLayersControl(
    overlayGroups = c(
      "Desa Dalam KK", "Status Kawasan",
      "Desa Intervensi", "Kontrol Rank 1",
      "Kontrol Rank 2", "Kontrol Rank 3"
    ),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%

  # Tambahkan tombol full screen
  addFullscreenControl(position = "topleft")

Referensi

Lonita, B. I., Prasetyo, Y., & Hani’ah. (2015). Analisis perubahan luas dan kerapatan hutan menggunakan algoritma NDVI dan EVI pada citra Landsat 7 ETM+ tahun 2006, 2009, dan 2012 (Studi kasus: Kabupaten Kendal, Provinsi Jawa Tengah). Jurnal Geodesi UNDIP, 4(3), 112–120. https://doi.org/10.14710/jgundip.2015.8965
(https://ejournal3.undip.ac.id/index.php/geodesi/article/view/8965/8712)

Sofro, A., Riani, R., Khikmah, K., Romadhonia, R., & Ariyanto, D. (2024). Analysis of rainfall in Indonesia using a time series‑based clustering approach. BAREKENG: Journal of Mathematics & Applied Science, 18(2), 837–848. https://ojs3.unpatti.ac.id/index.php/barekeng/article/view/11124/7822 Studi ini menggunakan DTW sebagai metode untuk mengukur kemiripan pola deret waktu antar wilayah geografis (provinsi).

Wikipedia contributors. (2024). Dynamic time warping. Wikipedia. https://en.wikipedia.org/wiki/Dynamic_time_warping

For More Information

Arif Darmawan
MEL Technical Supervisor
Yayasan Planet Indonesia

🌐 www.planetindonesia.org
📧
📧