R Basic for Spatial Data

MDM Manessa

Table of Contents

1. Basic of R

1.1 Introduction

1.1.1 R Environment

1.1.2 R Package

1.2 Basic of R

1.2.1 Basic Math

1.2.2 Variable

1.2.3 Calling Function

1.2.4 Data type

1.3 Logical function

1.3.1 if else

1.3.2 for loop

1.4 Data Visualizating

1.4.1 Base Graphics

1.4.2 ggplot2

2.Spatial data in R

Prerequisites

2.1 Introduction

2.2 Vector data

2.3 Raster data

2.4 Coordinate Reference System

3. Spatial Function

Prerequisites

3.1 Spatial operation on vector data

3.1.1 Spatial subsetting

3.1.2 Topological relations

3.1.3 Spatial joining

3.1.4 Spatial data aggregation

3.1.5 Spatial distance

3.2 Spatial operation on raster data

3.2.1 Spatial subsetting

3.2.2 Map algebra

3.2.3 Local operations

3.2.4 Focal operations

3.2.5 Zonal operation

3.2.6 Global operation

3.2.7 Aggregation and disaggregation

3.5 Raster-vector interactions

3.5.1 Raster cropping

3.5.2 Raster extraction

3.5.3 Rasterization

3.5.4 Spatial vectorization

3.6 Reproject Spatial Data

3.7 Geographic data I/O

3.7.1 Retrieving open data

3.7.2 Geographic data packages

3.7.3 Geographic web services

3.7.4 Data input (I)

3.7.5 Data output (0)

4. Making maps in R

Prerequisites

4.1 Introduction

4.2 Static maps

4.2.1 tmap basics

4.2.2 Map objects

4.2.3 Aesthetics

4.2.4 Color settings

4.2.5 Layouts

4.2.6 Faceted maps

4.3 Interactive maps

4.3.1 tmap

4.3.2 leaflet

4.4 Mapping applications

4.5 Other mapping packages

4.5.1 ggplot2

4.5.2 cartogram

1. Basic of R

1.1 Introduction

R [https://www.r-project.org/] adalah bahasa pemprograman statistik untuk analisis dan visualisasi yang cukup mudah untuk dipelajari. Dalam perkembangannya bahasa pemprograman R telah digunakan di segala sektor seperti perbankan, tech starup, universitas, NGO, dan masih banyak lagi. Saat ini R telah di kembangkan oleh peneliti dan praktisi menjadi tools untuk analisa Machine Learning untuk beragam kebutuhan seperti pemodelan lingkungan, prediksi saham, analisa ekologi, genetik, dan farmasi.

1.1.1 R Environment

Instalasi R cukup mudah dibandingkan bahasa pemograman lainnya, program dapat di unduh dengan mudah melalui Comprehensive R Archive Network (CRAN) [http:cran.r-project.org/]. Pada halaman CRAN terdapat pilihan download untuk Windows, Mac OS X, dan Linux.

Saat module ini dibuat, versi dari R terbaru adalah R version 3.5.3 (Great Truth). Versi terbaru secara otomatis akan mengikutkan fungsi-fungsi dari versi sebelumnya. Meskipun terkadang ada beberapa fungsi yang hilang saat versi R terbaru di keluarkan. Untuk itu sebelum melakukan instalasi atau update perlu membaca informasi detail perihal keterbaharuan dari versi R.

Untuk mempermudah proses penggunaan bahasa pemprograman R, disarankan untuk menjalankannya ada IDEs seperti RStudio [https://www.rstudio.com/].

Ada beberapa kemudahan-kemudahan yang diberikan oleh RStudio, yaitu:

  1. Bantuan code completion yang akan menampilkan daftar perintah saat kita telah mengetik karakter pada command prompt.
  2. Integrasi bantuan dalam satu layar, sehingga informasi hasil perintah help() dapat langsung ditampilkan.
  3. Tab environment yang menampilkan daftar objek yang dibuat.
  4. Tab display Files, Plots, Package yang juga terintegrasi pada satu layar.
  5. Penggunakan shortcut untuk posisi cursor masih dapat digunakan pada command prompt. Seperti shortcut tombol keyboard Ctrl dan panah kiri atau kanan untuk pindah cursor per kata.

Interface R Studio

Interface R Studio

[Selalu mengatur setting dari working directory dari menu -> Session -> Set Working Directory]

1.1.2 R Package

Alasan utama kenapa R banyak digunakan adalah ketersediaan package yang sangat beragam (per bulan Maret 2019 ada 15,000 additional packages) dan dapat membantu hampir semua analisa data (text, numeric, vector, raster). Package adalah

Instalasi Package dapat dilakukan dengan dua cara. Pertama melalui GUI yang tersedia di RStudio (Gambar 1).

[Tips: Ctrl+7]. Atau dengan menginputkan perintah pada Console:

install.packages("sp")

Instal Package - Tools - Instal Package

Instal Package - Tools - Instal Package

Package yang telah diinstall pada lingkungan R juga dikenal dengan istilah library. Untuk mengetahui lokasi penyimpanan library maka dapat digunakan perintah berikut ini.

.libPaths()

.libPaths("C:/Users/Manessa/Documents/R/win-library/3.5")

Untuk melihat daftar library yang telah terpasang pada lingkurangan R dapat menggunakan fungsi berikut ini.

library()

1.2 Basic of R

R adalah powerfull tools untuk perhitungan, data manipulasi, dan scientific computation. Untuk itu di perlukan pemahaman yang baik perihal operasi standar dan bentuk data. Seperti bahasa pemprograman pada umumnya R memiliki fungsi matematik, variable, fungsi, dan type data.

1.2.1 Basic Math

Sebagai bahasa pemprograman tentu saja R mampu melakukan operasi matematika seperti berikut:

9 + 9

## [1] 18

1 * 3 * 5

## [1] 15

9/3

## [1] 3

4/5

## [1] 0.8

2^3

## [1] 8

4 * 6 + 2

## [1] 26

(4 * 6) + 2

## [1] 26

4 * (6 + 2)

## [1] 32

1.2.2 Variable

Variabel adalah input yang disimpan untuk digunakan kembalo . Variable tidak hanya memuat angka atau karakter tetapi juga fungsi/formula, hasil analisis, dan plot.

Ada beberapa cara memerintahkan sebuah nilai sebagai variabel, tipe operasi yang valid adalah <- dan =.

[Tips: <- lebih disarankan].

x <- 6
x

## [1] 6

y <- 1
y

## [1] 1

z <- 3
z

## [1] 3

Penamaan variable dapat berupa gabungan huruf (sensitif terhadap capslock), angka, dan karakter seperti . atau **_**, akan tetapi angka tidak dapat disusun pada urutan pertama.

theVariable = 9
theVariable

## [1] 9

TheVariable

## Error in eval(expr, envir, enclos):  オブジェクト 'TheVariable' がありません

kesalahan penulisan adalah error yang paling sering terjadi, sehingga diperlukan ketelitian dalam penulisan.

[Tips: untuk pemula gunakan type huruf yang sama sepanjang script untuk meminimalisir error].

1.2.3 Calling Function

Fungsi sangat bermanfaat di R karena menjadikan code dapat diulangi dengan mudah. Fungsi merupakan sebuh perintah operasi di ditujukan kepada variable untuk membentuk variable baru atau memberikan informasi tentang variable.

print(theVariable)

## [1] 9

mean(y)

## [1] 1

length(y)

## [1] 1

assign("m", 2)
m

## [1] 2

1.2.4 Data type

Ada banyak tipe data di R yang dapat menyimpan beragam data/informasi, yakni sebagai berikut:

  1. Numeric adalah angka dapat berupa float(numeric) atau integer.

is.numeric(m)

## [1] TRUE

is.integer(m)

## [1] FALSE

class(m)

## [1] "numeric"

  1. Character adalah karakter huruf yang dapat berupa string atau factor

infimap <- "Spatial"
infimap

## [1] "Spatial"

infimaps <- as.factor(infimap)
infimaps

## [1] Spatial
## Levels: Spatial

nchar(infimap)

## [1] 7

class(infimap)

## [1] "character"

class(infimaps)

## [1] "factor"

nchar(infimaps)

  1. Dates (time-based) dapat berupa date atau POSIXct sebagai berikut:

waktu1 <- as.Date("2019-04-08")
waktu1

## [1] "2019-04-08"

class(waktu1)

## [1] "Date"

as.numeric(waktu1)

## [1] 17994

atau

waktu2 <- as.POSIXct("2019-04-08 10:30")
waktu2

## [1] "2019-04-08 10:30:00 +07"

class(waktu2)

## [1] "POSIXct" "POSIXt"

as.numeric(waktu2)

## [1] 1554694200

  1. Logical adalah cara untuk merepresentasikan data sebagai benar atau salah (True/False). Logical dapat di interpretasikan sebagai 1 apabila True dan 0 untuk False.

TRUE * 2

## [1] 2

FALSE * 2

## [1] 0

Learning = TRUE
class(Learning)

## [1] "logical"

is.logical(Learning)

## [1] TRUE

  1. Vector adalah kumpulan elemen dengan tipe yang sama dan tidak memiliki dimensi.

Fungsi c() digunakan untuk membuat data vector

contoh data vector yang berupa number

x <- c(2, 4, 6, 8, 10)
x

## [1]  2  4  6  8 10

contoh pengoprasian data vector yang berupa number

x + 1

## [1]  3  5  7  9 11

sqrt(x)

## [1] 1.414214 2.000000 2.449490 2.828427 3.162278

contoh data vector yang berupa character

y <- c("Manggis", "Mangga", "Duku", "Jeruk", "Rambutan", "Manggis", "Mangga", "Duku", "Jeruk", "Rambutan")
y

##  [1] "Manggis"  "Mangga"   "Duku"     "Jeruk"    "Rambutan" "Manggis"
##  [7] "Mangga"   "Duku"     "Jeruk"    "Rambutan"

Untuk mengakses suatu elemen dari suatu data vector nmengunakan fungsi berikut

#mengakses elemen ke-1 pada variable x
x[c(1)]

## [1] 2

#mengkakses elemen ke 2 pada variable y
y[c(2)]

## [1] "Mangga"

#mengkakses elemen ke 2 dan ke 5 dari variable y
y[c(2,5)]

## [1] "Mangga"   "Rambutan"

#jika ingin menganti suatu elemen pada sebuah data vector, misal character "Manggis" diganti dengan "Buah Naga"

y[c(1)] = "Buah Naga"
y

##  [1] "Buah Naga" "Mangga"    "Duku"      "Jeruk"     "Rambutan"
##  [6] "Manggis"   "Mangga"    "Duku"      "Jeruk"     "Rambutan"

set data vector ke factor jika ada kebutuhan untuk mengkategorikan nilai-nilai yang dimiliki oleh suatu variable maka factor dapat memberikan kemudahan

y.factor <- as.factor(y)
y.factor

##  [1] Buah Naga Mangga    Duku      Jeruk     Rambutan  Manggis   Mangga  
##  [8] Duku      Jeruk     Rambutan
## Levels: Buah Naga Duku Jeruk Mangga Manggis Rambutan

fungsi summary() dapat digunakan untuk mendapatkan detail dari data factor

summary(y.factor)

## Buah Naga      Duku     Jeruk    Mangga   Manggis  Rambutan
##         1         2         2         2         1         2

  1. Data frame adalah tabel yang memiliki kolom dan baris. Jumlah dari elemen seluruh kolom dalam data frame tetap harus sama

i <- 1:5
n <- c(0.1, 0.3, 0.5, -0.7, -0.9)
f <- c("Manggis", "Mangga", "Duku", "Jeruk", "Rambutan")
m <- data.frame(i,n,f)
m

##   i    n        f
## 1 1  0.1  Manggis
## 2 2  0.3   Mangga
## 3 3  0.5     Duku
## 4 4 -0.7    Jeruk
## 5 5 -0.9 Rambutan

Merubah header data frame

m <- data.frame(First = i, Second = n, Last = f)
m

##   First Second     Last
## 1     1    0.1  Manggis
## 2     2    0.3   Mangga
## 3     3    0.5     Duku
## 4     4   -0.7    Jeruk
## 5     5   -0.9 Rambutan

Memeriksa Jumlah Row pada Data Frame

nrow(m)

## [1] 5

memeriksa Jumlah kolom pada Data Frame

ncol(m)

## [1] 3

Memerikasa dimensi variabel Data Frame

dim(m)

## [1] 5 3

Memeriksa nama kolom

colnames(m)

## [1] "First"  "Second" "Last"

Memeriksa nama baris

rownames(m)

## [1] "1" "2" "3" "4" "5"

Merubah value pada nama baris

rownames(m) <- c("satu", "dua", "tiga", "empat", "lima")
rownames(m)

## [1] "satu"  "dua"   "tiga"  "empat" "lima"

memeriksa tipe data pembentuk objek

class(m)

## [1] "data.frame"

  1. List adalah kumpulan variable baik dengan format yang sama ataupun tidak. Sebagai contoh di dalam sebuat variable yang dibentuk dari fungsi list() dapat berisi atas tipe data skalar, vector, matrix dan lain-lain

Berikut adalah contoh dari variable yang dibentuk dari fungsi list().

infi <- list(m, "hi", c(1:4))
infi

## [[1]]
##       First Second     Last
## satu      1    0.1  Manggis
## dua       2    0.3   Mangga
## tiga      3    0.5     Duku
## empat     4   -0.7    Jeruk
## lima      5   -0.9 Rambutan
##
## [[2]]
## [1] "hi"
##
## [[3]]
## [1] 1 2 3 4

memeriksa tipe data pembentuk objek

class(infi)

## [1] "list"

Jika ingin mengakses nilai pertama pada variable mahasiswa maka digunakan kode seperti berikut.

#contoh fungsu mengakses elemen ke 2 dari fungsi list
infi[[2]]

## [1] "hi"

  1. Matrices adalah struktur data matematik yang mirip dengan data frame. Dalam membuat matrix harus diperhatikan bahwa:
  2. Seluruh kolom pada matrix harus berisi tipe data yang sama (numeric, char, boolean).
  3. Seluruh kolom harus memiliki panjang yang sama.

Berikut ini adalah contoh membuat matrix. Pada contoh di bawah ini akan dibuat dua buah matrix sederhana tetapi dengan menggunakan nilai parameter byrow yang berbeda. Secara default nilai parameter byrow adalah FALSE.

#byrow=FALSE
A <- matrix(1:10, nrow=5)
A

##      [,1] [,2]
## [1,]    1    6
## [2,]    2    7
## [3,]    3    8
## [4,]    4    9
## [5,]    5   10

#byrow=TRUE
B <- matrix(1:10, nrow=5, byrow=TRUE)
B

##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4
## [3,]    5    6
## [4,]    7    8
## [5,]    9   10

memeriksa dimensi variable

dim(A)

## [1] 5 2

merubah nama kolom

colnames(A) <- c("satu", "dua")
A

##      satu dua
## [1,]   1   6
## [2,]    2   7
## [3,]    3   8
## [4,]    4   9
## [5,]    5  10

memeriksa kelas data objek variable

class(A)

## [1] "matrix"

1.3 Logical function

1.3.1 if else

Statement if, merupakan fungsi yang sering digunakan dalam logical fungtion. Jika “kondisi” bernilai benar maka “fungsi” akan dijalankan, jika “kondisi” salah maka “fungsi” tidak akan dijalankan. Sebagai contoh:

if(5 > 4){print("5 lebih besar daripada 4")}

## [1] "5 lebih besar daripada 4"

Karena kondisi 5 > 4 adalah benar maka fungsi sprintf di dalam tanda { } dieksekusi.

Sintaks penggunaan statement ini yang lebih rumit dapat dilihat di bawah ini.

#if(kondisi 1) { ...} else if(kondisi 2) {...} else {...}

Untuk mencoba sintaks di atas maka dibuat fungsi sederhana dengan kode di bawah ini.

HitungNilai <- function(nilai){
if(nilai > 100 | nilai < 0) {sprintf("Nilai harus 0-100")} else {if(nilai >= 80){sprintf("A")} else if (nilai >= 70){sprintf("B")} else if (nilai >= 60){sprintf("C")} else if (nilai >= 50) {sprintf("D") } else {sprintf("E")}}}

HitungNilai(200)

## [1] "Nilai harus 0-100"

HitungNilai(-1)

## [1] "Nilai harus 0-100"

HitungNilai(40)

## [1] "E"

HitungNilai(70)

## [1] "B"

HitungNilai(90)

## [1] "A"

1.3.2 for loop

Statement for

Berikut ini adalah sintaks dari statement for.

# for(variable in sequence) {...}

Berikut adalah contoh penggunaan statement ini.

for(i in 1:5) {print(i)}

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5

data <- 1:10
for(i in data) {print(data[i])}

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10

Statement while

Cara lain untuk melakukan pengulangan adalah dengan menggunakan statement while. Sintaks dari statement ini adalah sebagai berikut.

i = 1;
while(i <= 5){
print(i);
i = i+1;
}

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5

Berikut ini adalah contoh penggunaan pengulangan while() untuk mengakses data pada vector

data_mahasiswa = c("ali", "mohammad", "reza", "faisal", "budi", "iwan", "wati","juna")
i=1
while(i <= 5){
print(data_mahasiswa[i])
i = i+1
}

## [1] "ali"
## [1] "mohammad"
## [1] "reza"
## [1] "faisal"
## [1] "budi"

1.4 Data Visualizating

1.4.1 Base Graphics

Berikut adalah grafik standar yang tidak begitu menarik untuk di gunakan.

data(airquality)
head(airquality)

##   Ozone Solar.R Wind Temp Month Day
## 1    41     190  7.4   67     5   1
## 2    36     118  8.0   72     5   2
## 3    12     149 12.6   74     5   3
## 4    18     313 11.5   62     5   4
## 5    NA      NA 14.3   56     5   5
## 6    28      NA 14.9   66     5   6

hist(airquality$Temp)

plot(Ozone ~ Wind, data = airquality)

boxplot(airquality)

1.4.2 ggplot2

R Package untuk membuat visualisasi data yang paling banyak di gunakan adalah ggplot2, karena memiliki tampilan yang menarik. Akan tetapi fungsi-fungsi di ggplot2 cukup rumit dalam pengaplikasiannya. Dibutuhkan banyak berlatih untuk dapat membuat grafik-grafik yang menarik dan mudah di pahami.

install.packages("ggplot2")

## package 'ggplot2' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
##  C:\Users\Manessa\AppData\Local\Temp\RtmpQbwzev\downloaded_packages

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.5.3

data(diamonds)
head(diamonds)

## # A tibble: 6 x 10
##   carat cut       color clarity depth table price     x     y     z
##   <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.23  Ideal     E     SI2      61.5    55   326  3.95  3.98  2.43
## 2 0.21  Premium   E     SI1      59.8    61   326  3.89  3.84  2.31
## 3 0.23  Good      E     VS1      56.9    65   327  4.05  4.07  2.31
## 4 0.290 Premium   I     VS2      62.4    58   334  4.2   4.23  2.63
## 5 0.31  Good      J     SI2      63.3    58   335  4.34  4.35  2.75
## 6 0.24  Very Good J     VVS2     62.8    57   336  3.94  3.96  2.48

g <- ggplot(data = diamonds, mapping = aes(x = carat, y = price))
g <- g + geom_point(aes(color = color))
g <- g + facet_wrap(~color)
g

g <- g + facet_grid(cut ~ clarity)
g

2.Spatial data in R

Prerequisites

Package yang digunakan pada sesi ini dapat di install dengan perintah berikut:

install.packages("devtools")
devtools::install_github("Nowosad/spDataLarge")
install.packages("sf")
install.packages("raster")
install.packages("spData")
install.packages("spDataLarge")

## Warning: package 'spDataLarge' is not available (for R version 3.5.1)

kemudian aktifasi package dengan perintah:

library(sf)          # classes and functions for vector data
library(raster)      # classes and functions for raster data
library(spData)        # load geographic data
library(spDataLarge)   # load larger geographic data

spData [https://github.com/Nowosad/spData] dan spDataLarge [https://github.com/Nowosad/spDataLarge] adalah database spasial online yang akan digunakan pada sesi ini.

2.1 Introduction

Sesi ini akan memberikan penjelasan singkat tentang data spasial vektor dan raster, serta operasi sederhana R untuk data spasial.

2.2 Vector data

Package sf dapat mewakili semua tipe geometri vektor yang umum (kelas data raster tidak didukung oleh sf): titik, garis, poligon, dan masing-masing versi ‘multi’ (Gambar 1) yang mengelompokkan fitur-fitur dari tipe yang sama menjadi satu fitur.

Gambar 1. Feature vektor

Gambar 1. Feature vektor

sf juga mendukung koleksi geometri, yang dapat berisi beberapa tipe geometri dalam satu objek. sf sebagian besar menggantikan ekosistem sp, yang terdiri dari sp , rgdal untuk data baca / tulis ) dan rgeos untuk operasi spasial. Isi dari package sf dapat dilihat dengan perintah operasi sebagai berikut:

vignette(package = "sf") # see which vignettes are available
vignette("sf1")          # an introduction to the package

objek fitur spasial dalam R disimpan dalam data frame, dengan informasi geografis menempati kolom khusus, biasanya dinamai geom atau 遯カ?eometri遯カ?. Kita akan menggunakan dataset dunia yang disediakan oleh spData, yang dimuat di awal bab ini (lihat [nowosad.github.io/spData] untuk daftar dataset yang dimuat oleh Package). file World adalah objek spasial yang berisi kolom atribut dan spasial, detail attribute data spasial dapat dilihat dengan operasi berikut dimana kolom terakhir berisi informasi geografis:

names(world)

Isi kolom geom pada file, memberikan informasi spasial object: world$geom adalah ‘kolom daftar’ yang berisi semua koordinat poligon negara. Paket sf menyediakan metode plot() untuk memvisualisasikan data geografis:

plot(world)

world$geom

summary(world)

Pemilihan informasi tertentu untuk di plot kan dapat dengan operasi perintah sebagai berikut:

plot(world["subregion"])

unique(world$region_un)

atau mengabungkan informasi attribut lain pada peta

plot(world["continent"], reset = FALSE)
cex = sqrt(world$pop) / 10000
world_cents = st_centroid(world, of_largest = TRUE)
plot(st_geometry(world_cents), add = TRUE, cex = cex)

2.3 Raster data

Data raster biasanya terdiri dari header raster dan matriks (dengan baris dan kolom) yang mewakili sel dengan interval yang sama (Gambar 2)

Gambar 2. Jenis data raster

Gambar 2. Jenis data raster

Peta raster biasanya mewakili fenomena berkelanjutan (continuous) seperti ketinggian, suhu, kepadatan populasi atau data spektral (Gambar 3A). Tentu saja, kita dapat merepresentasikan fitur-fitur yang terpisah (categorical) seperti kelas tanah atau tutupan lahan juga dengan bantuan model data raster (Gambar 3B).

Gambar 3. Examples of continuous and categorical rasters

Gambar 3. Examples of continuous and categorical rasters

Package yang digunakan untuk data raster adalah Raster. Gunakan operasi berikut untuk detail dari package

vignette("functions", package = "raster")

Untuk ilustrasi konsep raster, kami akan menggunakan dataset dari spDataLarge (perhatikan paket-paket ini dimuat pada awal bab). Ini terdiri dari beberapa objek raster dan satu objek vektor yang meliputi area Taman Nasional Zion (Utah, AS). Misalnya, srtm.tif adalah model elevasi digital area ini. Pertama, mari kita membuat objek RasterLayer bernama new_raster:

raster_filepath <- system.file("raster/srtm.tif", package = "spDataLarge")
new_raster <- raster(raster_filepath) #raster function only for single layer raster data

Mengetik nama raster ke konsol, akan mencetak header raster (luas, dimensi, resolusi, CRS) dan beberapa informasi tambahan (kelas, nama sumber data, ringkasan dari nilai-nilai raster):

new_raster

plot(new_raster)

Untuk data raster dengan multiple layer digunakan fungsi yang berbeda brick atau stack.

multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
r_brick = brick(multi_raster_file)
r_brick

2.4 Coordinate Reference System

Tipe data spasial vektor dan raster berbagi konsep yang intrinsik dengan data spasial. Mungkin yang paling mendasar dari ini adalah Sistem Referensi Koordinat (CRS), yang mendefinisikan bagaimana elemen spasial data berhubungan dengan permukaan Bumi. CRS bersifat geografis atau proyeksi.

Sistem koordinat geografis mengidentifikasi setiap lokasi di permukaan bumi menggunakan dua nilai - bujur dan lintang. Bujur adalah lokasi di arah Timur-Barat dalam jarak sudut dari bidang Meridian. Lintang adalah jarak sudut Utara atau Selatan bidang khatulistiwa. Jarak dalam CRS geografis diukur dalam meter. Daftar datum yang dapat digunakan dapat dilihat dengan perintah berikut

st_proj_info(type = "datum")

Proyeksi CRS didasarkan pada koordinat Cartesius pada permukaan datar (sistem kordinat proyeksi). Mereka memiliki asal, sumbu x dan y, dan unit pengukuran linier seperti meter. Sistem proyeksi yang tersedia dapat dilihat dengan

head(st_proj_info(type = "proj"), 6)

Dua cara utama untuk menggambarkan CRS dalam R adalah dengan kode epsg atau definisi proj4string.Kedua pendekatan ini memiliki kelebihan dan kekurangan. Kode epsg biasanya lebih pendek, dan karenanya lebih mudah diingat. Kode ini juga merujuk hanya satu sistem referensi koordinat yang terdefinisi dengan baik. Di sisi lain, definisi proj4string memungkinkan Anda lebih fleksibel ketika menentukan parameter yang berbeda seperti jenis proyeksi, datum dan ellipsoid. Dengan cara ini Anda dapat menentukan banyak proyeksi yang berbeda, dan memodifikasi yang sudah ada. Ini juga membuat pendekatan proj4string lebih rumit. epsg menunjuk ke tepat satu CRS tertentu.

vector_filepath = system.file("vector/zion.gpkg", package = "spDataLarge")
new_vector = st_read(vector_filepath)

st_crs(new_vector) # get CRS

new_vector = st_set_crs(new_vector, 4326)

3. Spatial Function

Prerequisites

Package yang perlu diinstall:

install.packages("rcartocolor")

## package 'rcartocolor' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
##  C:\Users\Manessa\AppData\Local\Temp\RtmpQbwzev\downloaded_packages

install.packages("grid")
install.packages("tmap")

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

install.packages("RQGIS")

## package 'RQGIS' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
##  C:\Users\Manessa\AppData\Local\Temp\RtmpQbwzev\downloaded_packages

install.packages("ows4R")

## package 'ows4R' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
##  C:\Users\Manessa\AppData\Local\Temp\RtmpQbwzev\downloaded_packages

install.packages("dplyr")

## package 'dplyr' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
##  C:\Users\Manessa\AppData\Local\Temp\RtmpQbwzev\downloaded_packages

Package yang perlu diaktifkan:

library(sf)
library(raster)
library(dplyr)
library(spData)
library(spDataLarge)
library(rcartocolor)
library(grid)
library(tmap)
library(RQGIS)
library(utils)

3.1 Spatial operation on vector data

3.1.1 Spatial subsetting

Spatial subsetting adalah proses memilih fitur objek spasial berdasarkan pada apakah object berkaitan atau tidak dengan objek lain. Contoh dari Spatial subsetting terhadap dataset nz dan nz_height dari Package spData. Data ini berisi wilayag Selandia Baru yang terdiri dari 16 wilayah utama dan 101 titik tertinggi.

plot(nz)
plot(nz_height)

canterbury = nz %>% filter(Name == "Canterbury")
canterbury_height = nz_height[canterbury, ]

3.1.2 Topological relations

Topological relations menggambarkan hubungan spasial antar objek. Untuk memahaminya, ada baiknya menggunakan contoh sederhana sebagai berikut. Gambar 4.2 berisi poligon (a), garis (l) dan beberapa titik (p), yang dibuat dalam kode di bawah ini.

# create a polygon
a_poly = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1))))
a = st_sfc(a_poly)
# create a line
l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2))
l = st_sfc(l_line)
# create points
p_matrix = matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2)
p_multi = st_multipoint(x = p_matrix)
p = st_cast(st_sfc(p_multi), "POINT")

query sederhana seperti: point p mana yang intersect dengan poligon a? Pertanyaan ini dapat dijawab menggunakan spatial predicate seperti do the objects intersect? Dapat di implementasikan pada fungsi sf sebagai berikut:

st_intersects(p, a)

Kebalikan dari st_intersects() adalah st_disjoint(), dimana hanya objek yang secara spasial tidak berkaitan akan terpilih:

st_disjoint(p, a, sparse = FALSE)[, 1]

st_within() returns TRUE hanya untuk objek yang berada didalam objek terpilih, sebagai contoh:

st_within(p, a, sparse = FALSE)[, 1]

Pada fungsi st_within() objek yang berada di batas tidak terindentifikasi, untuk itu diperlukan fungsi st_touches() untuk mengextrak informasi hingga area batas (border):

st_touches(p, a, sparse = FALSE)[, 1]

Kemudian dapat juga di peroleh informasi pada jarak tertentu dari area, mengunakan perintah sebagai berikut:

sel = st_is_within_distance(p, a, dist = 0.9) # can only return a sparse matrix
lengths(sel) > 0

3.1.3 Spatial joining

Mengabungkan dua data spasial lebih umum dikenal dengan istilah spatial overlay dengan fungsi st_join. Kode berikut adalah ilustrasinya:

set.seed(2018) # set seed for reproducibility
(bb_world = st_bbox(world)) # the world's bounds
random_df = tibble(
  x = runif(n = 10, min = bb_world[1], max = bb_world[3]),
  y = runif(n = 10, min = bb_world[2], max = bb_world[4])
)
random_points = random_df %>%
  st_as_sf(coords = c("x", "y")) %>% # set coordinates
  st_set_crs(4326) # set geographic CRS

world_random = world[random_points, ]

source("S:/Google/07_Training/05_RSpatial/Script/04-spatial-join.R")

## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

tmap_arrange(jm1, jm2, jm3, jm4, nrow = 2, ncol = 2)

3.1.4 Spatial data aggregation

Data aggregate menujukkan hasil statistik dari variable (umumnya rerata atau total) dan kaitannya dengan grouping variable. Dengan menggunakan sampel data New Zealand, kode berikut akan mereta ketinggian titik tinggi di tiap region.

nz_avheight = aggregate(x = nz_height, by = nz, FUN = mean)

Ditampilkan sebagai berikut:

library(tmap)
tm_shape(nz_avheight) +
  tm_fill("elevation", breaks = seq(27, 29, by = 0.5) * 1e2) +
  tm_borders()

3.1.5 Spatial distance

Jarak antara dua objek spasial di hitung menggunakan fungsi st_distance(). Kode berikut menunjukkan jarangk antara titik tertinggu dengan geographic centroid dari Canterbury region:

nz_heighest = nz_height %>% top_n(n = 1, wt = elevation)
canterbury_centroid = st_centroid(canterbury)

st_distance(nz_heighest, canterbury_centroid)

3.2 Spatial operation on raster data

3.2.1 Spatial subsetting

Spatial subsetting untuk data raster dapat menggunakan tiga macam fungsi [], mask(), dan overlay.

# create raster mask
rmask = elev
values(rmask) = sample(c(NA, TRUE), 36, replace = TRUE)

# spatial subsetting
maskElevOp = elev[rmask, drop = FALSE]            # with [ operator
maskElev = mask(elev, rmask)                      # with mask()
maskElevOvl = overlay(elev, rmask, fun = "max")   # with overlay

plot(maskElev)

3.2.2 Map algebra

Map algebra (or cartographic modeling) mengelompokkan operasi matematik raster menjadi 4 subckass sebagai berikut: 1. Local adalah per-cell operations. 2. Focal adalah neighborhood operations. (3 x 3 input cell block) 3. Zonal operations are similar to focal operations, but the surrounding pixel grid on which new values are computed can have irregular sizes and shapes. 4. Global or per-raster operations; that means the output cell derives its value potentially from one or several entire rasters.

3.2.3 Local operations

Clasifikasi data elev, dimana rentang 0–12, 12–24 and 24–36 di reclassified menjadi kelas 1, 2 and 3, sebagai berikut:

rcl = matrix(c(0, 12, 1, 12, 24, 2, 24, 36, 3), ncol = 3, byrow = TRUE)
recl = reclassify(elev, rcl = rcl)

Selain klasifikasi, raster algebra local operation dapat berupa penjumlahan, perkalian, atau operasi logika seperti berikut:

elev + elev
elev^2
log(elev)
elev > 5

3.2.4 Focal operations

Nama lain dari focal operation adalah kernel atau filter atau moving window yang dapat berupa 3 x 3 pixel menggunakan fungsi focal.

r_focal = focal(elev, w = matrix(1, nrow = 3, ncol = 3), fun = min)

3.2.5 Zonal operation

Zonal operations mirip dengan focal operations.Tetapi filter yang digunakan dapat berupa bentuk apa saja tidak khusus berupa jendela persegi. Fungsi yang digunakan adalah zonal dengan contoh aplikasi kode sebagai berikut:

z = zonal(elev, grain, fun = "mean") %>%
  as.data.frame()
z

3.2.6 Global operation

Global operations adalag kasus spesifik dari zonal operation dengan seluruh raster merepresentasikan satu zonal (satu informasi atau operasi), seperti perhitungan deskriptip statistik seperti minimum, maximum, atau menghitung jarang dan pembobotan raster.

min(elev)

## class       : RasterLayer
## dimensions  : 6, 6, 36  (nrow, ncol, ncell)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names       : layer
## values      : 1, 36  (min, max)

3.2.7 Aggregation and disaggregation

Dataset raster dapat di bedakan berdasrkan resolusinya (Spasial). Untuk mengatur resolusi spasial raster dapat digunakan fungsi aggregate() untuk mengurangi resolusi atau disaggregate() meningkatkan resolusi. Berikut contoh kode untuk memodifikasi resolusi spasial dari data dem

data("dem", package = "RQGIS")
dem_agg = aggregate(dem, fact = 5, fun = mean)

p_ar1 = tm_shape(dem) + tm_raster(style = "cont", legend.show = FALSE) +
  tm_layout(main.title = "Original", frame = FALSE)
p_ar2 = tm_shape(dem_agg) + tm_raster(style = "cont", legend.show = FALSE) +
  tm_layout(main.title = "Aggregated", frame = FALSE)

tmap_arrange(p_ar1, p_ar2, ncol = 2)

dem_disagg = disaggregate(dem_agg, fact = 5, method = "bilinear")
identical(dem, dem_disagg)

data(elev, package = "spData")
elev_agg = aggregate(elev, fact = 2, fun = mean)
plot(extend(elev, 1, 0), col = NA, legend = FALSE)
plot(elev_agg, add = TRUE)
xy_1 = xyFromCell(elev, 1)
xy_2 = xyFromCell(elev_agg, c(1, 2, 4, 5))
arrows(xy_2[, 1], xy_2[, 2], xy_1[, 1], xy_1[, 2], length = 0.1)
points(xy_1, pch = 16, col = "black")
points(xy_2, pch = 21, bg = "salmon", col = "black")
plot(rasterToPolygons(elev[1, drop = FALSE]), add = TRUE)

# spplot(elev_agg, col.regions = rev(terrain.colors(9)),
#        colorkey = FALSE, at = c(0, values(elev_agg + 1)),
#         sp.layout = list(
#           list("sp.polygons", rasterToPolygons(elev[1, drop = FALSE]),
#                first = FALSE),
#           list("sp.points", SpatialPoints(xy_1), pch = 21, first = FALSE),
#           list("sp.points", SpatialPoints(xy_2), pch = 21, bg = "salmon",
#                col = "black")
#           ))

3.5 Raster-vector interactions

3.5.1 Raster cropping

Pada analisa data spasial, sanggat umum beragam sumber digunakan, seperti data citra penginderan jauh (raster) dan batas administratif (vektor). Dimana sering kali data spasial lebih luas dari data admistrasi atau area yang akan dipetakan. untuk kasus seperti ini fungsi cropping and masking sangat berguna.

srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = st_read(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion = st_transform(zion, projection(srtm))

fungsi crop akan memotong citra sesuai area vektor

srtm_cropped = crop(srtm, zion)

Sedangkan fungsi mask akan menghilangkan pixel citra didalam area vektor.

srtm_masked = mask(srtm, zion)

srtm_inv_masked = mask(srtm, zion, inverse = T)

# TODO: split into reproducible script, e.g. in code/09-cropmask.R
library(tmap)
library(rcartocolor)
terrain_colors = carto_pal(7, "TealRose")
pz1 = tm_shape(srtm) + tm_raster(palette = terrain_colors, legend.show = FALSE, style = "cont") +
  tm_shape(zion) + tm_borders(lwd = 2) +
  tm_layout(main.title = "A. Original")

pz2 = tm_shape(srtm_cropped) + tm_raster(palette = terrain_colors, legend.show = FALSE, style = "cont") +
  tm_shape(zion) + tm_borders(lwd = 2) +
  tm_layout(main.title = "B. Crop")

pz3 = tm_shape(srtm_masked) + tm_raster(palette = terrain_colors, legend.show = FALSE, style = "cont") +
  tm_shape(zion) + tm_borders(lwd = 2) +
  tm_layout(main.title = "C. Mask")

pz4 = tm_shape(srtm_inv_masked) + tm_raster(palette = terrain_colors, legend.show = FALSE, style = "cont") +
  tm_shape(zion) + tm_borders(lwd = 2) +
  tm_layout(main.title = "D. Inverse mask")

tmap_arrange(pz1, pz2, pz3, pz4, ncol = 4)

3.5.2 Raster extraction

Raster extraction adalah proses untuk mengindentifikasi dan mengembalikan nilai pixel rasrer yang berasosiasi dengan feature vektor. Hasilnya tergantung pada feature vektor (points, lines or polygons)yang digunakan dan argument yang digunakan pada fungsi raster::extract().

data("zion_points", package = "spDataLarge")
zion_points$elevation = raster::extract(srtm, zion_points)

# Aim: demonstrate buffer arg in raster extract
elev_b1 = raster::extract(srtm, zion_points, buffer = 1000)

library(tmap)
library(grid)

rast_poly_point = tm_shape(srtm) +
  tm_raster(palette = terrain_colors, title = "Elevation (m)",
            legend.show = TRUE, style = "cont") +
  tm_shape(zion) +
  tm_borders(lwd = 2) +
  tm_shape(zion_points) +
  tm_symbols(shape = 1, col = "black") +
  tm_layout(legend.frame = TRUE, legend.position = c("right", "top"))

rast_poly_point

3.5.3 Rasterization

Rasterization adalah konversi objek vektor menjadi objek raster. Umumnya digunakan untuk pemodelan atau analisa kuantitatif (e.g., analysis of terrain).

install.packages("osmdata")

## package 'osmdata' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
##  C:\Users\Manessa\AppData\Local\Temp\RtmpQbwzev\downloaded_packages

library(osmdata)

q = add_osm_feature(opq = opq("London"), key = "network", value = "tfl_cycle_hire")
lnd_cycle_hire = osmdata_sf(q)
cycle_hire_osm = lnd_cycle_hire$osm_points
cycle_hire_osm_projected = st_transform(cycle_hire_osm, 27700)
raster_template = raster(extent(cycle_hire_osm_projected), resolution = 1000,
                         crs = st_crs(cycle_hire_osm_projected)$proj4string)

ch_raster1 = rasterize(cycle_hire_osm_projected, raster_template, field = 1)
plot(ch_raster1)

ch_raster2 = rasterize(cycle_hire_osm_projected, raster_template,
                       field = 1, fun = "count")
plot(ch_raster2)

3.5.4 Spatial vectorization

Spatial vectorization adalah kebalikan dari rasterization yang merubah objek raster menjadi vektor seperti points, lines or polygons.

data(dem, package = "RQGIS")
cl = rasterToContour(dem)
plot(dem, axes = FALSE)
plot(cl, add = TRUE)

3.6 Reproject Spatial Data

Ada dua tipe dari CRSs: geographic (窶詫on/lat?, with units in degrees longitude and latitude) dan projected (typically with units of meters from a datum). Ada beberapa atribut dari CRS, seperti projection, datum, and ellipsoid, Informasi detail mengenai CSR dapat dilihat pada [https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf]

Reproject umumnya digunakan bila ingin melakukan analisis yang memerlukan informasi jarak seperti sf_buffer atau sf_distance, apabila data spasial memiliki sistem proyeksi geographic maka analisis tersebut tidak dapat dilakukan.

Kode untuk proyeksi data vektor adalah sebagai berikut

jakarta_proj = data.frame(x = 106, y = -6.8) %>%
  st_as_sf(coords = 1:2, crs = 4326)

fungsi untuk mengecek informasi crs st_crs

st_crs(jakarta_proj)

fungsi transformasi dengan st_transform

jakartaUTM = st_transform(jakarta_proj, 32746)
st_crs(jakartaUTM)

Sedangkan code untuk transformasi proyeksi data raster dicontohka n menggunakan data Land cover (categorical maps) sebagian area di Utah, USA dari National Land Cover Database 2011 in the NAD83 / UTM zone 12N CRS. Prosesnya yang terdiri dari dua tahan membaca objek raster dengan fungsi raster, mengecek crs dengan fungsi crs, dan konfersi dengan fungsi projectRaster().

cat_raster = raster(system.file("raster/nlcd2011.tif", package = "spDataLarge"))
crs(cat_raster)

wgs84 = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
cat_raster_wgs84 = projectRaster(cat_raster, crs = wgs84)

3.7 Geographic data I/O

3.7.1 Retrieving open data

Saat ini ada beragam databased spatial online yang dapat diakses secara umum, seperti

Sedangkan di Indonesia databased yang menyediakan basedata adalaha Inageoportal, tetapi di perlukan pembuatan akun untuk mengakses data. Sehingga kode yang diperlukan menjadi tidak begitu simpel (akan di coba pada sesi studi kasus)

Sebagai contoh data yang dapat dengan mudah diakses tanpa login seperti data US National Parks data from: catalog.data.gov/dataset/national-parks. Kode untuk mengakses data tersebut adalah sebagai berikut :

download.file(url = "http://nrdata.nps.gov/programs/lands/nps_boundary.zip",
              destfile = "nps_boundary.zip")
unzip(zipfile = "nps_boundary.zip")
usa_parks = st_read(dsn = "nps_boundary.shp")

summary(usa_parks)

3.7.2 Geographic data packages

R packages yang di kembangkang untuk mengakses data geografik di tampilkan pada table. Package ini menyediakan interface untuk library spasial atau geoportals sehingga akses data menjadi lebih mudah.

datapackages = tibble::tribble(
  ~`Package`, ~Description,
  "getlandsat", "Provides access to Landsat 8 data.",
  "osmdata", "Download and import of OpenStreetMap data.",
  "raster", "getData() imports administrative, elevation, WorldClim data.",
  "rnaturalearth", "Access to Natural Earth vector and raster data.",
  "rnoaa", "Imports National Oceanic and Atmospheric Administration (NOAA) climate data.",
  "rWBclimate", "Access World Bank climate data."
)

Batas negara terkadang menjadi hal yang penting, yang dapat di akses dari fungsi ne_countries() dari package rnaturalearth sebagai berikut:

library(rnaturalearth)
Ina = ne_countries(scale = 50, country = "Indonesia")
class(Ina)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

# alternative way of accessing the data, with raster::getData()
#raster::getData("GADM", country = "idn", level = 0)

Setting standar dari rnaturalearth memyediakan objek spasial dalam bentuk class Spatial. Objek ini dapat dikonversi menjadi objeksf dengan fungsi st_as_sf() sebagai berikut:

Ina_sf = st_as_sf(Ina)

Ina_unit <- ne_download(scale = 10, type = 'map_units', category = 'cultural')

## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\Manessa\AppData\Local\Temp\RtmpQbwzev", layer: "ne_10m_admin_0_map_units"
## with 295 features
## It has 94 fields
## Integer64 fields read as strings:  NE_ID

airport <- ne_download(scale = 10, type = 'airports', category = 'cultural')

## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\Manessa\AppData\Local\Temp\RtmpQbwzev", layer: "ne_10m_airports"
## with 891 features
## It has 35 fields
## Integer64 fields read as strings:  ne_id

plot(airport)

Unit_ina = crop(Ina_unit , extent(Ina))
airport_Ina = crop(airport, extent(Ina))

Contoh selanjutnya adalah series dari data raster yang berisi global monthly precipitation dengan resolusi sepuluh menit.The result is a multilayer object of class RasterStack.

#library(raster)
worldclim_prec = getData(name = "worldclim", var = "prec", res = 10)
#class(worldclim_prec)
plot(worldclim_prec)

3.7.3 Geographic web services

untuk mengakse web API dengan data spasial, In an effort to standardize web APIs for accessing spatial data, Open Geospatial Consortium (OGC) telah menentukan beberapa spesifikasi dengan nama OWS (OGC Web Services). Ketentuan ini meliputi Web Feature Service (WFS), Web Map Service (WMS), Web Map Tile Service (WMTS), Web Coverage Service (WCS) dan Wep Processing Service (WPS). Map servers seperti PostGIS sudah mengadopsi protocol ini, OWS API memiliki query yang sama dengan web API lainnya. Berikut adalah contoh geoserver yang di sediakan FAO:

library(ows4R)
wfs = WFSClient$new("http://www.fao.org/figis/geoserver/wfs",
                      serviceVersion = "1.0.0", logger = "INFO")

## [ows4R][INFO] OWSGetCapabilities - Fetching http://www.fao.org/figis/geoserver/wfs?service=WFS&version=1.0.0&request=GetCapabilities

fao_areas = wfs$getFeatures("area:FAO_AREAS")

## [ows4R][INFO] WFSClient - Fetching features for 'area:FAO_AREAS' ...
## [ows4R][INFO] WFSGetFeature - Fetching http://www.fao.org/figis/geoserver/wfs?service=WFS&version=1.0.0&typeName=area:FAO_AREAS&logger=INFO&request=GetFeature
## [ows4R][INFO] WFSDescribeFeatureType - Fetching http://www.fao.org/figis/geoserver/wfs?service=WFS&version=1.0.0&typeName=area:FAO_AREAS&request=DescribeFeatureType

# not shown as too verbose an example already
area_27 = wfs$getFeatures("area:FAO_AREAS",
                          cql_filter = URLencode("F_CODE= '27'"))

## [ows4R][INFO] WFSClient - Fetching features for 'area:FAO_AREAS' ...
## [ows4R][INFO] WFSGetFeature - Fetching http://www.fao.org/figis/geoserver/wfs?service=WFS&version=1.0.0&typeName=area:FAO_AREAS&logger=INFO&cql_filter=F_CODE=%20'27'&request=GetFeature

plot(area_27)

3.7.4 Data input (I)

  1. Untuk vector data berupa data tabular dengan posisi koordinat dapat menggunakan perintah sebagai berikut:

#ibrary(sf)

world_txt = system.file("misc/world_wkt.csv", package = "spData")
world_wkt = read_sf(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT")
# the same as
world_wkt = st_read(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT",
                    quiet = TRUE, stringsAsFactors = FALSE)

atau

store_csv = read.csv("S:/Google/07_Training/05_RSpatial/Script/store.csv", sep=";")
summary(store_csv)

##        X                                               X1   
##  Min.   : 1.00   Amx                                    : 1 
##  1st Qu.:15.75   Anticimex Pest Management Pte Ltd      : 1 
##  Median :30.50   Awox Pte Ltd                           : 1 
##  Mean   :30.50   B.i.G. Bioengineered in Germany Pte Ltd: 1 
##  3rd Qu.:45.25   BorgWarner Beru Systems GmbH           : 1 
##  Max.   :60.00   C-Mar Asia Pte. Ltd.                   :
##                  (Other)                                :54 
##        X2              X3      
##  Min.   :1.321   Min.   :103.7 
##  1st Qu.:1.325   1st Qu.:103.7 
##  Median :1.326   Median :103.7 
##  Mean   :1.326   Mean   :103.7 
##  3rd Qu.:1.327   3rd Qu.:103.7 
##  Max.   :1.330   Max.   :103.8 
##                                
##                                                                  X4   
##  3 International Business Park Road                               : 5 
##  25 International Business Park Road, Singapore                   : 4 
##  3 International Business Park Road, Singapore                    : 4 
##  25 International Business Park Road                              : 3 
##  29 International Business Park Road                              : 3 
##  7 International Business Park Road, Techquest Building, Singapore: 3 
##  (Other)                                                          :38

coordinates(store_csv) <- c("X3", "X2")
store_csv

## class       : SpatialPointsDataFrame
## features    : 60
## extent      : 103.7395, 103.754, 1.32051, 1.330237  (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## variables   : 3
## names       :  X,                    X1,                                                                             X4
## min values  :  1,                   Amx, #03-28, Nordic European Centre,, 3 International Business Park Road, Singapore
## max values  : 60, YKK AP FACADE PTE LTD,                      The Strategy, Tower 2, 2 International Business Park Road

untuk data kml

u = "https://developers.google.com/kml/documentation/KML_Samples.kml"
download.file(u, "KML_Samples.kml")
st_layers("KML_Samples.kml")

## Driver: KML
## Available layers:
##              layer_name  geometry_type features fields
## 1            Placemarks       3D Point        3      2
## 2      Highlighted Icon       3D Point        1      2
## 3                 Paths 3D Line String        6      2
## 4         Google Campus     3D Polygon        4      2
## 5      Extruded Polygon     3D Polygon        1      2
## 6 Absolute and Relative     3D Polygon        4      2

kml = read_sf("KML_Samples.kml", layer = "Placemarks")

untuk data shapefile dapat menggunakan kode berikut

shpjakarta <- shapefile("S:/Google/07_Training/05_RSpatial/Input/Jakarta.shp")

plot(shpjakarta)

untuk membuka data atribut table buka shapefile sebagai dataframe

atributshp <- as.data.frame(shpjakarta)

atributshp

##   PROVNO KABKOTNO KODE2010    PROVINSI           KABKOT
## 0     31       01     3101 DKI JAKARTA KEPULAUAN SERIBU
## 1     31       71     3171 DKI JAKARTA  JAKARTA SELATAN
## 2     31       72     3172 DKI JAKARTA    JAKARTA TIMUR
## 3     31       73     3173 DKI JAKARTA    JAKARTA PUSAT
## 4     31       74     3174 DKI JAKARTA    JAKARTA BARAT
## 5     31       75     3175 DKI JAKARTA    JAKARTA UTARA

  1. Raster data

Sama dengan data vektor, data raster dapat berupa beragam format yang terkadang terdiri dari multilayer.Fungsi yang digunakan untuk membaca file raster adalah raster().

raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
single_layer = raster(raster_filepath)

Apabila data yang ingin di baca adalah In case you want to read in a single band from a multilayer file, use the band parameter to indicate a specific layer.

multilayer_filepath = system.file("raster/landsat.tif", package = "spDataLarge")
band3 = raster(multilayer_filepath, band = 3)

Untuk menggunakan semua band gunakan fungsi brick() atau stack().

multilayer_brick = brick(multilayer_filepath)
multilayer_stack = stack(multilayer_filepath)

3.7.5 Data output (0)

  1. Shapfile to excel file

Packages yang digunakan untuk membaca data excel adalah “openxlsx”

install.packages("openxlsx")

## package 'openxlsx' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
##  C:\Users\Manessa\AppData\Local\Temp\RtmpQbwzev\downloaded_packages

konversi shapfile menjadi tabel excel adalah sebagai berikut

library(openxlsx)

## Warning: package 'openxlsx' was built under R version 3.5.3

library(shapefiles)

## Warning: package 'shapefiles' was built under R version 3.5.2

## Loading required package: foreign

##
## Attaching package: 'shapefiles'

## The following objects are masked from 'package:foreign':
##
##     read.dbf, write.dbf

shpjakarta <- shapefile("S:/Google/07_Training/05_RSpatial/Input/Jakarta.shp")
atributshp <- as.data.frame(shpjakarta)

data_shp_to_xlsx <- write.xlsx(atributshp, "S:/Google/07_Training/05_RSpatial/Output/data_shp_to_xlsx.xlsx")

## Note: zip::zip() is deprecated, please use zip::zipr() instead

  1. Shapfile to CSV

shpjakarta <- shapefile("S:/Google/07_Training/05_RSpatial//Input/Jakarta.shp")
atributshp <- as.data.frame(shpjakarta)

data_shp_to_csv <- write.csv(atributshp, "S:/Google/07_Training/05_RSpatial/data_shp_to_csv.csv")

4. Making maps in R

Prerequisites

Package yang perlu diinstall:

listPackage <- c("sf", "raster", "dplyr", "spData", "spDataLarge",
                 "tmap", "leaflet", "mapview", "ggplot2", "shiny",
                 "cartogram")

for (i in listPackage){
  print(i)
  if(!requireNamespace(paste(i), quietly = TRUE))
    install.packages(paste(i), quiet = TRUE)
}

## [1] "sf"
## [1] "raster"
## [1] "dplyr"
## [1] "spData"
## [1] "spDataLarge"
## [1] "tmap"
## [1] "leaflet"
## [1] "mapview"
## [1] "ggplot2"
## [1] "shiny"
## [1] "cartogram"

Package yang perlu diaktifasi:

library(tmap)   # for static and interactive maps
library(leaflet) # for interactive maps

## Warning: package 'leaflet' was built under R version 3.5.3

library(mapview) # for interactive maps

## Warning: package 'mapview' was built under R version 3.5.3

library(ggplot2) # tidyverse vis package
library(shiny)   # for web applications

4.1 Introduction

Salah satu hal yang penting dari penggunaan data spasial adalah visualisasi dari setiap proses analisa data. Dimana aspek kartografi menjadi esensial untuk dapat mengkomunikasikan informasi dari data spasial. R menyediakan beragam package untuk mendukung transfer kreatifitas dalam proses visualisasi data spasial. Training merekondasikan penggunaan package tmap dan ggplot2 karena fungsi yang di miliki cukup untuk dapat membuat peta yang informatif dan menarik.

4.2 Static maps

Static maps adalah peta yang paling umum di buat. Dimana peta hasil akan di cetak dengan format Gambar meliputi .jpg, .png and .pdf. Meskipun peta interaktif adalah hal yang menarik dan inovatif tetapi untuk beberapa kebutuhan peta static tetap menjadi pelihan terbaik. Di awal training fungsi standar plot() digunakan untuk menampilkan data raster dan vektor menjadi peta statik, fungsi ini sangat simple tetapi hal yang didapatkan kurang elegan dan menarik. Package tmap merupakan package yang fleksible dan powerful dengan hasil yang baik. Kode yang digunakan cukup simple tetapi tetap memerlukan latihan yang rutin untuk mendapatkan hasil yang baik. Hal yang menarik dan cukup berbeda dengan package kartografi lainnya adalah kemampuan untuk membuat peta statik dan interaktif hanya dengan menyetel pengaturan via tmap_mode().

4.2.1 tmap basics

Seperti fungsi di package ggplot2, tmap berbasis ‘grammar of graphics’. Ini meliputi pemilahan antara data input dan layer astetiknya meliputi geometry, warna, dan variable visual lainnya. Basic fungsi nya adalah tm_shape() (mendefinisikan input data, raster dan vector objects), diikuti oleh satu atau lebih elemen layer seperti tm_fill() and tm_dots(). Kode berikut menunjukkan proses layering:

# Add fill layer to nz shape
tm_shape(Ina) +
  tm_fill()
# Add border layer to nz shape
tm_shape(Ina) +
  tm_borders()
# Add fill and border layers to nz shape
tm_shape(Ina) +
  tm_fill() +
  tm_borders()

4.2.2 Map objects

Keungulan dari tmap adalah kemampuannya untuk menyimpan object yang merepresentasikan peta. Kode berikut menunjukkan peta dapat disimpan sebagai object variable yang dapat digunakan pada kode berikutnya :

map_Ina = tm_shape(Ina) + tm_polygons() # tm_polygon() adalah tm_fill() + tm_borders()
class(map_Ina)

## [1] "tmap"

map_Ina

map_Ina dapat di tampilkan lagi atau di modifikasi dengan menambahkan layer lain seperti cerikut:

map_Ina = map_Ina +
  tm_shape(airport_Ina) + tm_bubbles(size = 0.5) #coba tm_dot()

map_Ina

4.2.3 Aesthetics

Peta statik yang telah di buat menunjukkan standar setting visualisasi data. Dimana gray shade dan black line digunakan untuk semua fungsi. Tentu saja setting standar ini ini dapat di modifikasi sesuai dengan kebutuhan. fungsi umum yang digunakan untuk meningkatkan estetika tampilan adalah warna, transparency, ketebalan dan tipe garis yang diatur dengan argumen col, alpha, lwd, dan lty.

ma1 = tm_shape(Ina) + tm_fill(col = "red")
ma2 = tm_shape(Ina) + tm_fill(col = "red", alpha = 0.3)
ma3 = tm_shape(Ina) + tm_borders(col = "blue")
ma4 = tm_shape(Ina) + tm_borders(lwd = 3)
ma5 = tm_shape(Ina) + tm_borders(lty = 2)
ma6 = tm_shape(Ina) + tm_fill(col = "red", alpha = 0.3) +
  tm_borders(col = "blue", lwd = 3, lty = 2)
tmap_arrange(ma1, ma2, ma3, ma4, ma5, ma6)

4.2.4 Color settings

Pengaturan warna dilakukan dengan ketentuan sebagai berikut - Setting warna standar menggunakan ‘pretty’ breaks. - breaks memungkinkan untuk mengatur rentang nilai untuk tiap warna - n adalah kumpulan nilai yang di konvesi dalam jumlah kategori/kelas - palette adalah kumpulan/gradasi warna seperti BuGn.

ma1 + tm_shape(airport_Ina) + tm_bubbles(col = "scalerank")

breaks = c(0, 1, 3, 5)
ma1 + tm_shape(airport_Ina) + tm_bubbles(col = "scalerank", breaks = breaks)

ma1 + tm_shape(airport_Ina) + tm_bubbles(col = "scalerank", n = 10)

ma1 + tm_shape(airport_Ina) + tm_bubbles(col = "scalerank", palette = "BuGn")

MapIna2a = tm_shape(Unit_ina) + tm_polygons("SOVEREIGNT", palette = "Blues")

MapIna2b = tm_shape(Unit_ina) + tm_polygons("SOVEREIGNT", palette = "YlOrBr")

tmap_arrange(MapIna2a, MapIna2b)

4.2.5 Layouts

Layout peta merujuk pada pengabungan combinasi dari elemen peta menajdi peta utuh secara kartografis. Elemen peta meliputi judul, skala, direction, dan inset. Fungsi yang digunakan adalah seperti tm_compass() and tm_scale_bar()..

MapIna2a +
  tm_compass(type = "8star", position = c("left", "bottom")) +
  tm_scale_bar(breaks = c(0, 100, 200), size = 1)

map_Ina + tm_layout(title = "Indonesia")

 

map_Ina + tm_layout(scale = 5)

 

 

 

 

 

map_Ina + tm_layout(bg.color = "lightblue")

 

map_Ina + tm_layout(frame = FALSE)

fungsi lain dari (tm_layout)[https://www.rdocumentation.org/packages/tmap/versions/2.2/topics/tm_layout]. Sedangkan template yang dapat digunankan dari tmap adalah

MapIna2 + tm_style("bw")

 

MapIna2 + tm_style("classic")

#

MapIna2 + tm_style("cobalt")

#

MapIna2 + tm_style("col_blind")

#

4.2.6 Faceted maps

Faceted maps atay multiple theme map adalah beberapa peta yang disusun secara beraturan dengan aturan tertentu dalam penyusunannya. Contoh kode nya sebagai berikut:

urb_1970_2030 = urban_agglomerations %>%
  filter(year %in% c(1970, 1990, 2010, 2030))
tm_shape(world) + tm_polygons() +
  tm_shape(urb_1970_2030) + tm_symbols(col = "black", border.col = "white",
                                       size = "population_millions") +
  tm_facets(by = "year", nrow = 2, free.coords = FALSE)

4.3 Interactive maps

Peta interaktif dapat berupa banyak bentuk, dan mampu memberikan informasi spasial menjadi lebih menarik terutama yang menginginkan interaksi langsung dengan data. Hal yang paling penting dari peta interaktif adalah dapat memuatnya dalam applikasi web map. Dengan mengkombinasikan package seperti tmap, leaflet, dan shiny hal ini dapat dilakukan.

4.3.1 tmap

Untuk dapat menampilkan peta interaktif dengan package tmap, perlu memindakan setting dengan perintah tmap_mode("view").

tmap_mode("view")
map_Ina

kombinasi fungsi dari package leaflet dan tmap

tmap_mode("view")
m_tmview = map_Ina
tmap_leaflet(m_tmview)

Setelah mode interaktif aktif, maka perlu milih basemap yang digunakan dengan fungsi tm_basemap() (atau tmap_options()) sebagai berikut:

map_Ina + tm_basemap(server = "OpenTopoMap")

fungsi facet juga dapat digunakan dalam tampilan interaktif sebagai berikut

world_coffee = left_join(world, coffee_data, by = "name_long")
facets = c("coffee_production_2016", "coffee_production_2017")
tm_shape(world_coffee) + tm_polygons(facets) +
  tm_facets(nrow = 1, sync = TRUE)

Switch tmap menjadi plot untuk menghasilkan peta statik:

tmap_mode("plot")

4.3.2 leaflet

leaflet adalah interactive mapping package di R yang paling sering digunakan. Package ini menyediakan beragam Leaflet JavaScript library dan argumen yang memudahkan proses pengembangan leafletjs.com).

Peta leaflet di buat menggunakan fungsi leaflet(), hasilnya adalah leaflet map object yang dapat digunakan untuk leaflet functions lain. Hal ini memungkinkan multiple map layer di tampilkan secara interaktif seperti berikut rstudio.github.io/leaflet/:

pal = colorNumeric("RdYlBu", domain = cycle_hire$nbikes)
leaflet(data = cycle_hire) %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  addCircles(col = ~pal(nbikes), opacity = 0.9) %>%
  addPolygons(data = lnd, fill = FALSE) %>%
  addLegend(pal = pal, values = ~nbikes) %>%
  setView(lng = -0.1, 51.5, zoom = 12) %>%
  addMiniMap()

4.4 Mapping applications

Untuk aplikasi web, shiny package menyediakan opsi yang relatif mudah dalam penyusunan aplikasi. Dengan aplikasi ini ini peta interaktif yang dibangun dengan package leaflet dapat digabungkan dengan fungsi lain, fungsi renderLeaflet() adalah penyatu antara leaflet dan shiny, dokumentasi detail dapat di akses di Shiny integration. Sedangkan untuk shiny dapat dilihat secara detail pada shiny.rstudio.com.

aplikasi shiny dibangun dengan dua elemen yang merefleksikan web application development: ‘front end’ (the bit the user sees) and ‘back end’ code. Elemen ini dinamakan sebagai ui and server dengan nama R script app.R, di simpan dalam ‘app folder’. Sebagai contoh Error! Hyperlink reference not valid.

Contoh yang lain adalah sebagai berikut

library(shiny)    # for shiny apps
library(leaflet)  # renderLeaflet function
library(spData)   # loads the world dataset
ui = fluidPage(
  sliderInput(inputId = "life", "Life expectancy", 49, 84, value = 80),
      leafletOutput(outputId = "map")
  )
server = function(input, output) {
  output$map = renderLeaflet({
    leaflet() %>% addProviderTiles("OpenStreetMap.BlackAndWhite") %>%
      addPolygons(data = world[world$lifeExp < input$life, ])})
}
shinyApp(ui, server)

4.5 Other mapping packages

4.5.1 ggplot2

Package ggplot2 tidak hanya dapat di pergunakan untuk membuat grafik data non spasial. Package banyak digunakan karena ketersedian add-on function untuk menyelesaikan beragam permasalahan seperti pembuatan peta. Penulisan code dengan ggplot2 package sama dengan tmap berbasis grammar of graphic*. Kode awal dari gglot2 adalah ggplot() yang diikuti dengan pemangilan layer-layer berikutnya, yakni + geom_*(), diman * adalah bentuk dari data type data geom_sf() (for sf objects), geom_points() (for points), geom_polygon(); labs() untuk setting penamaan judul, axis.

library(ggplot2)
g1 = ggplot() + geom_sf(data = nz, aes(fill = Median_income)) +
  geom_sf(data = nz_height) +
  scale_x_continuous(breaks = c(170, 175))
g1

4.5.2 cartogram

cartogram package menyediakan fungsi untuk membuat modifikasi geometry untuk merepresentasikan variable pemetaan.

library(cartogram)
us_states2163 = st_transform(us_states, 2163)
us_states2163_ncont = cartogram_ncont(us_states2163, "total_pop_15")
us_states2163_dorling = cartogram_dorling(us_states2163, "total_pop_15")

carto_map3 = tm_shape(us_states2163_ncont) +
  tm_polygons("total_pop_15", title = "Total population", palette = "BuPu") +
  tm_layout(legend.show = FALSE)
carto_map4 = tm_shape(us_states2163_dorling) +
  tm_polygons("total_pop_15", title = "Total population", palette = "BuPu") +
  tm_layout(legend.show = FALSE)
carto_map_34legend = tm_shape(us_states2163_dorling) +
  tm_polygons("total_pop_15", title = "Total population", palette = "BuPu") +
  tm_layout(legend.only = TRUE)
tmap_arrange(carto_map3, carto_map4, carto_map_34legend, ncol = 3)

Machine Learning in R

Materi mengenai Machine Learning in R ada di link berikut

next session