Praktikum 3 - Data Spasial Menggunakan R

options(repos = c(CRAN = "https://cloud.r-project.org/"))

remotes::install_github("yutannihilation/ggsflabel")

## Skipping install of 'ggsflabel' from a github remote, the SHA1 (846ede01) has not changed since last install.
##   Use `force = TRUE` to force installation

install.packages("viridis")

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

## package 'viridis' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Hanung Safrizal\AppData\Local\Temp\RtmpKApP8X\downloaded_packages

install.packages("mapview")

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

## package 'mapview' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Hanung Safrizal\AppData\Local\Temp\RtmpKApP8X\downloaded_packages

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 titik

Import data spasial

CholeraDeaths <- st_read(dsn = "SnowGIS_SHP",layer = "Cholera_Deaths")

## Reading layer `Cholera_Deaths' from data source 
##   `D:\Kuliah\S2\Sains Data\Responsi Sains Data_R Project\SnowGIS_SHP' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 250 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 529160.3 ymin: 180857.9 xmax: 529655.9 ymax: 181306.2
## Projected CRS: OSGB36 / British National Grid

CholeraDeaths <- st_read(dsn = "SnowGIS_SHP/Cholera_Deaths.shp")

## Reading layer `Cholera_Deaths' from data source 
##   `D:\Kuliah\S2\Sains Data\Responsi Sains Data_R Project\SnowGIS_SHP\Cholera_Deaths.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 250 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 529160.3 ymin: 180857.9 xmax: 529655.9 ymax: 181306.2
## Projected CRS: OSGB36 / British National Grid

CholeraDeaths

Visualisasi data spasial

Visualisasi statis

ggplot(data = CholeraDeaths)+
  geom_sf()+
  #mengubah tema
  theme_bw()

ggplot(CholeraDeaths) + 
  annotation_map_tile(type = "osm", zoomin = 0) + 
  geom_sf(aes(size = Count), alpha = 0.7)+
  theme_bw()

## 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: 17

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(3857)$wkt

## [1] "PROJCRS[\"WGS 84 / Pseudo-Mercator\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"Popular Visualisation Pseudo-Mercator\",\n        METHOD[\"Popular Visualisation Pseudo Mercator\",\n            ID[\"EPSG\",1024]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"False easting\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting (X)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"northing (Y)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Web mapping and visualisation.\"],\n        AREA[\"World between 85.06°S and 85.06°N.\"],\n        BBOX[-85.06,-180,85.06,180]],\n    ID[\"EPSG\",3857]]"

st_crs(CholeraDeaths)$proj4string

## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"

st_crs(CholeraDeaths)$wkt

## [1] "PROJCRS[\"OSGB36 / British National Grid\",\n    BASEGEOGCRS[\"OSGB36\",\n        DATUM[\"Ordnance Survey of Great Britain 1936\",\n            ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4277]],\n    CONVERSION[\"British National Grid\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",49,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-2,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",400000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",-100000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.\"],\n        BBOX[49.75,-9.01,61.01,2.01]],\n    ID[\"EPSG\",27700]]"

cholera_latlong <- CholeraDeaths %>%
  st_transform(3857)

st_crs(cholera_latlong)$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(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"

ggplot(cholera_latlong) + 
  annotation_map_tile(type = "osm", zoomin = 0) + 
  geom_sf(aes(size = Count),color="#0288D1", alpha = 0.7)+
  theme_bw()

## Zoom: 17

Visualisasi Interaktif

Visualisasi spasial interaktif bisa dilakukan dengan menggunakan fungsi mapview package mapview.

cholera_latlong %>% 
  mapview::mapview(zcol = "Count")

Data Spasial dengan geometri poligon

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:\Kuliah\S2\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

kab_indo

Import data tabular

jum_nakes <- read_csv("jml_tenaga_kesehatan_msyrkt__kabupatenkota.csv")

## Rows: 162 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): nama_provinsi, nama_kabupaten_kota, satuan
## dbl (5): id, kode_provinsi, kode_kabupaten_kota, jumlah_nakesmas, 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_nakes)

## Rows: 162
## Columns: 8
## $ 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, 3202, 3203, 3204, 3205, 3206, 3207, 3208, 32…
## $ nama_kabupaten_kota <chr> "KABUPATEN BOGOR", "KABUPATEN SUKABUMI", "KABUPATE…
## $ jumlah_nakesmas     <dbl> 156, 62, 19, 9, 53, 13, 24, 25, 58, 40, 41, 62, 38…
## $ satuan              <chr> "ORANG", "ORANG", "ORANG", "ORANG", "ORANG", "ORAN…
## $ tahun               <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 20…

jum_nakes %>% 
  count(tahun)

jum_nakes2021 <- jum_nakes %>% 
                    filter(tahun==2021) %>%
                    select(nama_kabupaten_kota,jumlah_nakesmas)
jum_nakes2021

Visualisasi Data Spasial

Visualisasi Statis

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

ggplot(data = kab_indo)+
  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.

jabar_kabkot <- kab_indo %>%
                    filter(ADM1_EN=="Jawa Barat")
jabar_kabkot

jabar_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

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

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

## Zoom: 9

Kemudian, langkah selanjutnya kita akan memasukan informasi jumlah tenaga kesehatan 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_nakes2021 %>% count(nama_kabupaten_kota)

kab_indo %>% as_tibble() %>% count(ADM2_EN)

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

jum_nakes2021fin <- jum_nakes2021 %>% 
                    mutate(nama_kabupaten_kota= nama_kabupaten_kota %>% 
                             #dibuang kata "KABUPATEN"
                             str_remove(pattern = "KABUPATEN ") %>% 
                             #diubah menjadi huruf kapital setiap awal kata
                             str_to_title()
                           )
jum_nakes2021fin

Setelah nama kabupaten/kota sama, langkah selanjutnya adalah menggabungkan data spasial dan data jumlah tenaga kesehatan. 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_nakes <- jabar_kabkot %>% 
                  full_join(y = jum_nakes2021fin,
                            by = join_by(ADM2_EN==nama_kabupaten_kota))
jabar_nakes

p0 <- jabar_nakes %>% 
  ggplot()+
  geom_sf(aes(fill=jumlah_nakesmas))+
  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_nakesmas dalam geom_sf(aes(fill=jumlah_nakesmas)) 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_nakes %>% 
  ggplot()+
  annotation_map_tile(type = "osm", zoomin = 0)+
  geom_sf(aes(fill=jumlah_nakesmas))+
  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_nakes %>% 
  mapview::mapview(zcol="jumlah_nakesmas",
                   at = seq(0,300,50),
                   col.regions=viridis::magma(n = 7,direction = -1))