Spasial dengan R
Tugas STA500 Metode Penelitian Kuantitatif Materi Spasial
Pendahuluan
Jenis Data Spasial
Data spasial terdiri dari data vector dan data raster. Contoh data vektor adalah data peta administrasi. Data vektor terdiri dari point, line, polygon. Contoh data raster adalah data pada google eart. Data raster terdiri dari Rasterlayer dan Rasterstack.
Packages
Packages yang digunakan adalah sp, rgdal, spData, readxl, raster, spdep
library(sp)
library(rgdal) # membaca shapefile
library(spData) #data columbus
library(readxl) # membaca file excel
library(raster) # untuk menambahkan text pada peta
library(spdep) # pembobot data spasialData
Data dalam package R
Dengan memanfaatkan datasets yang ada di R pada package yang sesuai tema. Sebagai contoh untuk data spasial ini dapat menggunakan data yang berada pada package spData.
ls() # melihat data set yang sedang dijalankan pada R kota## character(0)data(columbus) # menggunakan data columbus yang berada di package spData
ls() # panggil list lagi untuk melihat berhasil atau tiDak## [1] "bbs" "col.gal.nb" "columbus" "coords" "polys"Untuk melihat struktur data dapat menggunakan fungsi str diikuti oleh nama datasets nya.
str(columbus)## 'data.frame': 49 obs. of 22 variables:
## $ AREA : num 0.3094 0.2593 0.1925 0.0838 0.4889 ...
## $ PERIMETER : num 2.44 2.24 2.19 1.43 3 ...
## $ COLUMBUS. : int 2 3 4 5 6 7 8 9 10 11 ...
## $ COLUMBUS.I: int 5 1 6 2 7 8 4 3 18 10 ...
## $ POLYID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ NEIG : int 5 1 6 2 7 8 4 3 18 10 ...
## $ HOVAL : num 80.5 44.6 26.4 33.2 23.2 ...
## $ INC : num 19.53 21.23 15.96 4.48 11.25 ...
## $ CRIME : num 15.7 18.8 30.6 32.4 50.7 ...
## $ OPEN : num 2.851 5.297 4.535 0.394 0.406 ...
## $ PLUMB : num 0.217 0.321 0.374 1.187 0.625 ...
## $ DISCBD : num 5.03 4.27 3.89 3.7 2.83 3.78 2.74 2.89 3.17 4.33 ...
## $ X : num 38.8 35.6 39.8 36.5 40 ...
## $ Y : num 44.1 42.4 41.2 40.5 38 ...
## $ AREA : num 10.39 8.62 6.98 2.91 16.83 ...
## $ NSA : num 1 1 1 1 1 1 1 1 1 1 ...
## $ NSB : num 1 1 1 1 1 1 1 1 1 1 ...
## $ EW : num 1 0 1 0 1 1 0 0 1 1 ...
## $ CP : num 0 0 0 0 0 0 0 0 0 0 ...
## $ THOUS : num 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 ...
## $ NEIGNO : num 1005 1001 1006 1002 1007 ...
## $ PERIM : num 2.44 2.24 2.19 1.43 3 ...Dapat terlihat bahwa data tersebut berbentuk data.frame terdiri dari 49 observasi dengan 22 peubah. Masing- masing variabel pun diberikan penjelasan jenis/tipe datanya. Untuk angka yang muncul adalah cuplikan data dari masing-masing peubah.
Untuk melihat deskripsi dari data dapat menggunakan fungsi help
help(columbus)Columbus OH spatial analysis data set. The columbus data frame has 49 rows and 22 columns. Unit of analysis: 49 neighbourhoods in Columbus, OH, 1980 data. In addition the data set includes a polylist object polys with the boundaries of the neighbourhoods, a matrix of polygon centroids coords, and col.gal.nb, the neighbours list from an original GAL-format file. The matrix bbs is DEPRECATED, but retained for other packages using this data set.
Jika akan dilakukan analisis dapat dilihat dari deksripsi setiap peubah ini.
External Data Input
Selain menggunakan data yang berada di dalam package R, kita juga dapat memasukkan data dari luar package. Sebagai contoh, akan digunakan data gizi buruk berupa file excel.
setwd("D:\\Kuliah S2 IPB\\Bahan Kuliah\\Semester 1 SSD 2020\\01 Metode Penelitian Kuantitatif\\uts\\spasial\\SPASIAL DENGAN R\\Spatial with R")
# untuk melakukan pengaturan working directory
data.giziburuk <- read_excel("data giziburuk jabar.xlsx", sheet = "Sheet1")
# memasukkan file excel ke dalam objek data.giziburukUntuk melihat dimensi dari data dapat menggunakan fungsi dim. Terlihat bahwa dimensi dari data tersebut adalah 27 observasi dan 5 variabel/peubah.
dim(data.giziburuk)## [1] 27 5Untuk melihat struktur data dapat menggunakan fungsi str. Terlihat bahwa data gizi buruk ini berbentuk tibble terdiri dari 27 observasi dengan 5 peubah. Sebagai contoh peubah kode kabko terdiri dari kode kabupaten kota berupa tipe data numerik.
str(data.giziburuk)## tibble [27 x 5] (S3: tbl_df/tbl/data.frame)
## $ kode kabko: num [1:27] 3201 3202 3203 3204 3205 ...
## $ KABKOT : chr [1:27] "BOGOR" "SUKABUMI" "CIANJUR" "BANDUNG" ...
## $ Giziburuk : num [1:27] 211 117 133 186 132 110 104 111 159 129 ...
## $ Long : num [1:27] 107 107 107 108 108 ...
## $ Lat : num [1:27] -6.54 -7.07 -7.05 -7.07 -7.34 ...Untuk melihat data gizi buruk tersebut dalam bentuk tabel dapat menggunakan fungsi view. Akan memunculkan sebuah menu berisi tabel tersebut. Bentuknya sama seperti excel.
View(data.giziburuk)Plot
Ada beberapa alternatif dalam R untuk membuat peta. Dapat menggunakan base plot atau menggunakan levelplot yang dapat diakses dengan fungsi spplot. Untuk membuat peta diperlukan shape file. File ini adalah digital vector dari peta tersebut, berisi beberapa file yaitu .shp, .shx, dan .dbf. Untuk membaca peta pada R menggunakan fungsi readOGR yang berada di library rgdal.
Menampilkan peta dalam bentuk dasar berupa garis-garis wilayah.
columbusMap <- readOGR(system.file("shapes/columbus.shp", package="spData"))## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\User\Documents\R\win-library\4.0\spData\shapes\columbus.shp", layer: "columbus"
## with 49 features
## It has 20 fields
## Integer64 fields read as strings: COLUMBUS_ COLUMBUS_I POLYIDplot(columbusMap) # plot dari base map
Untuk mewarnai peta dengan data yang tersedia dalam data columbus tersebut dapat menggunakan fungsi spplot. Sehingga peta yang tadinya hanya batas wilayah sekarang bisa memmiliki arti. Gradasi warna akan menunjukkan harga-harga perumahan pada data tersebut.Range warna akan menyesuaikan dengan nilai terendah dan tertinggi dari peubah yang kita gunakan.
spplot(columbusMap, "HOVAL", sub="Map of Hoval") 
# HOVAL adalah salah satu peubah di dalam peta Columbus
# dapat diganti dengan peubah lain.
# Sub untuk menambahkan judul di bawahspplot(columbusMap, "HOVAL", sub="Map of Hoval",cut=5) 
# cut untuk memperkecil range warna. Impor Data SHP dari External
setwd("D:\\Kuliah S2 IPB\\Bahan Kuliah\\Semester 1 SSD 2020\\01 Metode Penelitian Kuantitatif\\uts\\spasial\\SPASIAL DENGAN R\\Spatial with R")
# untuk melakukan pengaturan working directory
jabar <- readOGR(dsn = "petajabar27", layer = "Peta Jabar 27") ## OGR data source with driver: ESRI Shapefile
## Source: "D:\Kuliah S2 IPB\Bahan Kuliah\Semester 1 SSD 2020\01 Metode Penelitian Kuantitatif\uts\spasial\SPASIAL DENGAN R\Spatial with R\petajabar27", layer: "Peta Jabar 27"
## with 27 features
## It has 2 fields# dsn adalah nama folder dari peta tersebut
# layer adalah nama petanya
plot(jabar) # peta dasar dalam bentuk batas wilayah
Membuat peta dengan text
plot(jabar)
text(jabar, 'KABKOT', cex=0.5) 
# 'KABKOT' adalah peubah kabkot pada data
# cex untuk menentukan ukuran fontUntuk memberikan warna peta menggunakan pallete untuk membatasi warna apa yang akan dipakai.
palette(rainbow(6))
plot(jabar, col=jabar$IDKAB2) 
#setiap warna sesuai dengan ID-nya
# karena dibatasi hanya 6 warna sehingga ada beberapa warna yang sama pada suatu daerahMengisi Peta dengan Data
colfunc <- colorRampPalette(c("green", "yellow","red"))
#colorRampPalette menentukan warna apa saja yang kita inginkan
jabar$gizi <- data.giziburuk$Giziburuk
spplot(jabar, "gizi",
col.regions = colfunc(6),
cuts = 5,
main = "Peta Sebaran")
Pembobot Spasial
Pembobot spasial terdiri dari pembobot Area/poligon dan Titik. Pembobot untuk area/poligon adalah ketetanggan (contiguity) terdiri dari beberapa cara bishop, rook, dan queen. Pembobot untuk titik adalah jarak(distance), terdiri dari beberapa cara seperti K-nn, Inverse, dan Eksponential.
Penimbang berdasarkan contiguity
Package yang digunakan untuk membuat penimbang spasial adalah spdep.
Queen’s Case
columbus.map <- st_read(system.file("shapes/columbus.shp",
package="spData"),
quiet=TRUE)
queen.w <- poly2nb(as(columbus.map, "Spatial"),
queen = TRUE)
summary(queen.w)## Neighbour list object:
## Number of regions: 49
## Number of nonzero links: 236
## Percentage nonzero weights: 9.829238
## Average number of links: 4.816327
## Link number distribution:
##
## 2 3 4 5 6 7 8 9 10
## 5 9 12 5 9 3 4 1 1
## 5 least connected regions:
## 1 6 42 46 47 with 2 links
## 1 most connected region:
## 20 with 10 linksQueen’s Case – W matrix
queen.w1 <- nb2mat(queen.w,style = "B")Plot the links between the polygons
plot(st_geometry(columbus.map), border="white",col='gray')
coords <- st_coordinates(st_centroid(st_geometry(columbus.map)))
plot(queen.w, coords, add = TRUE, col = "red")
Rook’s Case
rook.w <- poly2nb(as(columbus.map, "Spatial"), queen=FALSE)
plot(st_geometry(columbus.map), border="white",col='gray')
coords <- st_coordinates(st_centroid(st_geometry(columbus.map)))
plot(rook.w, coords, add = TRUE, col = "red")
Plot the links between the polygons
Plot Queen dan Rook W Matrices
plot(st_geometry(columbus.map), border="white",col='grey')
plot(queen.w, coords, add = TRUE, col = "red")
plot(rook.w, coords, add = TRUE, col = "blue")
Penimbang berdasarkan Distance
Pertama, akan mengekstrak koordinat, kemudian dibuat matriks berisi penimbang spasial.
#coordinates(columbus) <-~X+Y #x=long; y=lat
D <- as.matrix(dist(coordinates(columbus),method="euclidean"))
head (D)## 1005 1001 1006 1002 1007 1008 1004 1003
## 1005 0.00000 37.07994 56.56666 53.56437 68.59988 54.10403 24.17551 51.30504
## 1001 37.07994 0.00000 24.49558 25.86193 42.83320 24.66782 40.33898 25.55160
## 1006 56.56666 24.49558 0.00000 16.72850 23.95452 9.76021 59.63632 17.32145
## 1002 53.56437 25.86193 16.72850 0.00000 27.90103 20.07869 55.01842 12.78083
## 1007 68.59988 42.83320 23.95452 27.90103 0.00000 27.32557 72.99820 23.37633
## 1008 54.10403 24.66782 9.76021 20.07869 27.32557 0.00000 55.55063 20.01434
## 1018 1010 1038 1037 1039 1040 1009 1036
## 1005 42.02125 31.75656 98.19590 94.44389 80.82717 86.14240 75.32251 94.86906
## 1001 36.97166 59.52493 84.27394 79.81817 75.25821 81.17488 48.88561 79.79545
## 1006 36.60188 71.85081 66.40066 62.29815 64.33849 69.55758 28.71374 62.15548
## 1002 40.27883 67.13648 71.39838 68.10726 67.99644 72.87520 30.96508 67.38333
## 1007 41.51472 76.64141 58.10954 55.83231 62.21652 64.26446 21.17025 55.11953
## 1008 31.22484 68.67832 65.80145 61.47626 61.59284 67.49313 30.04163 60.61093
## 1011 1042 1041 1017 1043 1019 1012 1035
## 1005 53.91820 82.48801 93.80881 42.01996 99.11878 67.19808 49.74838 74.70231
## 1001 38.73196 84.73298 84.23977 58.95921 86.78280 49.34614 39.69784 73.70668
## 1006 30.58184 79.68763 70.13787 72.82125 72.82245 36.75166 41.12922 68.19724
## 1002 29.78647 83.27702 74.87221 74.73122 79.84725 41.36249 43.81283 70.72295
## 1007 33.87101 80.22679 64.89317 84.71421 68.22090 38.80258 50.89159 68.11019
## 1008 25.96574 77.31596 67.78798 66.37720 69.13314 31.85043 33.67985 63.67520
## 1032 1020 1021 1031 1033 1034 1045 1013
## 1005 98.33540 80.66349 78.80169 94.21592 93.01180 104.54333 96.68331 63.82410
## 1001 82.55910 60.50248 65.00193 79.84132 82.70211 91.28023 88.96477 51.32032
## 1006 64.84498 44.18285 52.36817 63.54707 68.63469 75.39718 81.30780 49.07368
## 1002 69.14075 47.73862 55.74638 67.29627 71.26517 79.31005 87.06008 52.02196
## 1007 54.95650 41.38094 47.17642 55.12256 60.48024 65.28459 85.33352 56.20647
## 1008 63.15269 40.91610 50.88765 60.95906 66.02982 73.84909 76.26345 41.84404
## 1022 1044 1023 1046 1030 1024 1047 1016
## 1005 85.06426 99.35983 83.73819 99.38709 82.67482 94.95609 101.99862 65.00571
## 1001 68.16005 90.32309 68.39659 93.64101 75.89641 79.28130 97.62375 65.75159
## 1006 54.10367 81.03969 56.24367 87.68411 68.09134 64.87389 91.80842 70.45546
## 1002 57.12011 86.70628 60.10677 93.17669 71.29517 67.39760 96.43631 72.72182
## 1007 51.21922 82.64499 54.21017 92.09241 65.90887 58.19684 94.49518 77.50555
## 1008 49.96878 75.99260 51.66738 82.37208 63.45836 61.84989 86.39894 63.88464
## 1014 1049 1029 1025 1028 1048 1015 1027
## 1005 73.45247 104.52995 95.40871 86.84862 93.81061 101.1269 81.46998 96.66114
## 1001 63.70802 102.20867 82.36111 75.88744 81.49586 109.2405 72.72558 84.11485
## 1006 61.98752 97.74379 71.21834 68.82206 72.44941 110.2554 69.25111 75.29042
## 1002 64.29584 103.14758 74.71309 72.26271 76.10283 112.2191 70.68618 78.33910
## 1007 67.97006 100.92889 69.88803 70.75961 72.41371 114.2364 71.03826 76.61968
## 1008 55.11663 91.92015 66.26125 62.71664 66.82763 104.7452 62.96282 69.89376
## 1026
## 1005 91.53183
## 1001 82.31526
## 1006 76.72837
## 1002 79.92787
## 1007 79.00506
## 1008 70.37637Invers Distance Weights
w<-1/D
#row-normalized
diag(w)<-0
rtot<-rowSums(w, na.rm =T)
w_std<-w/rtot
rowSums(w_std, na.rm=T)## 1005 1001 1006 1002 1007 1008 1004 1003 1018 1010 1038 1037 1039 1040 1009 1036
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1011 1042 1041 1017 1043 1019 1012 1035 1032 1020 1021 1031 1033 1034 1045 1013
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1022 1044 1023 1046 1030 1024 1047 1016 1014 1049 1029 1025 1028 1048 1015 1027
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1026
## 1invers.w<-mat2listw(w_std, style="W")
summary(invers.w)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 49
## Number of nonzero links: 2352
## Percentage nonzero weights: 97.95918
## Average number of links: 48
## Link number distribution:
##
## 48
## 49
## 49 least connected regions:
## 1005 1001 1006 1002 1007 1008 1004 1003 1018 1010 1038 1037 1039 1040 1009 1036 1011 1042 1041 1017 1043 1019 1012 1035 1032 1020 1021 1031 1033 1034 1045 1013 1022 1044 1023 1046 1030 1024 1047 1016 1014 1049 1029 1025 1028 1048 1015 1027 1026 with 48 links
## 49 most connected regions:
## 1005 1001 1006 1002 1007 1008 1004 1003 1018 1010 1038 1037 1039 1040 1009 1036 1011 1042 1041 1017 1043 1019 1012 1035 1032 1020 1021 1031 1033 1034 1045 1013 1022 1044 1023 1046 1030 1024 1047 1016 1014 1049 1029 1025 1028 1048 1015 1027 1026 with 48 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 49 2401 49 2.716383 197.1685plot(st_geometry(columbus.map), border="white",col='gray', sub="Inverse Distance")
plot(invers.w, coords, add = TRUE, col = "red")
Eksponential Distance Weights
alpha <- 2
w2 <- exp(-alpha*D)
diag(w2)<-0
rtot<-rowSums(w2, na.rm =T)
rtot## 1005 1001 1006 1002 1007 1008
## 1.003269e-21 1.002146e-21 3.329584e-09 7.922652e-12 4.155522e-19 3.329581e-09
## 1004 1003 1018 1010 1038 1037
## 1.003269e-21 7.920604e-12 1.429511e-23 2.609779e-28 3.809716e-07 4.627262e-07
## 1039 1040 1009 1036 1011 1042
## 3.934315e-10 3.935495e-10 4.090544e-19 8.186012e-08 4.713461e-19 2.736346e-27
## 1041 1017 1043 1019 1012 1035
## 1.329262e-13 3.024577e-34 5.482211e-25 5.901717e-14 3.564140e-14 4.797072e-24
## 1032 1020 1021 1031 1033 1034
## 7.133463e-09 2.160092e-11 3.356299e-22 7.216211e-09 8.276463e-11 4.810802e-15
## 1045 1013 1022 1044 1023 1046
## 2.951847e-09 3.837424e-13 2.197172e-08 1.485294e-09 2.195000e-08 1.907799e-09
## 1030 1024 1047 1016 1014 1049
## 1.341360e-19 1.778688e-13 4.681446e-10 7.055729e-21 1.998307e-12 2.713481e-11
## 1029 1025 1028 1048 1015 1027
## 1.308403e-10 1.478249e-08 5.026807e-09 7.520383e-31 1.650206e-12 4.081598e-09
## 1026
## 1.392096e-08w_std<-w2/rtot
rowSums(w_std, na.rm=T)## 1005 1001 1006 1002 1007 1008 1004 1003 1018 1010 1038 1037 1039 1040 1009 1036
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1011 1042 1041 1017 1043 1019 1012 1035 1032 1020 1021 1031 1033 1034 1045 1013
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1022 1044 1023 1046 1030 1024 1047 1016 1014 1049 1029 1025 1028 1048 1015 1027
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1026
## 1eksp.w=mat2listw(w_std, style="W")
summary(eksp.w)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 49
## Number of nonzero links: 2352
## Percentage nonzero weights: 97.95918
## Average number of links: 48
## Link number distribution:
##
## 48
## 49
## 49 least connected regions:
## 1005 1001 1006 1002 1007 1008 1004 1003 1018 1010 1038 1037 1039 1040 1009 1036 1011 1042 1041 1017 1043 1019 1012 1035 1032 1020 1021 1031 1033 1034 1045 1013 1022 1044 1023 1046 1030 1024 1047 1016 1014 1049 1029 1025 1028 1048 1015 1027 1026 with 48 links
## 49 most connected regions:
## 1005 1001 1006 1002 1007 1008 1004 1003 1018 1010 1038 1037 1039 1040 1009 1036 1011 1042 1041 1017 1043 1019 1012 1035 1032 1020 1021 1031 1033 1034 1045 1013 1022 1044 1023 1046 1030 1024 1047 1016 1014 1049 1029 1025 1028 1048 1015 1027 1026 with 48 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 49 2401 49 68.55016 227.8029plot(st_geometry(columbus.map), border="white",col='gray', sub="Eksponential Distance Weights")
plot(eksp.w, coords, add = TRUE, col = "red")
k-Nearest Neighbor Weights
columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
coords <- st_centroid(st_geometry(columbus), of_largest_polygon=TRUE)
w3=knearneigh(coords, k=1)
knn.w=nb2listw(knn2nb(w3))
summary(knn.w)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 49
## Number of nonzero links: 49
## Percentage nonzero weights: 2.040816
## Average number of links: 1
## Non-symmetric neighbours list
## Link number distribution:
##
## 1
## 49
## 49 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 with 1 link
## 49 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 with 1 link
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 49 2401 49 75 230plot(st_geometry(columbus.map), border="white",col='gray', sub="1-nn")
plot(knn.w, coords, add = TRUE, col = "red")
w2=knearneigh(coords, k=2)
knn.w2=nb2listw(knn2nb(w2))
summary(knn.w2)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 49
## Number of nonzero links: 98
## Percentage nonzero weights: 4.081633
## Average number of links: 2
## Non-symmetric neighbours list
## Link number distribution:
##
## 2
## 49
## 49 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 with 2 links
## 49 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 with 2 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 49 2401 49 42.5 208.5plot(st_geometry(columbus.map), border="white",col='gray', sub="2-nn")
plot(knn.w2, coords, add = TRUE, col = "red")
TERIMA KASIH
Mahasiswa Pascasarjana Statistika dan Sains Data, IPB University, reniamelia@apps.ipb.ac.id↩︎