Analisis Spasial

Pendahuluan

Jenis Data Spasial

Data titik (point patterns): Lokasi suatu kejadian (lokasi pohon-pohon di hutan, lokasi sarang burung pada pohon)
Data kontinu (Continus data): kejadian-kejadian yang bersifat kontinu dalam suatu ruang (pola curah hujan, pola salinitas air tanah di Jawa Barat)
Data Area: kejadian dipisahkkan dalam zona (kepadatan penduduk kabupaten/ kota, angka kesakitan malaria kab/kota di Indonesia)

Tipe 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 earth. Data raster terdiri dari Rasterlayer dan Rasterstack.

Packages

Packages yang digunakan adalah sp, rgdal, spData, readxl, raster, spdep

library(sp)
library(rgdal) # membaca shape file
library(spData)
library(readxl) # membaca file excel
library(raster) # text pada peta
library(spdep) # pembobot data spasial

Data

Data dalam `package` R

Untuk melakukan uji coba program dapat memanfaatkan dataset yang tersedia di dalam package R sesuai tema latihan. Sebagai contoh, untuk data spasial ini dapat menggunakan data yang berada pada package spData.

ls() # melihat dataset yang sedang dijalankan pada R

character(0)

Jika dikembalikan dengan character(0) berarti tidak ada dataset yang sedang dijalankan.

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"

Jika fungsi ls sudah memberikan output tandanya data sudah berhasil dijalankan pada R. Untuk melihat struktur data dapat menggunakan fungsi str.

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 merupakan class data.frame terdiri dari 49 observasi dengan 22 peubah. Diberikan penjelasan untuk masing-masing variabel berupa tipe datanya. Untuk angka angka yang muncul adalah cuplikan data dari masing-masing peubah. Sebagai contoh peubah AREA tipe datanya adalah NUMERIC dengan cuplikan sebagian datanya adalah 10.39, 8.62, dan seterusnya.

Untuk melihat deskripsi dari dataset dapat menggunakan fungsi help

help(columbus)

Berikut adalah penjelasan dari dataset columbus yang akan digunakan untuk latihan.

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.

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:\\SD\\") # untuk melakukan pengaturan working directory
data.giziburuk <- read_excel("data giziburuk jabar.xlsx", sheet = "Sheet1") # memasukkan file excel ke dalam objek data.giziburuk

Untuk melihat dimensi dari data dapat menggunakan fungsi dim. Terlihat bahwa dimensi dari data tersebut adalah 27x5.

dim(data.giziburuk)

[1] 27  5

Untuk melihat struktur data dapat menggunakan fungsi str. Terlihat bahwa data gizi buruk ini memiliki class 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. Fungsi tersebut akan memunculkan sebuah layar baru berisi tabel tersebut. Bentuk tabel tersebut menyerupai tampilan tabel pada excel.

View(data.giziburuk)

Plot

Ada beberapa alternatif dalam R untuk membuat peta, seperti menggunakan base plot atau menggunakan levelplot yang dapat diakses dengan fungsi spplot. Untuk membuat peta diperlukan shape file. Shape file merupakan 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\ANNISARAHMI\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 POLYID

plot(columbusMap) # plot dari base map

Untuk mewarnai peta dengan data yang tersedia dalam data columbus tersebut dapat menggunakan fungsi spplot. Peta yang sebelumnya hanya terdiri dari batas wilayah saja bisa memmiliki arti ketika diberikan perbedaan warna sesuai dengan data yang berada di dalam peta tersebut. Sebagai contoh akan dibuat peta yang menampilkan informasi harga-harga perumahan. Gradasi warna akan menunjukkan harga-harga perumahan pada data tersebut.

spplot(columbusMap, "HOVAL", sub="Map of Hoval") # HOVAL adalah salah satu peubah di dalam peta, dapat diganti dengan peubah lain. Sub untuk menambahkan judul di bawah

Impor Data SHP dari External
Jika sebelumnya kita menggunakan shp yang sudah tersedia dalam package, selanjutnya kita juga dapat menggunakan file shp yang berada di luar package. Untuk membaca file tersebut, kita lakukan setwd terlebih dahulu jika belum dilakukan, jika sudah maka tinggal menuliskan dsn-nya. dsn adalah nama folder dari peta tersebut. layer adalah nama shp dari peta yang akan ditampilkan.

jabar <- readOGR(dsn = "petajabar27", layer = "Peta Jabar 27") # dsn adalah nama folder dari oeta tersebut, layer adalah nama petanya

OGR data source with driver: ESRI Shapefile 
Source: "D:\SD\petajabar27", layer: "Peta Jabar 27"
with 27 features
It has 2 fields

plot(jabar) # peta dasar dalam bentuk batas wilayah

Membuat peta dengan text
Untuk menambahkan text di dalam peta, dapat menggunakan fungsi text. Untuk mengatur ukuran dari text yang ditampilkan dapat menggunakan fungsi cex.

plot(jabar)
text(jabar, 'KABKOT', cex=0.5) # 'KABKOT' adalah peubah kabkot pada data, cex untuk menentukan ukuran fonr

Memberi warna pada peta
Untuk memberikan warna dapat menggunakan fungsi pallete. Fungsi ini dapat membatasi warna apa yang akan dipakai.

palette(rainbow(6))
plot(jabar, col=jabar$IDKAB2) #setiap warna sesuai dengan IDnya, karena dibatasi hanya 6 warna sehingga ada beberapa warna yang sama pada suatu daerah

Mengisi Peta dengan Data
Untuk menampilan peta yang berisi data mengenai peta tersebut maka dapat menggunakan fungsi

colfunc <- colorRampPalette(c("green", "yellow","red"))
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. Package yang digunakan untuk membuat penimbang spasial adalah spdep.

Penimbang berdasarkan contiguity

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 links

Queen’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
         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
         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
         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
         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
         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
         1026
1005 91.53183
 [ reached getOption("max.print") -- omitted 5 rows ]

Invers 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 
   1

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

plot(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-08

w_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 
   1

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

plot(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 230

plot(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.5

plot(st_geometry(columbus.map), border="white",col='gray', sub="2-nn")
plot(knn.w2, coords, add = TRUE, col = "red")

Referensi

Djuraidah, Anik. (2021). STA500 Metode Penelitan Kuantitatif: Spasial. Retrieved from https://newlms.ipb.ac.id/

Hujjah A, Al. (2021). Spatial with R. Retrieved from https://newlms.ipb.ac.id/