Percobaan Pertama: Kasus Jurnal Pertama

Latar belakang

Salah satu sumber makanan pokok yang dikonsumsi masyarakat Indonesia adalah berasal dari sektor pertanian yaitu padi. Sebagian besar penduduk Indonesia, terutama di daerah pedesaan, telah lama bergantung pada pertanian sebagai sumber pendapatan utama. Kabupaten Majalengka sebagai salah satu wilayah agraris di Jawa Barat memiliki potensi besar dalam sektor pertanian. Sejalan dengan ini, Open Data Majalengka hadir sebagai platform terbuka yang menyajikan informasi akurat dari berbagai organisasi di Kabupaten Majalengka. Permasalahan yang dihadapi terkait sulitnya memahami data pertanian yang banyak dan bervariasi. Hal ini menciptakan hambatan dalam analisis dan pemahaman mendalam terhadap situasi pertanian di Kabupaten Majalengka. Data produksi yang terpisah- pisah berdasarkan jenis komoditas menyulitkan proses penggalian informasi yang efisien dan efektif. Oleh karena itu diusulkan penelitian dengan judul “Visualisasi Data Produksi Pertanian Melalui Penerapan Media Interaktif Menggunakan Tableau“

Persiapan Data Penelitian

Visualisasi pertama

Persiapan library yang digunakan terlebih dahulu

# Library yang digunakan dalam kasus ketersediaan pangan 
library(readxl)
library(dplyr)
## 
## 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(tidyr)
library(ggplot2)
library(patchwork)
library(sf)
## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
library(geodata)
## Loading required package: terra
## terra 1.9.27
## 
## Attaching package: 'terra'
## The following object is masked from 'package:patchwork':
## 
##     area
## The following object is masked from 'package:tidyr':
## 
##     extract
library(stringr)
library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:terra':
## 
##     rescale

Langkah berikutnya membaca data excel terlebih dahulu

## membaca data excel yang didapat
data <- readxl::read_excel(
  "/Users/M.Fabian.R.D/Desktop/SEMESTER 7/KETEHANAN PANGAN BERBASIS DATA]/KETERSEDIAAN PANGAN/DATA EXCEL/Data_tugas_ketersediaan_pangan.xlsx"
)

# Cek struktur data
names(data)
str(data)

# Filter data Kabupaten Majalengka
data_majalengka <- data %>%
  dplyr::filter(
    bps_nama_kabupaten_kota == "KABUPATEN MAJALENGKA",
    !is.na(bps_nama_kecamatan)
  ) %>%
  dplyr::mutate(
    bps_nama_kecamatan = stringr::str_to_title(
      bps_nama_kecamatan
    )
  )

# Cek jumlah kecamatan
dplyr::n_distinct(data_majalengka$bps_nama_kecamatan)

# Cek nama kecamatan
sort(unique(data_majalengka$bps_nama_kecamatan))

# Cek ringkasan produksi pertanian
summary(data_majalengka$jumlah_produksi_pertanian)

Menyiapkan data komoditas majalengka

## menyiapkan data komoditas
kolom_komoditas <- c(
  "jumlah_produksi_padi",
  "jumlah_produksi_kopi_arabika",
  "jumlah_produksi_kopi_robusta",
  "jumlah_produksi_kacang_kedelai",
  "jumlah_produksi_cabai_besar",
  "jumlah_produksi_ubi_jalar",
  "jumlah_produksi_ubi_kayu",
  "jumlah_produksi_tomat",
  "jumlah_produksi_bawang_merah",
  "jumlah_produksi_jagung",
  "jumlah_produksi_cabai_rawit",
  "jumlah_produksi_pertanian"
)

# Pastikan kolom produksi berbentuk numerik
data_majalengka <- data_majalengka %>%
  mutate(
    across(
      all_of(kolom_komoditas),
      as.numeric
    )
  )

Visualisasi peta produksi bahan pangan di majalengka

# VISUALISASI PETA PRODUKSI PADI
indo <- geodata::gadm(
  country = "IDN",
  level = 3,
  path = tempdir()
)
## Cached as: /var/folders/tg/mlw6kz8515df8t8pytbhzf3r0000gp/T//RtmpY1ylWM/gadm/gadm41_IDN_3_pk.rds
indo_sf <- st_as_sf(indo)

## Ambil wilayah Kabupaten Majalengka
majalengka_map <- indo_sf %>%
  filter(NAME_2 == "Majalengka") %>%
  mutate(
    kecamatan = str_to_title(
      str_squish(NAME_3)
    )
  )

## Cek nama kecamatan yang belum cocok
nama_excel <- data_majalengka %>%
  distinct(bps_nama_kecamatan)

nama_peta <- majalengka_map %>%
  st_drop_geometry() %>%
  distinct(kecamatan)

print(
  anti_join(
    nama_excel,
    nama_peta,
    by = c(
      "bps_nama_kecamatan" = "kecamatan"
    )
  )
)

Menghitung nilai produksi pertanian di majalengka

## Hitung produksi padi per kecamatan
produksi_pertanian <- data_majalengka %>%
  group_by(bps_nama_kecamatan) %>%
  summarise(
    total_produksi_pertanian = sum(
      jumlah_produksi_pertanian,
      na.rm = TRUE
    ),
    .groups = "drop"
  )

# Menggabungkan data produksi dengan peta
peta_data <- majalengka_map %>%
  left_join(
    produksi_pertanian,
    by = c("kecamatan" = "bps_nama_kecamatan")
  )

# Membuat peta persebaran produksi pertanian
plot_peta <- ggplot(peta_data) +
  geom_sf(
    aes(fill = total_produksi_pertanian),
    color = "white",
    linewidth = 0.3
  ) +
  scale_fill_gradient(
    low = "#E5F5E0",
    high = "#238443",
    na.value = "grey85",
    labels = scales::label_number(big.mark = "."),
    name = "Produksi pertanian"
  ) +
  labs(
    title = "Persebaran Produksi Pertanian",
    subtitle = "Kabupaten Majalengka"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(
      face = "bold",
      size = 12
    ),
    axis.title = element_blank(),
    axis.text = element_blank(),
    panel.grid = element_blank(),
    legend.position = "right"
  )

plot_peta
## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both '.', which could be confusing

terlihat visualisasi peta kecamatan di Kabupaten Majalengka dengan jumlah produksi pertanian. Kecamatan berwarna hijau tebal menunjukkan tingginya produksi pertanian di wilayah tersebut. Jumlah produksi pertanian tersebut dihasilkan dari jumlah total jenis produksi yang bervariasi yaitu padi, jagung, bawang merah, cabai rawit, cabai besar, ubi jalar, ubi kayu, kenang, tomat, kacang kedelai, kopi robusta, dan kopi arabika.

Visualisasi produksi pertanian per kecamatan di majalengka

## VISUALISASI PRODUKSI PERTANIAN PER KECAMATAN
plot_pertanian <- produksi_pertanian %>%
  ggplot(
    aes(
      x = reorder(
        bps_nama_kecamatan,
        total_produksi_pertanian
      ),
      y = total_produksi_pertanian
    )
  ) +
  geom_col(
    fill = "#79BD67",
    width = 0.75
  ) +
  coord_flip() +
  scale_y_continuous(
    labels = scales::label_number(big.mark = ".")
  ) +
  labs(
    title = "Produksi Pertanian per Kecamatan",
    subtitle = "Kabupaten Majalengka",
    x = NULL,
    y = "Total Produksi Pertanian"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(
      face = "bold",
      size = 12
    ),
    plot.subtitle = element_text(size = 10),
    panel.grid.minor = element_blank()
  )

print(plot_pertanian)
## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both '.', which could be confusing

Kecamatan Kertajati merupakan daerah dengan produksi pertanian tertinggi, dari total jenis produksi pertanian pada tahun 2022. Data ini menunjukkan kontribusi signifikan Kecamatan Kertajati dalam hasil pertanian di Kabupaten Majalengka.

Visual per komoditas produksi

## VISUALISASI TOTAL PRODUKSI PER KOMODITAS
produksi_komoditas <- data_majalengka %>%
  summarise(
    across(
      all_of(kolom_komoditas),
      ~ sum(.x, na.rm = TRUE)
    )
  ) %>%
  pivot_longer(
    cols = everything(),
    names_to = "komoditas",
    values_to = "produksi"
  ) %>%
  mutate(
    komoditas = str_remove(
      komoditas,
      "^jumlah_produksi_"
    ),
    komoditas = str_replace_all(
      komoditas,
      "_",
      " "
    ),
    komoditas = str_to_title(komoditas)
  ) %>%
  arrange(produksi)

plot_komoditas <- ggplot(
  produksi_komoditas,
  aes(
    x = reorder(komoditas, produksi),
    y = produksi
  )
) +
  geom_col(
    fill = "#80B96A",
    width = 0.75
  ) +
  coord_flip() +
  scale_y_continuous(
    labels = label_number(big.mark = ".")
  ) +
  labs(
    title = "Total Produksi per Komoditas",
    x = NULL,
    y = "Jumlah produksi"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(
      face = "bold",
      size = 12
    ),
    panel.grid.minor = element_blank()
  )

print(produksi_komoditas)
## # A tibble: 12 × 2
##    komoditas      produksi
##    <chr>             <dbl>
##  1 Kopi Arabika      1231.
##  2 Kopi Robusta      4740.
##  3 Kacang Kedelai   13248 
##  4 Tomat            16597.
##  5 Ubi Kayu         23179 
##  6 Cabai Besar      32803.
##  7 Cabai Rawit      38437.
##  8 Ubi Jalar        51890 
##  9 Bawang Merah    173693.
## 10 Jagung          679038 
## 11 Padi           3592015 
## 12 Pertanian      4626872.
print(plot_komoditas)
## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both '.', which could be confusing

produksi padi menjadi komoditas dengan produksi tertinggi pada tahun 2022. Data ini mencerminkan kontribusi signifikan jenis produksi padi terhadap hasil total pada tahun 2022.

Membuat visualisasi gabungan

## MEMBUAT DASHBOARD GABUNGAN
dashboard <- (
  plot_peta | plot_pertanian
) /
  plot_komoditas +
  plot_layout(
    heights = c(1.2, 0.8)
  ) +
  plot_annotation(
    title = paste(
      "DASHBOARD PRODUKSI PERTANIAN",
      "KABUPATEN MAJALENGKA"
    ),
    subtitle = paste(
      "Visualisasi ketersediaan pangan",
      "berdasarkan produksi komoditas pertanian"
    ),
    caption = "Sumber: Data tugas ketersediaan pangan"
  ) &
  theme(
    plot.title = element_text(
      face = "bold",
      size = 16
    ),
    plot.subtitle = element_text(
      size = 11
    ),
    plot.caption = element_text(
      size = 9
    )
  )

# Menampilkan dashboard
print(dashboard)
## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both '.', which could be confusing
## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both '.', which could be confusing
## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both '.', which could be confusing

# Menyimpan dashboard sebagai gambar
ggsave(
  filename = "Dashboard_Produksi_Majalengka.png",
  plot = dashboard,
  width = 16,
  height = 11,
  units = "in",
  dpi = 300,
  bg = "white"
)
## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both '.', which could be confusing
## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both '.', which could be confusing
## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' and 'decimal.mark' are both '.', which could be confusing

Berdasarkan pembahasan sebelumnya, dapat disimpulkan bahwa transformasi data yang menggabungkan data terkait menjadi satu kesatuan memungkinkan visualisasi data yang informatif melalui grafik yang menarik. Visualisasi data produksi pertanian, seperti peta persebaran produksi, grafik jumlah produksi per kecamatan, dan grafik jenis komoditas, memberikan gambaran yang jelas.

Percobaan Kedua: Kasus Jurnal Lainnya

latar belakang

Proporsi pengeluaran per kapita merupakan indikator penting kesejahteraan rumah tangga karena mencerminkan kapasitas ekonomi dan pola konsumsi masyarakat, sebagaimana dijelaskan dalam Hukum Engel (Akita & Miyata, 2017). Di wilayah dengan karakteristik geografis yang beragam seperti Provinsi Papua, analisis spasial diperlukan untuk memahami variasi pengeluaran antar kabupaten/kota serta perbedaan antara wilayah perkotaan dan perdesaan. Studi ini bertujuan menganalisis distribusi spasial persentase pengeluaran per kapita untuk barang makanan dan non-makanan di sembilan kabupaten/kota di Provinsi Papua selama periode 2022– 2024.

Visualisasi Kedua

Persiapan library yang digunakan terlebih dahulu

# library yang digunakan
library(readxl)
library(dplyr)
library(ggplot2)
library(tidyr)
library(sf)
library(scales)
library(stringr)
library(readr)
## 
## Attaching package: 'readr'
## The following object is masked from 'package:scales':
## 
##     col_factor
## pengecekam data
## Import data Excel
data <- read_excel("/Users/M.Fabian.R.D/Desktop/SEMESTER 7/KETEHANAN PANGAN BERBASIS DATA]/KETERSEDIAAN PANGAN/DATA EXCEL/gabungan_data_persetase_perkapita_papua.xlsx")

## Rapikan nama kolom
names(data) <- c(
  "Kabupaten_Kota",
  "Tahun",
  "Makanan",
  "Bukan_Makanan"
)

## Pastikan tipe datanya benar
data <- data %>%
  mutate(
    Kabupaten_Kota = trimws(Kabupaten_Kota),
    Tahun = as.integer(Tahun),
    Makanan = as.numeric(Makanan),
    Bukan_Makanan = as.numeric(Bukan_Makanan)
  )
# Kabupaten_Kota, Tahun, Makanan, Bukan_Makanan
data_bersih <- data %>%
  mutate(
    Tahun = as.integer(as.character(Tahun)),
    Makanan = parse_number(
      as.character(Makanan),
      locale = locale(
        decimal_mark = ",",
        grouping_mark = "."
      )
    ),
    Bukan_Makanan = parse_number(
      as.character(Bukan_Makanan),
      locale = locale(
        decimal_mark = ",",
        grouping_mark = "."
      )
    )
  ) %>%
  filter(Tahun %in% c(2022, 2023, 2024))

# 2. Membuat data per kabupaten/kota dan tahun
data_wide <- data_bersih %>%
  select(Kabupaten_Kota, Tahun, Makanan, Bukan_Makanan) %>%
  pivot_longer(
    cols = c(Makanan, Bukan_Makanan),
    names_to = "Jenis_Pengeluaran",
    values_to = "Persentase"
  ) %>%
  mutate(
    Variabel = paste0(
      Tahun, "_",
      ifelse(
        Jenis_Pengeluaran == "Makanan",
        "MKN",
        "NON"
      )
    )
  ) %>%
  select(Kabupaten_Kota, Variabel, Persentase) %>%
  pivot_wider(
    names_from = Variabel,
    values_from = Persentase
  )

# 3. Menentukan urutan variabel
kolom_korelasi <- c(
  "2022_MKN", "2023_MKN", "2024_MKN",
  "2022_NON", "2023_NON", "2024_NON"
)

# 4. Memastikan seluruh kolom numerik
data_numerik <- data_wide %>%
  select(all_of(kolom_korelasi)) %>%
  mutate(
    across(
      everything(),
      as.numeric
    )
  )

# 5. Menghitung korelasi Pearson
matriks_korelasi <- cor(
  data_numerik,
  use = "pairwise.complete.obs",
  method = "pearson"
)

# 6. Mengubah matriks korelasi ke format long
data_heatmap <- as.data.frame(
  as.table(matriks_korelasi)
) %>%
  rename(
    Tahun_X = Var1,
    Tahun_Y = Var2,
    Korelasi = Freq
  ) %>%
  mutate(
    Tahun_X = factor(
      Tahun_X,
      levels = kolom_korelasi
    ),
    Tahun_Y = factor(
      Tahun_Y,
      levels = rev(kolom_korelasi)
    )
  )

# 7. Membuat heatmap
plot_heatmap <- ggplot(
  data_heatmap,
  aes(
    x = Tahun_X,
    y = Tahun_Y,
    fill = Korelasi
  )
) +
  geom_tile(
    color = "white",
    linewidth = 0.4
  ) +
  scale_fill_gradientn(
    colours = c(
      "#440154",
      "#3B528B",
      "#21918C",
      "#5EC962",
      "#FDE725"
    ),
    limits = c(-1, 1),
    name = "Korelasi"
  ) +
  labs(
    title = "Stabilitas Temporal Pola Pengeluaran",
    subtitle = paste(
      "Korelasi pengeluaran makanan dan bukan makanan",
      "di Provinsi Papua, 2022–2024"
    ),
    x = NULL,
    y = NULL,
    caption = "Sumber: BPS, data Susenas"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(
      face = "bold",
      size = 13
    ),
    plot.subtitle = element_text(size = 10),
    axis.text.x = element_text(
      angle = 35,
      hjust = 1,
      size = 9
    ),
    axis.text.y = element_text(size = 9),
    panel.grid = element_blank(),
    legend.position = "right"
  )

# 8. Menampilkan heatmap
print(plot_heatmap)

asil uji korelasi Pearson antara pengeluaran makanan dan non-makanan antar tahun menunjukkan nilai r yang tinggi, berkisar antara 0,85 hingga 0,93. Nilai korelasi ini menunjukkan dua hal penting: Polanya relatif stabil antar tahun, artinya distribusi pengeluaran antar kabupaten/kota tidak mengalami perubahan signifikan meskipun terjadi dinamika makroekonomi. Struktur ekonomi rumah tangga di Papua cenderung kaku, di mana kapasitas konsumsi masyarakat dipengaruhi oleh faktor struktural seperti geografi, aksesibilitas, dan pengembangan infrastruktur. Heatmap korelasi mengindikasikan bahwa pola pengeluaran tahun 2023 memiliki kemiripan tertinggi dengan tahun 2024, sebagaimana tercermin dari nilai korelasi yang paling besar di antara seluruh pasangan tahun, sementara perubahan terbesar terjadi antara 2022 dan 2024 (BPS Indonesia, 2024). Namun, semua nilai korelasi tetap berada dalam kategori sangat kuat, sehingga pola ini dapat dianggap persisten.

# choropleth visualisasi
# Baca shapefile
# Baca shapefile batas wilayah provinsi
# Mengatur tampilan nilai NA
options(na.print = "NA")

# Ambil shapefile kabupaten/kota Indonesia
shp_kab <- geodata::gadm(
  country = "IDN",
  level = 2,
  path = tempdir()
) |>
  sf::st_as_sf()
## Cached as: /var/folders/tg/mlw6kz8515df8t8pytbhzf3r0000gp/T//RtmpY1ylWM/gadm/gadm41_IDN_2_pk.rds
# Cek kolom dan nama provinsi
names(shp_kab)
##  [1] "GID_2"     "GID_0"     "COUNTRY"   "GID_1"     "NAME_1"    "NL_NAME_1"
##  [7] "NAME_2"    "VARNAME_2" "NL_NAME_2" "TYPE_2"    "ENGTYPE_2" "CC_2"     
## [13] "HASC_2"    "geometry"
unique(shp_kab$NAME_1)
##  [1] "Aceh"                "Bali"                "Bangka Belitung"    
##  [4] "Banten"              "Bengkulu"            "Gorontalo"          
##  [7] "Jakarta Raya"        "Jambi"               "Jawa Barat"         
## [10] "Jawa Tengah"         "Jawa Timur"          "Kalimantan Barat"   
## [13] "Kalimantan Selatan"  "Kalimantan Tengah"   "Kalimantan Timur"   
## [16] "Kalimantan Utara"    "Kepulauan Riau"      "Lampung"            
## [19] "Maluku"              "Maluku Utara"        "Nusa Tenggara Barat"
## [22] "Nusa Tenggara Timur" "Papua"               "Papua Barat"        
## [25] "Riau"                "Sulawesi Barat"      "Sulawesi Selatan"   
## [28] "Sulawesi Tengah"     "Sulawesi Tenggara"   "Sulawesi Utara"     
## [31] "Sumatera Barat"      "Sumatera Selatan"    "Sumatera Utara"     
## [34] "Yogyakarta"
# Ambil kabupaten/kota di Provinsi Papua
shp_papua <- shp_kab |>
  dplyr::filter(NAME_1 == "Papua") |>
  dplyr::mutate(
    Kabupaten_Kota = stringr::str_to_lower(
      stringr::str_squish(as.character(NAME_2))
    )
  )

# Bersihkan data Excel
data_peta <- data |>
  dplyr::mutate(
    Kabupaten_Kota = stringr::str_to_lower(
      stringr::str_squish(as.character(Kabupaten_Kota))
    ),
    Tahun = as.integer(as.character(Tahun))
  ) |>
  dplyr::filter(Tahun %in% c(2022, 2023, 2024))

# Cek nama wilayah pada data Excel yang tidak cocok dengan shapefile
wilayah_tidak_cocok <- dplyr::anti_join(
  data_peta |> dplyr::distinct(Kabupaten_Kota),
  sf::st_drop_geometry(shp_papua) |>
    dplyr::distinct(Kabupaten_Kota),
  by = "Kabupaten_Kota"
)

knitr::kable(
  wilayah_tidak_cocok,
  caption = "Wilayah pada Data Excel yang Tidak Cocok dengan Shapefile"
)
Wilayah pada Data Excel yang Tidak Cocok dengan Shapefile
Kabupaten_Kota
# Gabungkan data dengan shapefile
peta_semua <- dplyr::left_join(
  shp_papua,
  data_peta,
  by = "Kabupaten_Kota"
)

# Cek hasil penggabungan
tabel_peta <- peta_semua |>
  sf::st_drop_geometry() |>
  dplyr::select(
    Kabupaten_Kota,
    Tahun,
    Makanan,
    Bukan_Makanan
  )

knitr::kable(
  head(tabel_peta, 30),
  caption = "Hasil Penggabungan Data dengan Shapefile"
)
Hasil Penggabungan Data dengan Shapefile
Kabupaten_Kota Tahun Makanan Bukan_Makanan
asmat NA NA NA
biak numfor 2022 50.42 49.58
biak numfor 2023 47.32 52.68
biak numfor 2024 48.62 51.38
boven digoel NA NA NA
deiyai NA NA NA
dogiyai NA NA NA
intan jaya NA NA NA
jayapura 2022 48.58 51.42
jayapura 2023 55.30 44.70
jayapura 2024 52.93 47.07
jayawijaya NA NA NA
keerom 2022 55.42 44.58
keerom 2023 54.13 45.87
keerom 2024 54.79 45.21
kepulauan yapen 2022 55.42 44.58
kepulauan yapen 2023 54.60 45.40
kepulauan yapen 2024 58.16 41.84
kota jayapura 2022 43.84 56.16
kota jayapura 2023 42.16 57.84
kota jayapura 2024 41.05 58.95
lanny jaya NA NA NA
mamberamo raya 2022 68.63 31.37
mamberamo raya 2023 68.64 31.36
mamberamo raya 2024 68.96 31.04
mamberamo tengah NA NA NA
mappi NA NA NA
merauke NA NA NA
mimika NA NA NA
nabire NA NA NA
# Nama wilayah dari shapefile
nama_peta <- shp_papua |>
  sf::st_drop_geometry() |>
  dplyr::distinct(Kabupaten_Kota)

# Nama wilayah dari Excel
nama_excel <- data_peta |>
  dplyr::distinct(Kabupaten_Kota)

# Wilayah yang hanya terdapat pada shapefile
wilayah_shapefile_saja <- dplyr::anti_join(
  nama_peta,
  nama_excel,
  by = "Kabupaten_Kota"
)

knitr::kable(
  wilayah_shapefile_saja,
  caption = "Wilayah Shapefile yang Tidak Ada di Excel"
)
Wilayah Shapefile yang Tidak Ada di Excel
Kabupaten_Kota
asmat
boven digoel
deiyai
dogiyai
intan jaya
jayawijaya
lanny jaya
mamberamo tengah
mappi
merauke
mimika
nabire
nduga
paniai
pegunungan bintang
puncak
puncak jaya
tolikara
yahukimo
yalimo
# Wilayah yang hanya terdapat pada Excel
wilayah_excel_saja <- dplyr::anti_join(
  nama_excel,
  nama_peta,
  by = "Kabupaten_Kota"
)

knitr::kable(
  wilayah_excel_saja,
  caption = "Wilayah Excel yang Tidak Ada di Shapefile"
)
Wilayah Excel yang Tidak Ada di Shapefile
Kabupaten_Kota
# choropleth makanan ketiga tahun 
# Pastikan hanya tahun 2022, 2023, dan 2024
peta_plot <- peta_semua %>%
  filter(!is.na(Tahun)) %>%
  filter(Tahun %in% c(2022, 2023, 2024)) %>%
  mutate(Tahun = factor(Tahun, levels = c(2022, 2023, 2024)))

# Buat peta per tahun dengan legend masing-masing
for (thn in c(2022, 2023, 2024)) {

  p <- peta_plot %>%
    filter(Tahun == thn) %>%
    ggplot() +
    geom_sf(
      aes(fill = Makanan),
      color = "white",
      linewidth = 0.4
    ) +
    scale_fill_gradient(
      low = "#FFFF00",
      high = "#FF0000",
      name = "Persentase (%)",
      na.value = "grey85"
    ) +
    labs(
      title = "Peta Persentase Pengeluaran Makanan",
      subtitle = paste("Provinsi Papua - Tahun", thn),
      caption = "Sumber: BPS, Susenas"
    ) +
    coord_sf(datum = NA) +
    theme_void() +
    theme(
      plot.title = element_text(
        face = "bold", size = 16
      ),
      plot.subtitle = element_text(size = 13),
      plot.caption = element_text(size = 10),
      legend.position = "right",
      legend.title = element_text(size = 11),
      legend.text = element_text(size = 10)
    )

  print(p)
}

Pengeluaran makanan tertinggi terdapat di Mamberamo Raya, Waropen, dan Kepulauan Yapen, sedangkan Kota Jayapura memiliki proporsi terendah. Pola ini menunjukkan perbedaan konsumsi antara wilayah perkotaan dan pedalaman yang dipengaruhi oleh tingkat pendapatan, aksesibilitas, dan biaya distribusi pangan

# choropleth bukan makanan ketiga tahun 

for (thn in c(2022, 2023, 2024)) {

  p <- peta_plot %>%
    filter(Tahun == thn) %>%
    ggplot() +
    geom_sf(
      aes(fill = Bukan_Makanan),
      color = "white",
      linewidth = 0.4
    ) +
    scale_fill_gradient(
      low = "#FFFF00",
      high = "#FF0000",
      name = "Persentase (%)",
      na.value = "grey85"
    ) +
    labs(
      title = "Peta Persentase Pengeluaran Bukan Makanan",
      subtitle = paste("Provinsi Papua - Tahun", thn),
      caption = "Sumber: BPS, Susenas"
    ) +
    coord_sf(datum = NA) +
    theme_void() +
    theme(
      plot.title = element_text(
        face = "bold", size = 16
      ),
      plot.subtitle = element_text(size = 13),
      plot.caption = element_text(size = 10),
      legend.position = "right",
      legend.title = element_text(size = 11),
      legend.text = element_text(size = 10)
    )

  print(p)
}

Pengeluaran non-pangan tertinggi terdapat di Jayapura dan Keerom, sedangkan Mamberamo Raya, Waropen, dan sebagian wilayah pesisir utara berada pada kelas terendah. Pola ini menunjukkan adanya kesenjangan konsumsi antara wilayah perkotaan dan pedesaan yang dipengaruhi oleh aksesibilitas dan perkembangan ekonomi lokal