Pengenalan Geospasial dengan Simple Features `sf`

Geospasial

R Programming

Data Visualization

Author

Gerry Alfa Dito

Published

August 17, 2023

Pengantar Data Spasial

Data spasial adalah data yang memberi informasi tentang di mana letak (lokasi) suatu objek. Informasi yang diberikan dapat berupa lokasi geografis atau posisi objek tertentu di permukaan Bumi. Data Spasial ini dapat berbentuk titik-titik (points), garis-garis (lines), poligon-poligon (polygons), atau bahkan citra-citra (images). Dalam tulisan ini, kita akan hanya akan mengenalkan data spasial sebagai data yang berisi informasi lokasi geografis.

Sebagai ilustrasi, Di peta, kita dapat melihat lokasi kota-kota, sungai-sungai, gunung-gunung, dan objek-objek lainnya. Setiap elemen ini mewakili data spasial. Namun, data spasial tidak hanya terbatas pada peta. Setiap informasi yang terhubung dengan suatu tempat tertentu bisa dianggap sebagai data spasial. Hal ini bisa mencakup hal-hal seperti pengukuran suhu di lokasi berbeda, distribusi pohon di hutan, atau jumlah penduduk di berbagai lingkungan.

Data spasial sangat penting karena membantu kita memahami bagaimana objek-objek tersebar di suatu ruang dan bagaimana interaksi di antara objek-objek tersebut berdasarkan lokasi objek. Dengan menganalisis data spasial, kita dapat mengungkap pola dan hubungan yang mungkin tidak terlihat dari data biasa. Di sinilah statistik spasial berperan, yakni memberikan alat dan teknik untuk menganalisis dan memahami jenis data ini.

Coordinate Reference Systems (CRS)

Coordinate Reference System, sering disingkat sebagai CRS, adalah cara untuk menentukan bagaimana koordinat yakni lintang (latitude) dan bujur (longitude) diberikan kepada lokasi objek di permukaan bumi. Kita dapat memandang CRS sebagai kumpulan aturan matematis yang memungkinkan kita untuk dengan tepat mendefinisikan di mana lokasi suatu objek terletak di peta atau dalam ruang geografis.

Bumi itu bulat, seperti bola, tetapi peta yang digunakan berupa permukaan datar. Ketika kita merepresentasikan permukaan Bumi yang melengkung pada peta datar, terjadi sedikit distorsi. CRS membantu kita mengelola distorsi ini dan memastikan bahwa peta yang terbentuk mewakili lokasi-lokasi di dunia nyata secara akurat.

CRS mencakup dua komponen utama: datum dan proyeksi.

datum adalah sebuah parameter atau sekumpulan parameter yang menentukan posisi asal, skala, dan orientasi sistem koordinat.
Proyeksi adalah metode yang digunakan untuk meratakan permukaan Bumi ke dalam peta 2D. Proyeksi yang berbeda digunakan untuk tujuan yang berbeda, karena setiap proyeksi mengubah permukaan Bumi dengan cara yang unik.

Ada dua jenis utama datum yaitu:

Datum Geografis: Datum geografis mendefinisikan pusat, orientasi, dan skala bentuk elipsoid permukaan Bumi. Ini seperti sistem koordinat yang menempatkan titik-titik di permukaan Bumi. Datum ini membantu menentukan koordinat lintang dan bujur. Jenis-jenis datum Geografis antara lain:

WGS84 (World Geodetic System 1984): Ini adalah datum yang banyak digunakan untuk sistem GPS dan aplikasi pemetaan. Ini menyediakan kerangka referensi global dan umumnya digunakan untuk navigasi satelit.
NAD83 (North American Datum 1983): Datum ini digunakan di Amerika Utara dan menyediakan kerangka kerja untuk pemetaan dan survei di Amerika Serikat, Kanada, dan Meksiko.
ETRS89 (European Terrestrial Reference System 1989): ETRS89 adalah datum geosentris yang digunakan di Eropa untuk berbagai tujuan pemetaan dan penentuan posisi.
GDA94 (Geocentric Datum of Australia 1994): Datum ini digunakan untuk tujuan pemetaan dan survei di Australia dan wilayah-wilayahnya.
Tokyo Datum 2000: Digunakan di Jepang, datum ini menyediakan kerangka referensi untuk pemetaan dan survei di kepulauan Jepang.
Indian 1975: Digunakan di India, datum ini berfungsi sebagai referensi untuk pemetaan dan kegiatan geodetik di negara tersebut.
HARN (High Accuracy Reference Network): Ini adalah sekelompok datum yang digunakan di Amerika Serikat untuk meningkatkan akurasi pemetaan dan survei di wilayah-wilayah tertentu.
Beijing 1954: Digunakan di Tiongkok, datum ini berfungsi sebagai referensi untuk pemetaan dan pekerjaan geodetik di negara tersebut.

Datum Vertikal: Sementara datum geografis fokus pada permukaan Bumi, datum vertikal digunakan untuk mengukur ketinggian dan kedalaman, terutama dalam hubungannya dengan permukaan laut. Mereka mendefinisikan titik awal untuk mengukur elevasi. Jenis-jenis datum Vertikal antara lain:

NGVD29 (National Geodetic Vertical Datum of 1929): Datum ini banyak digunakan di Amerika Serikat dan menjadi acuan untuk mengukur ketinggian berdasarkan titik ukur pasang surut tertentu.
NAVD88 (North American Vertical Datum of 1988): Versi yang diperbarui dari NGVD29, NAVD88 adalah datum vertikal standar saat ini di Amerika Utara. Ini menggunakan jaringan titik ukur dan teknologi canggih untuk memberikan pengukuran ketinggian yang akurat.
EGM96 (Earth Gravitational Model 1996): Meskipun bukan datum vertikal tradisional, EGM96 adalah model matematika yang mendefinisikan medan gravitasi Bumi. Sering digunakan bersama dengan model geoid untuk menentukan ketinggian ortometrik.
AHD (Australian Height Datum): Digunakan di Australia, AHD memberikan acuan untuk mengukur ketinggian di berbagai lanskap yang beragam di negara tersebut.
GEOID12B: Model geoid yang dikembangkan oleh National Geodetic Survey di Amerika Serikat untuk memberikan ketinggian ortometrik yang tepat menggunakan EGM96.
CVD2014 (Canadian Vertical Datum 2014): Datum ini digunakan di Kanada dan memberikan referensi yang diperbarui untuk pengukuran ketinggian.

Transformasi antar datum dalam CRS ini dapat menggunakan PROJ4 dan WKT-4.

Objek Geometri data spasial

type	description
POINT	single point geometry
MULTIPOINT	set of points
LINESTRING	single linestring (two or more points connected by straight
lines)
MULTILINESTRING	set of linestrings
POLYGON	exterior ring with zero or more inner rings, denoting holes
set of polygons
MULTIPOLYGON	set of points
GEOMETRYCOLLECTION	set of the geometries above

Data

Data yang digunakan untuk praktikum kali ini adalah

Cholera Death
Peta Administratif tingkat kabupaten/kota seluruh indonesia (Update 2020)
Jumlah tenaga kesehatan masyarakat di Jawa Barat Tahun 2016-2021

Ketiga dataset ini dapat didownload melalui link berikut Download Dataset

Silahkan extract zip terlebih dahulu sebelum menggunakan data-data tersebut.

Cholera Death

John Snow produced a famous map in 1854 showing the deaths caused by a cholera outbreak in Soho, London, and the locations of water pumps in the area. By doing this he found there was a significant clustering of the deaths around a certain pump – and removing the handle of the pump stopped the outbreak.

Package

Silahkan install package berikut dari github.

Catatan: jika saat menjalankan sintaks ada output

These packages have more recent versions available. It is recommended to update all of them. Which would you like to update? * 1: All
* 2: CRAN packages only
* 3: None

Silahkan pilih None

remotes::install_github("yutannihilation/ggsflabel")

Silahkan install package berikut tanpa perlu dipanggil

install.packages("viridis")
install.packages("mapview")

library(sf)

Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.4.4     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ 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 dalam R

Di dalam R terdapat package-package yang tergabung dalam ekosistem untu menangani data spasial. Secara umum, package dasar yang dapat menangani data spasial adalah package sp dan sf. Dalam ilustrasi kali ini kita akan berfokus pada penggunaan package sf.

Package `sf`

Package sf dimaksudkan untuk menggantikan package R sp, rgeos dan bagian vektor dari rgdal, package R sf (Pebesma 2018) dikembangkan untuk memindahkan analisis data spasial dalam R yang lebih dekat dengan pendekatan berbasis standar dalam industri dan proyek open-source, untuk membangun versi yang lebih modern dari software geospasial open-source, dan untuk memungkinkan integrasi software spasial R dengan tidyverse (Wickham dkk. 2019).

Untuk melakukannya, package R sf menyediakan akses simple features (Herring et al. 2011), secara native, ke R. Package ini dapat dihubungkan ke beberapa package tidyverse, khususnya ke ggplot2, dplyr, dan tidyr. Package ini dapat membaca dan menyimpan data melalui GDAL, menjalankan operasi geometri menggunakan GEOS (untuk koordinat yang diproyeksikan) atau s2geometry (untuk koordinat ellipsoidal), dan menjalankan transformasi koordinat atau konversi menggunakan PROJ. Library C++ eksternal dihubungkan dengan package R Rcpp (Eddelbuettel 2013).

Paket sf merepresentasikan kumpulan simple features dalam objek sf, sebuah sub-kelas dari data.frame atau tibble. Objek sf berisi setidaknya satu geometri list-column dari kelas sfc, yang mana untuk setiap elemennya berisi geometri sebagai objek R dari kelas sfg. Geometri list-column bertindak sebagai sebuah variabel dalam data.frame atau tibble, tetapi memiliki struktur yang lebih kompleks daripada vektor dasar seperti variabel numerik atau karakter.

Sebuah objek sf memiliki metadata sebagai berikut: - nama kolom geometri (aktif), yang disimpan dalam atribut sf_column - untuk setiap variabel non-geometri, hubungan atribut-geometri yang disimpan dalam atribut agr

Suatu sfc geometri list-column diekstrak dari objek sf dengan st_geometry dan memiliki metadata sebagai berikut: - sistem referensi koordinat yang disimpan dalam atribut crs - kotak pembatas yang disimpan dalam atribut bbox - presisi disimpan dalam atribut presisi - jumlah geometri kosong yang disimpan dalam atribut n_empty

Atribut-atribut ini paling mudah diakses atau diatur dengan menggunakan fungsi-fungsi seperti st_bbox, st_crs, st_set_crs, st_agr, st_set_agr, st_precision, dan st_set_precision.Kolom geometri pada objek sf dapat diatur atau diganti dengan menggunakan st_geometry<- atau st_set_geometry.

Data Spasial dengan geometri titik

Import data spasial

Data spasial dapat berupa ekstensi .shp atau dbf. Untuk mengimport data spasial bisa menggunakan fungsi st_read dari package sf.

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

Reading layer `Cholera_Deaths' from data source 
  `D:\Campus\MyWebsite\Quarto Web\gerrydito.github.io\posts\2023-8-17-Pengenalan-Geospasial-dengan-sf\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

argumen layer pada sintaks diatas adalah nama file .shp atau .dbf yang akan diimport. Sementara itu, dsn adalah nama folder tempat file .shp atau .dbf berada. Selain dengan cara diatas, cara lain untuk import data spasial menggunakan sintaks st_read langsung menuliskan direktori atau nama file .shp atau .dbf. Berikut ilustrasinya.

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

Reading layer `Cholera_Deaths' from data source 
  `D:\Campus\MyWebsite\Quarto Web\gerrydito.github.io\posts\2023-8-17-Pengenalan-Geospasial-dengan-sf\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

Setelah proses import selesai, maka akan tampil beberapa informasi tentang data spasial yang diimport seperti

250 Features dan 2 fields berarti terdapat 250 objek spasial dan 2 kolom data yang bersesuaian dengan objek spasial tersebut.
Tipe Geometri pada data ini adalah POINT yang berarti bahwa data spasial yang diimport berupa titik-titik amatan yang tersebar dalam ruang spasial.
Bounding box yang merupakan batas minimum dan maximum di latitude (sumbu \(x\)) atau longitude (sumbu \(y\)).
CRS yang digunakan dalam data ini adalah British National Grid

Kemudian jika kita panggil objek CholeraDeaths akan muncul overview objek-objek spasial (baris) sebagai berikut:

CholeraDeaths

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
First 10 features:
   Id Count                  geometry
1   0     3 POINT (529308.7 181031.4)
2   0     2 POINT (529312.2 181025.2)
3   0     1 POINT (529314.4 181020.3)
4   0     1 POINT (529317.4 181014.3)
5   0     4 POINT (529320.7 181007.9)
6   0     2   POINT (529336.7 181006)
7   0     2 POINT (529290.1 181024.4)
8   0     2   POINT (529301 181021.2)
9   0     3   POINT (529285 181020.2)
10  0     2 POINT (529288.4 181031.8)

Visualisasi data spasial

Visualisasi statis

Kemudian kita bisa mevisualisasikan data spasial diatas dengan bantuan package ggplot2, secara khusus dengan memanfaatkan fungsi geom_sf.

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

Grafik diatas sulit diinterpretasikan karena hanya berupa titik-titik spasial saja. Kita dapat menambahkan peta yang bersesuaian dengan menambahkan fungsi annotation_map_tile dari package ggspatial. Fungsi annotation_map_tile mengcapture peta yang berasal dari OpenStreetMap sesuai dengan posisi objek spasial. Kita dapat mengontrol tingkat zoom, serta type-nya.

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

Loading required namespace: raster

Zoom: 17

Hal menarik lain yang terlihat dari map diatas adalah ada beberapa titik yang masuk dalam map tiles. Hal ini dikarenakan adanya perbedaan CRS pada data Cholera Death dan map tiles yang berasal dari Open Street Map. CRS pada data Cholera Death adalah OSGB36 / British National Grid. Semetara itu, CRS pada Open Street Map atau maps yang lain adalah WGS84 atau sering juga dikenal sebagai EPSG:3857. Fungsi st_crs() akan menerjemahkan dari kode EPSG ke string PROJ.4 dan WKT. Berikut sintaks untuk menampilkan string PROJ.4 dan WKT

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]]"

Jika kita bandingkan dengan PROJ.4 dan WKT dari CholeraDeaths

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]]"

Kemudian dengan st_transform kita melakukan proyeksi (transformasi) data CholeraDeaths dari CRS OSGB36 ke EPSG:3857

cholera_latlong <- CholeraDeaths %>%
  st_transform(3857)

karena WKT agak sulit untuk dibandingkan kita menggunakan Proj4 saja

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"

PROJ4 untuk data CholeraDeaths hasil transformasi sudah sama dengan EPSG:3857.

Kemudian jika kita buatkan kembali map, maka hasilnya adalah

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")

argumen zcol digunakan untuk menampilkan informasi tertentu dalam bentuk warna.

Data Spasial dengan geometri poligon

Untuk ilustrasi data spasial yang berbentuk poligon, kita akan mencoba melakukan visualisasi spasial Jumlah tenaga kesehatan masyarakat setiap Kabupaten/kota di Jawa Barat tahun 2021.

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:\Campus\MyWebsite\Quarto Web\gerrydito.github.io\posts\2023-8-17-Pengenalan-Geospasial-dengan-sf\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

Setelah proses import selesai, maka akan tampil beberapa informasi tentang data spasial yang diimport seperti

522 Features dan 14 fields berarti terdapat 522 objek spasial (poligon) dan 14 kolom data yang bersesuaian dengan objek spasial tersebut. Dalam kasus ini banyaknya poligon adalah sebanyak kabupaten/kota di seluruh indonesia.
Tipe Geometri pada data ini adalah MULTIPOLYGON.
Bounding box yang merupakan batas minimum dan maximum di latitude (sumbu \(x\)) atau longitude (sumbu \(y\)).
CRS yang digunakan dalam data ini adalah 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....

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…

Data tabular yang kita import memliki informasi lebih dari satu tahun. Sehingga kita perlu filter terlebih dahulu untuk Tahun 2021.

jum_nakes %>% 
  count(tahun)

# A tibble: 6 × 2
  tahun     n
  <dbl> <int>
1  2016    27
2  2017    27
3  2018    27
4  2019    27
5  2020    27
6  2021    27

Berikut ilustrasi filternya:

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

# A tibble: 27 × 2
   nama_kabupaten_kota   jumlah_nakesmas
   <chr>                           <dbl>
 1 KABUPATEN BOGOR                   291
 2 KABUPATEN SUKABUMI                 74
 3 KABUPATEN CIANJUR                  66
 4 KABUPATEN BANDUNG                 135
 5 KABUPATEN GARUT                   134
 6 KABUPATEN TASIKMALAYA               1
 7 KABUPATEN CIAMIS                  105
 8 KABUPATEN KUNINGAN                 86
 9 KABUPATEN CIREBON                  97
10 KABUPATEN MAJALENGKA               40
# ℹ 17 more rows

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

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

Kemudian, peta Jawa Barat bisa kita dapatkan

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

Fetching 16 missing tiles


  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |======================================================================| 100%

...complete!

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)

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

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

# A tibble: 27 × 2
   nama_kabupaten_kota jumlah_nakesmas
   <chr>                         <dbl>
 1 Bogor                           291
 2 Sukabumi                         74
 3 Cianjur                          66
 4 Bandung                         135
 5 Garut                           134
 6 Tasikmalaya                       1
 7 Ciamis                          105
 8 Kuningan                         86
 9 Cirebon                          97
10 Majalengka                       40
# ℹ 17 more rows

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

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_nakesmas                       geometry
1     <NA>             135 MULTIPOLYGON (((107.7331 -6...
2     <NA>              50 MULTIPOLYGON (((107.4095 -6...
3     <NA>              77 MULTIPOLYGON (((107.0757 -5...
4     <NA>             291 MULTIPOLYGON (((106.9708 -6...
5     <NA>             105 MULTIPOLYGON (((108.36 -7.0...
6     <NA>              66 MULTIPOLYGON (((107.2302 -6...
7     <NA>              97 MULTIPOLYGON (((108.685 -6....
8     <NA>             134 MULTIPOLYGON (((107.9182 -6...
9     <NA>             102 MULTIPOLYGON (((108.2003 -6...
10    <NA>              44 MULTIPOLYGON (((107.1126 -5...

Karena banyaknya features (baris) sudah sama seperti sebelum digabung, maka hasil penggabungan sudah sesuai. Kemudian kita akan langsung membuat visualisasinya

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

argumen zcol digunakan untuk menampilkan informasi tertentu dalam bentuk warna.
at digunakan untuk membagi angka-angka yang yang akan dimasukan ke argumen zcol ke dalam beberapa interval.

Pengantar Data Spasial

Coordinate Reference Systems (CRS)

Objek Geometri data spasial

Data

Cholera Death

Package

Data Spasial dalam R

Package sf

Data Spasial dengan geometri titik

Import data spasial

Visualisasi data spasial

Visualisasi statis

Visualisasi Interaktif

Data Spasial dengan geometri poligon

Import data spasial

Import data tabular

Visualisasi Data Spasial

Visualisasi Statis

Visualisasi Interaktif

Package `sf`