Visualisasi Data Spasial Persentase Kemiskinan DIY 2024

Author

Freditasari Purwa Hidayat

Panggil library

#panggil library
library(sf)
Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   4.0.0     ✔ tibble    3.3.0
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.1.0     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggspatial)
library(prettymapr)
library(dplyr)
#Data menggunakan persentase kemiskinan DIY tahun 2024 BPS
#set directory
setwd("D:/FREDITA MAGISTER STAT/Sains Data/Praktikum Sains Data/data-spasial-sains-data")

Import data spasial yaitu kabupaten/kota di Seluruh Indonesia dengan format shp

#Import data spasial
kab_indo <- st_read("Admin2Kabupaten/idn_admbnda_adm2_bps_20200401.shp")
Reading layer `idn_admbnda_adm2_bps_20200401' from data source 
  `D:\FREDITA MAGISTER STAT\Sains Data\Praktikum Sains Data\data-spasial-sains-data\Admin2Kabupaten\idn_admbnda_adm2_bps_20200401.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 522 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
Geodetic CRS:  WGS 84
kab_indo
Simple feature collection with 522 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 95.01079 ymin: -11.00762 xmax: 141.0194 ymax: 6.07693
Geodetic CRS:  WGS 84
First 10 features:
   Shape_Leng Shape_Area         ADM2_EN ADM2_PCODE ADM2_REF ADM2ALT1EN
1    2.360029  0.2289681      Aceh Barat     ID1107     <NA>       <NA>
2    1.963994  0.1541359 Aceh Barat Daya     ID1112     <NA>       <NA>
3    4.590182  0.2363958      Aceh Besar     ID1108     <NA>       <NA>
4    3.287754  0.3161611       Aceh Jaya     ID1116     <NA>       <NA>
5    4.448584  0.3430383    Aceh Selatan     ID1103     <NA>       <NA>
6    4.907219  0.1534414    Aceh Singkil     ID1102     <NA>       <NA>
7    2.593385  0.1745672    Aceh Tamiang     ID1114     <NA>       <NA>
8    3.676889  0.3834894     Aceh Tengah     ID1106     <NA>       <NA>
9    3.473021  0.3374562   Aceh Tenggara     ID1104     <NA>       <NA>
10   5.037148  0.4434042      Aceh Timur     ID1105     <NA>       <NA>
   ADM2ALT2EN ADM1_EN ADM1_PCODE   ADM0_EN ADM0_PCODE       date    validOn
1        <NA>    Aceh       ID11 Indonesia         ID 2019-12-20 2020-04-01
2        <NA>    Aceh       ID11 Indonesia         ID 2019-12-20 2020-04-01
3        <NA>    Aceh       ID11 Indonesia         ID 2019-12-20 2020-04-01
4        <NA>    Aceh       ID11 Indonesia         ID 2019-12-20 2020-04-01
5        <NA>    Aceh       ID11 Indonesia         ID 2019-12-20 2020-04-01
6        <NA>    Aceh       ID11 Indonesia         ID 2019-12-20 2020-04-01
7        <NA>    Aceh       ID11 Indonesia         ID 2019-12-20 2020-04-01
8        <NA>    Aceh       ID11 Indonesia         ID 2019-12-20 2020-04-01
9        <NA>    Aceh       ID11 Indonesia         ID 2019-12-20 2020-04-01
10       <NA>    Aceh       ID11 Indonesia         ID 2019-12-20 2020-04-01
   validTo                       geometry
1     <NA> MULTIPOLYGON (((96.26836 4....
2     <NA> MULTIPOLYGON (((96.80559 3....
3     <NA> MULTIPOLYGON (((95.20544 5....
4     <NA> MULTIPOLYGON (((95.58431 4....
5     <NA> MULTIPOLYGON (((97.59461 2....
6     <NA> MULTIPOLYGON (((97.39178 2....
7     <NA> MULTIPOLYGON (((98.27612 4....
8     <NA> MULTIPOLYGON (((96.64762 4....
9     <NA> MULTIPOLYGON (((97.82461 3....
10    <NA> MULTIPOLYGON (((98.01772 4....

Kemudian, import data persentase kemiskinan di Provinsi DIY Tahun 2024, data yang digunakan diambil dari BPS

#import data tabular
pov_diy <- read_csv("poverty-2024.csv")
Rows: 552 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Nama Wilayah
dbl (1): Persentase

ℹ 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.
glimpse(pov_diy)
Rows: 552
Columns: 2
$ `Nama Wilayah` <chr> "ACEH", "Simeulue", "Aceh Singkil", "Aceh Selatan", "Ac…
$ Persentase     <dbl> 14.23, 17.69, 19.06, 12.02, 11.99, 13.26, 14.27, 17.60,…

Kita akan mengambil hanya data kabupaten/kota dari Provinsi DIY saya, sehingga gunakan code berikut untuk mengambil baris tersebut, kemudian melihat tipe datanya dengan fungsi glimpse()

diy_data <- pov_diy %>%
  slice(237:241)

glimpse(diy_data)
Rows: 5
Columns: 2
$ `Nama Wilayah` <chr> "Kulon Progo", "Bantul", "Gunung Kidul", "Sleman", "Kot…
$ Persentase     <dbl> 15.62, 11.66, 15.18, 7.46, 6.26

Melakukan filter kembali, kali ini kita akan filter hanya DIY saja yang diambil dari data spasial seluruh Indonesia

diy_kabkot <- kab_indo %>%
  filter(ADM1_EN=="Daerah Istimewa Yogyakarta")

Membuat peta untuk Provinsi DIY

diy_kabkot %>%
  ggplot()+
  geom_sf()+
  ggsflabel::geom_sf_label_repel(aes(label=str_wrap(ADM2_EN,width = 1)),size=1.75)+
  theme_bw()+
  theme(axis.title = element_blank())
Warning in st_point_on_surface.sfc(data$geometry): st_point_on_surface may not
give correct results for longitude/latitude data

st_crs(3857)$proj4string
[1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
st_crs(diy_kabkot)$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"
diy_kabkot_osm <- diy_kabkot %>%
  st_transform(3857)
st_crs(diy_kabkot_osm)$proj4string
[1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"

Setelah sudah membuat visualisasi untuk Prov di DIY, selanjutnya kita akan menggabungkan data spasial dan data persentase kemsikinan

diy_p0 <- diy_kabkot %>%
  full_join(y = diy_data,
            by = join_by(ADM2_EN==`Nama Wilayah`))

Melakukan visualisasi persentase penduduk miskin DIY

p0 <- diy_p0 %>%
  ggplot()+
  geom_sf(aes(fill=Persentase))+
  ggsflabel::geom_sf_label_repel(aes(label=str_wrap(ADM2_EN,width = 1)),size=1.75)+
  theme_bw()+
  theme(axis.title = element_blank())
p0
Warning in st_point_on_surface.sfc(data$geometry): st_point_on_surface may not
give correct results for longitude/latitude data

Menambahkan peta “asli” dari OpenStreetMap

p1 <- diy_p0 %>%
  st_transform(3857) %>%
  ggplot()+
  annotation_map_tile(type = "osm", zoomin = 0)+
  geom_sf(aes(fill=Persentase))+
  ggsflabel::geom_sf_label_repel(aes(label=str_wrap(ADM2_EN, width=1)),size=1.75)+
  scale_fill_stepsn(
    colours = viridis::magma(n = 6, direction = -1),
    breaks = seq(0, 20, 5),   # ✅ updated range
    limits = c(0, 20)          # ✅ force scale between 0 and 100
  )+
  annotation_north_arrow(
    location = "tl", which_north = "true",
    pad_x = unit(0.2, "cm"),
    pad_y = unit(0.2, "cm"),
    style = north_arrow_nautical, width = unit(1.5, "cm"),
    height = unit(1.5, "cm")
  ) +
  theme_bw()+
  theme(axis.title = element_blank())
p1
Zoom: 10

Lalu, bisa juga membuat visualisasi interaktifnya

diy_p0 %>%
  mapview::mapview(
    zcol = "Persentase",
    at = seq(0, 20, 1),   # ✅ updated range
    col.regions = viridis::magma(n = 6, direction = -1)
  )
Warning: Found less unique colors (6) than unique zcol values (21)! 
Interpolating color vector to match number of zcol values.