Tugas Sains Data Spasial P3

install.packages("viridis")

## Installing package into 'C:/Users/Hanung Safrizal/AppData/Local/R/win-library/4.3'
## (as 'lib' is unspecified)

## Warning in download.file(url, destfile, method, mode = "wb", ...): downloaded
## length 0 != reported length 0

## Warning in download.file(url, destfile, method, mode = "wb", ...): URL
## 'https://cloud.r-project.org/bin/windows/contrib/4.3/viridis_0.6.4.zip':
## Timeout of 60 seconds was reached

## Error in download.file(url, destfile, method, mode = "wb", ...) : 
##   download from 'https://cloud.r-project.org/bin/windows/contrib/4.3/viridis_0.6.4.zip' failed

## Warning in download.packages(pkgs, destdir = tmpd, available = available, :
## download of package 'viridis' failed

install.packages("mapview")

## Installing package into 'C:/Users/Hanung Safrizal/AppData/Local/R/win-library/4.3'
## (as 'lib' is unspecified)

## Warning in download.file(url, destfile, method, mode = "wb", ...): downloaded
## length 0 != reported length 0

## Warning in download.file(url, destfile, method, mode = "wb", ...): URL
## 'https://cloud.r-project.org/bin/windows/contrib/4.3/mapview_2.11.0.zip':
## Timeout of 60 seconds was reached

## Error in download.file(url, destfile, method, mode = "wb", ...) : 
##   download from 'https://cloud.r-project.org/bin/windows/contrib/4.3/mapview_2.11.0.zip' failed

## Warning in download.packages(pkgs, destdir = tmpd, available = available, :
## download of package 'mapview' failed

library(sf)

## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2

## ── 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)

Data Spasial dengan geometri poligon

Import data spasial

kabkot<- st_read("Admin2Kabupaten/idn_admbnda_adm2_bps_20200401.shp")

## Reading layer `idn_admbnda_adm2_bps_20200401' from data source 
##   `D:\Kuliah\S2\Semester 1\Sains Data\Responsi Sains Data_R Project\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

kabkot

## 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....

Import data tabular

jum_puskesmas <- read_csv("jml_puskesmas__jaringan_puskesmas_kabupatenkota_data.csv")

## Rows: 432 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): nama_provinsi, nama_kabupaten_kota, jenis_puskesmas, satuan
## dbl (5): id, kode_provinsi, kode_kabupaten_kota, jumlah_puskesmas, tahun
## 
## ℹ 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(jum_puskesmas)

## Rows: 432
## Columns: 9
## $ id                  <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,…
## $ kode_provinsi       <dbl> 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32…
## $ nama_provinsi       <chr> "JAWA BARAT", "JAWA BARAT", "JAWA BARAT", "JAWA BA…
## $ kode_kabupaten_kota <dbl> 3201, 3201, 3201, 3201, 3202, 3202, 3202, 3202, 32…
## $ nama_kabupaten_kota <chr> "KABUPATEN BOGOR", "KABUPATEN BOGOR", "KABUPATEN B…
## $ jenis_puskesmas     <chr> "PUSKESMAS RAWAT INAP", "PUSKESMAS NON RAWAT INAP"…
## $ jumlah_puskesmas    <dbl> 27, 74, 131, 122, 10, 48, 116, 118, 8, 37, 45, 116…
## $ satuan              <chr> "UNIT", "UNIT", "UNIT", "UNIT", "UNIT", "UNIT", "U…
## $ tahun               <dbl> 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 20…

jum_puskesmas %>% 
  count(tahun)

## # A tibble: 4 × 2
##   tahun     n
##   <dbl> <int>
## 1  2019   108
## 2  2020   108
## 3  2021   108
## 4  2022   108

jum_puskesmas2021 <- jum_puskesmas %>% 
                    filter(tahun==2021) %>%
                    select(nama_kabupaten_kota,jumlah_puskesmas)
jum_puskesmas2021

## # A tibble: 108 × 2
##    nama_kabupaten_kota jumlah_puskesmas
##    <chr>                          <dbl>
##  1 KABUPATEN BOGOR                   30
##  2 KABUPATEN BOGOR                   71
##  3 KABUPATEN BOGOR                  154
##  4 KABUPATEN BOGOR                  122
##  5 KABUPATEN SUKABUMI                10
##  6 KABUPATEN SUKABUMI                48
##  7 KABUPATEN SUKABUMI                79
##  8 KABUPATEN SUKABUMI               159
##  9 KABUPATEN CIANJUR                  8
## 10 KABUPATEN CIANJUR                 39
## # ℹ 98 more rows

Visualisasi Data Spasial

Visualisasi Statis

Sama seperti sebelumnya, visualisasi data spasial juga menggunakan package ggplot2 dan fungsi geom_sf.

ggplot(data = kabkot)+
  geom_sf()+
  theme_bw()

Hasil visualisasi diatas menujukkan peta Indonesia yang sudah terbagi menjadi poligon-poligon berdasarkan kabupaten/kota.

Untuk mendapatkan hanya provinsi Jawa Barat saja, maka kita bisa menggunakan fungsi filter dari package dplyr untuk menyaring kolom ADM1_EN yang berisi nama-nama provinsi di Indoneisa.

jabarkabkot <- kabkot %>%
                    filter(ADM1_EN=="Jawa Barat")
jabarkabkot

## Simple feature collection with 28 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.3703 ymin: -7.820979 xmax: 108.8469 ymax: -5.918148
## Geodetic CRS:  WGS 84
## First 10 features:
##    Shape_Leng Shape_Area       ADM2_EN ADM2_PCODE ADM2_REF ADM2ALT1EN
## 1    3.093005 0.14360665       Bandung     ID3204     <NA>       <NA>
## 2    3.027563 0.10447586 Bandung Barat     ID3217     <NA>       <NA>
## 3    2.553369 0.10359739        Bekasi     ID3216     <NA>       <NA>
## 4    4.450342 0.24474031         Bogor     ID3201     <NA>       <NA>
## 5    3.130958 0.13047729        Ciamis     ID3207     <NA>       <NA>
## 6    4.741142 0.29407915       Cianjur     ID3203     <NA>       <NA>
## 7    2.752827 0.08738231       Cirebon     ID3209     <NA>       <NA>
## 8    3.477542 0.25300160         Garut     ID3205     <NA>       <NA>
## 9    3.148833 0.17028961     Indramayu     ID3212     <NA>       <NA>
## 10   2.898903 0.15641686      Karawang     ID3215     <NA>       <NA>
##    ADM2ALT2EN    ADM1_EN ADM1_PCODE   ADM0_EN ADM0_PCODE       date    validOn
## 1        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 2        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 3        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 4        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 5        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 6        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 7        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 8        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 9        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 10       <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
##    validTo                       geometry
## 1     <NA> MULTIPOLYGON (((107.7331 -6...
## 2     <NA> MULTIPOLYGON (((107.4095 -6...
## 3     <NA> MULTIPOLYGON (((107.0757 -5...
## 4     <NA> MULTIPOLYGON (((106.9708 -6...
## 5     <NA> MULTIPOLYGON (((108.36 -7.0...
## 6     <NA> MULTIPOLYGON (((107.2302 -6...
## 7     <NA> MULTIPOLYGON (((108.685 -6....
## 8     <NA> MULTIPOLYGON (((107.9182 -6...
## 9     <NA> MULTIPOLYGON (((108.2003 -6...
## 10    <NA> MULTIPOLYGON (((107.1126 -5...

jabarkabkot %>% 
  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

• Fungsi geom_sf_label_repel dari package ggsflabel berguna untuk memberikan label-label kabupaten/kota. Karena ada kata repel difungsi tersebut, menyebabkan label yang dihasilkan tidak akan tumpang tindih. Label yang berpotensi tumpang-tindih satu sama lain akan digeser kemudian diberi garis yang menandakan asal dari label tersebut.
• Fungsi str_wrap digunakan untuk memindahkan kata ke baris baru jika kata-kata tersebut terdiri dari dua kata seperti “Kota Bogor”.

kita juga bisa menambahkan peta “asli” dari OpenStreetMap dengan fungsi annotation_map_tile.

jabarkabkot %>% 
  ggplot()+
    annotation_map_tile(type = "osm", zoomin = 0)+
    geom_sf()+
    ggsflabel::geom_sf_label_repel(aes(label=str_wrap(ADM2_EN,width = 1)),size=1.75)+
    theme_bw()

## Warning in st_point_on_surface.sfc(data$geometry): st_point_on_surface may not
## give correct results for longitude/latitude data

## Loading required namespace: raster

## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)

## Zoom: 9

Kemudian, langkah selanjutnya kita akan memasukan informasi jumlah puskesmas tahun 2021 ke dalam peta diatas. Pertama kita harus memastikan bahwa nama-nama kabupaten/kota di data jumlah tenaga kesehatan dan di dalam data spasial kita (yang diimport dari .shp) itu sama. Berikut adalah hasilnya

jum_puskesmas2021 %>% count(nama_kabupaten_kota)

## # A tibble: 27 × 2
##    nama_kabupaten_kota         n
##    <chr>                   <int>
##  1 KABUPATEN BANDUNG           4
##  2 KABUPATEN BANDUNG BARAT     4
##  3 KABUPATEN BEKASI            4
##  4 KABUPATEN BOGOR             4
##  5 KABUPATEN CIAMIS            4
##  6 KABUPATEN CIANJUR           4
##  7 KABUPATEN CIREBON           4
##  8 KABUPATEN GARUT             4
##  9 KABUPATEN INDRAMAYU         4
## 10 KABUPATEN KARAWANG          4
## # ℹ 17 more rows

kabkot %>% as_tibble() %>% count(ADM2_EN)

## # A tibble: 519 × 2
##    ADM2_EN             n
##    <chr>           <int>
##  1 Aceh Barat          1
##  2 Aceh Barat Daya     1
##  3 Aceh Besar          1
##  4 Aceh Jaya           1
##  5 Aceh Selatan        1
##  6 Aceh Singkil        1
##  7 Aceh Tamiang        1
##  8 Aceh Tengah         1
##  9 Aceh Tenggara       1
## 10 Aceh Timur          1
## # ℹ 509 more rows

Berdasarkan output tersebut, dapat disimpulkan bahwa nama kabupaten/kota dari data spasial dan data jumlah puskesmas tidak sama, sehingga kita perlu memodifikasi terlebih dahulu supaya sama. Dalam kasus ini kita akan memodifikasi nama kabupaten/kota dari data jumlah puskesmas.

library(dplyr)

jum_puskesmas2021fin <- jum_puskesmas2021 %>%
  mutate(nama_kabupaten_kota = str_remove(nama_kabupaten_kota, "KABUPATEN ") %>% 
                                  str_to_title()) %>%
  group_by(nama_kabupaten_kota) %>%
  summarise(jumlah_puskesmas = sum(jumlah_puskesmas))

jum_puskesmas2021fin

## # A tibble: 27 × 2
##    nama_kabupaten_kota jumlah_puskesmas
##    <chr>                          <dbl>
##  1 Bandung                          188
##  2 Bandung Barat                    100
##  3 Bekasi                           142
##  4 Bogor                            377
##  5 Ciamis                           137
##  6 Cianjur                          163
##  7 Cirebon                          214
##  8 Garut                            290
##  9 Indramayu                        165
## 10 Karawang                         172
## # ℹ 17 more rows

Setelah nama kabupaten/kota sama, langkah selanjutnya adalah menggabungkan data spasial dan data jumlah puskesmas. Proses penggabungan ini bisa menggunakan join dari package dplyr. Dalam kasus ini kita coba menggunakan full_join untuk mengidentifikasi apakah nama kabupaten/kota di kedua sumber sudah sama.

jabar_puskesmas <- jabarkabkot %>% 
                  full_join(y = jum_puskesmas2021fin,
                            by = join_by(ADM2_EN==nama_kabupaten_kota))
jabar_puskesmas

## Simple feature collection with 28 features and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 106.3703 ymin: -7.820979 xmax: 108.8469 ymax: -5.918148
## Geodetic CRS:  WGS 84
## First 10 features:
##    Shape_Leng Shape_Area       ADM2_EN ADM2_PCODE ADM2_REF ADM2ALT1EN
## 1    3.093005 0.14360665       Bandung     ID3204     <NA>       <NA>
## 2    3.027563 0.10447586 Bandung Barat     ID3217     <NA>       <NA>
## 3    2.553369 0.10359739        Bekasi     ID3216     <NA>       <NA>
## 4    4.450342 0.24474031         Bogor     ID3201     <NA>       <NA>
## 5    3.130958 0.13047729        Ciamis     ID3207     <NA>       <NA>
## 6    4.741142 0.29407915       Cianjur     ID3203     <NA>       <NA>
## 7    2.752827 0.08738231       Cirebon     ID3209     <NA>       <NA>
## 8    3.477542 0.25300160         Garut     ID3205     <NA>       <NA>
## 9    3.148833 0.17028961     Indramayu     ID3212     <NA>       <NA>
## 10   2.898903 0.15641686      Karawang     ID3215     <NA>       <NA>
##    ADM2ALT2EN    ADM1_EN ADM1_PCODE   ADM0_EN ADM0_PCODE       date    validOn
## 1        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 2        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 3        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 4        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 5        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 6        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 7        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 8        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 9        <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
## 10       <NA> Jawa Barat       ID32 Indonesia         ID 2019-12-20 2020-04-01
##    validTo jumlah_puskesmas                       geometry
## 1     <NA>              188 MULTIPOLYGON (((107.7331 -6...
## 2     <NA>              100 MULTIPOLYGON (((107.4095 -6...
## 3     <NA>              142 MULTIPOLYGON (((107.0757 -5...
## 4     <NA>              377 MULTIPOLYGON (((106.9708 -6...
## 5     <NA>              137 MULTIPOLYGON (((108.36 -7.0...
## 6     <NA>              163 MULTIPOLYGON (((107.2302 -6...
## 7     <NA>              214 MULTIPOLYGON (((108.685 -6....
## 8     <NA>              290 MULTIPOLYGON (((107.9182 -6...
## 9     <NA>              165 MULTIPOLYGON (((108.2003 -6...
## 10    <NA>              172 MULTIPOLYGON (((107.1126 -5...

p0 <- jabar_puskesmas %>% 
  ggplot()+
  geom_sf(aes(fill=jumlah_puskesmas))+
  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,300,50))+
#colours =viridis::magma(n=10,direction = -1 ),nbreaks=5) +
  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

- Argumen fill=jumlah_puskesmas dalam geom_sf(aes(fill=jumlah_puskesmas)) menunjukkan bahwa informasi banyaknya tenaga kesehatan masyarakat akan ditampilkan dalam bentuk warna-warna. - fungsi scale_fill_stepsn akan memberi warna pada angka-angka pada yang ada dalam kolom jumlah_nakesmas. Sebelumnya angka-angka tersebut harus kita ubah ke dalam interval-interval. Dalam kasus ini intervalnya adalah

seq(0,300,50)

## [1]   0  50 100 150 200 250 300

• fungsi magma() dari package viridis merupakan fungsi untuk memanggil color pallete. Untuk keterangan lebih lanjut silahkan help dari fungsi magma help(viridis::magma).
• Argumen axis.title = element_blank() dari fungsi theme berarti kita ingin menghilangkan judul kolom dari plot.

kita juga bisa menambahkan peta “asli” dari OpenStreetMap dengan fungsi annotation_map_tile.

p1 <- jabar_puskesmas %>% 
  ggplot()+
  annotation_map_tile(type = "osm", zoomin = 0)+
  geom_sf(aes(fill=jumlah_puskesmas))+
  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,300,50))+
#colours =viridis::magma(n=10,direction = -1 ),nbreaks=5) +
  theme_bw()+
  theme(axis.title = element_blank())
p1

## Warning in st_point_on_surface.sfc(data$geometry): st_point_on_surface may not
## give correct results for longitude/latitude data

## Zoom: 9

Visualisasi Interaktif

jabar_puskesmas %>% 
  mapview::mapview(zcol="jumlah_puskesmas",
                   at = seq(0,300,50),
                   col.regions=viridis::magma(n = 7,direction = -1))