Pertemuan 9: Interpolasi Spasial

Author

Rosita Ria Rusesta

knitr::opts_chunk$set(warning = FALSE, message = FALSE)

Interpolasi Spasial

Interpolasi spasial merupakan metode yang digunakan untuk memperkirakan nilai peubah spasial pada lokasi yang tidak diketahui berdasarkan data yang tersedia pada lokasi tetangga yang diketahui (Cressie, 1993). Interpolasi ini digunakan dalam berbagai bidang seperti geografi, meteorologi, pertanian, dan lingkungan untuk menghasilkan peta atau permukaan kontinu dari data yang awalnya hanya tersedia pada titik-titik diskrit.

Konsep Dasar

Interpolasi spasial mengasumsikan bahwa nilai variabel yang diukur (misalnya, suhu, ketinggian, kelembaban tanah, dll.) di suatu lokasi akan serupa dengan nilai di lokasi-lokasi tetangga. Dengan kata lain, ada korelasi spasial, di mana lokasi yang lebih dekat cenderung memiliki nilai yang lebih mirip dibandingkan lokasi yang lebih jauh.

Data yang digunakan

Dalam bagian ini, kita akan menggunakan data harga properti dari paket spData (Bivand et al., 2022) yang berisi informasi tentang harga apartemen di Athena, Yunani, pada tahun 2017. Data yang digunakan adalah objek sf yang mencakup beberapa kolom seperti harga per meter persegi (prpsqm) dan geometri yang memuat koordinat geografis dari masing-masing apartemen.

Langkah-langkah:

Memuat data dari paket spData dan mengubahnya menjadi objek sf.
Membuat variabel baru vble untuk menyimpan harga properti per meter persegi.

# Memuat pustaka yang diperlukan
library(spData)
library(sf)
library(terra)
library(tmap)
library(viridis)
library(dismo)
library(dplyr)

# Memuat data properti
d <- properties

# Menambahkan kolom harga per meter persegi ke dalam dataset
d$vble <- d$prpsqm

# Membuat peta dari batas administratif Athena
map <- st_union(depmunic) %>% st_sf()

# Menampilkan peta harga properti
tmap_mode("plot")
tm_shape(map) + tm_polygons(alpha = 0.3) + tm_shape(d) +
tm_dots("vble", palette = "viridis")

Di sini, depmunic adalah objek yang memuat batas administratif Athena, yang nantinya akan digunakan untuk membatasi area perhitungan.

Untuk memprediksi harga properti di Athena, kita akan membuat grid raster yang mencakup seluruh area dan menggunakan titik pusat grid sebagai lokasi prediksi.

Langkah-langkah:

Membuat grid raster dengan dimensi yang cukup kecil (misalnya 100x100 sel) untuk menutupi area Athena.
Mengambil koordinat dari setiap sel pada grid raster.
Mengubah koordinat-koordinat sel grid menjadi objek sf yang akan digunakan untuk prediksi.

# Membuat grid raster 100x100 sel untuk menutupi Athena
grid <- rast(map, nrows = 100, ncols = 100)

# Mengambil koordinat titik pusat dari setiap sel dalam grid
xy <- xyFromCell(grid, 1:ncell(grid))

# Membuat objek sf untuk lokasi prediksi
coop <- st_as_sf(as.data.frame(xy), coords = c("x", "y"), crs = st_crs(map))

# Menyaring lokasi prediksi yang berada dalam batas Athena
coop <- st_filter(coop, map)

# Menampilkan peta lokasi prediksi
qtm(coop)

Metode Deteministik

Closest Observation

Metode ini menggunakan diagram Voronoi untuk membagi ruang menjadi poligon yang masing-masing berisi satu titik penghasil. Setiap titik dalam poligon lebih dekat ke titik penghasil daripada ke titik lainnya. Prediksi harga properti dilakukan berdasarkan nilai titik terdekat.

Langkah-langkah:

Membuat diagram Voronoi berdasarkan data titik properti.
Menggunakan st_intersection() untuk menentukan prediksi harga pada titik-titik yang terletak dalam setiap poligon Voronoi.

# Membuat diagram Voronoi dari lokasi properti
v <- voronoi(x = vect(d), bnd = map)
plot(v)
# Menampilkan diagram Voronoi dengan harga properti
points(vect(d), cex = 0.5)

# Prediction
v <- st_as_sf(v)
tm_shape(v) + tm_fill(col = "vble", palette = "viridis", alpha = 0.6)

Di sini, kita dapat melihat bagaimana setiap poligon dalam diagram Voronoi mewakili wilayah prediksi yang dipengaruhi oleh titik terdekat.

resp <- st_intersection(v, coop)
resp$pred <- resp$vble
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c1<-tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c1

Invers Distance Weight

Metode IDW (Inverse Distance Weighting) menggunakan bobot berdasarkan jarak antara lokasi yang belum teramati dan lokasi-lokasi yang sudah teramati. Semakin dekat jarak, semakin besar bobot yang diberikan pada nilai lokasi yang teramati.

Langkah-langkah:

Menerapkan IDW dengan menggunakan fungsi gstat().
Menggunakan predict() untuk mendapatkan prediksi harga pada lokasi yang tidak teramati.

library(gstat)

# Menerapkan IDW menggunakan semua lokasi properti
res <- gstat(formula = vble ~ 1, locations = d, nmax = nrow(d), set = list(idp = 1))

# Mendapatkan prediksi harga untuk lokasi yang belum teramati
resp <- predict(res, coop)

[inverse distance weighted interpolation]

resp$x <- st_coordinates(resp)[,1]
resp$y <- st_coordinates(resp)[,2]
resp$pred <- resp$var1.pred

# Membuat peta prediksi harga dengan metode IDW
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c2<-tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c2

Nearest Neighbors

Metode Nearest Neighbors menghitung rata-rata harga dari k titik terdekat. Setiap tetangga dianggap memiliki bobot yang sama.

Langkah-langkah:

Menggunakan fungsi gstat() untuk menerapkan nearest neighbors dengan parameter nmax yang menentukan jumlah titik tetangga.
Menggunakan fungsi predict() untuk memprediksi harga di lokasi yang tidak teramati.

res <- gstat(formula = vble ~ 1, locations = d, nmax = 5, set = list(idp = 0))

resp <- predict(res, coop)

[inverse distance weighted interpolation]

resp$x <- st_coordinates(resp)[,1]
resp$y <- st_coordinates(resp)[,2]

resp$pred <- resp$var1.pred

pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c3<-tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c3

Pendekatan Ensembel

Pendekatan Ensembel menggabungkan hasil prediksi dari berbagai metode interpolasi (seperti Closest Observation, IDW, dan Nearest Neighbors) dengan memberikan bobot pada setiap metode. Bobot ini bisa berdasarkan kinerja masing-masing metode.

Langkah-langkah:

Menghitung prediksi harga dengan masing-masing metode.
Menghitung prediksi ensembel dengan menggabungkan hasil prediksi setiap metode menggunakan bobot tertentu.

#Closestobservation(Voronoi)
v<-voronoi(x=vect(d),bnd= map)
v<-st_as_sf(v)
p1<-st_intersection(v,coop)$vble

#IDW
gs<-gstat(formula=vble~1,locations= d,nmax= nrow(d),
set= list(idp=1))
p2<-predict(gs,coop)$var1.pred

[inverse distance weighted interpolation]

#Nearestneighbors
nn<-gstat(formula=vble~1,locations= d,nmax= 5,
set= list(idp=0))
p3<-predict(nn,coop)$var1.pred

[inverse distance weighted interpolation]

#Ensemble(equalweights)
weights<-c(1/3,1/3,1/3)
p4<-p1*weights[1]+p2*weights[2]+p3*weights[3]

resp<-data.frame(x=st_coordinates(coop)[,1], y=st_coordinates(coop)[,2], pred=p4)

resp<-st_as_sf(resp,coords= c("x","y"),crs= st_crs(map))
pred<-rasterize(resp,grid, field= "pred",fun= "mean")
c4<-tm_shape(pred)+tm_raster(alpha= 0.6,palette="viridis")
c4

# Gabungkan semua plot dalam satu frame
tmap_arrange(c1, c2, c3, c4, ncol = 2)

Evaluasi Model

Untuk menilai kinerja setiap metode, kita melakukan K-fold cross-validation dan menghitung RMSE (Root Mean Square Error). Cross-validation digunakan untuk memastikan bahwa model tidak hanya cocok dengan data latih tetapi juga dapat memprediksi data baru dengan baik.

Langkah-langkah:

Membagi data menjadi K bagian dengan fungsi dismo::kfold().
Menghitung RMSE untuk setiap metode dengan membandingkan nilai prediksi dan nilai sebenarnya.

#install.packages("dismo")

set.seed(123)

# Fungsi untuk menghitung RMSE
RMSE <- function(observed, predicted) {
  sqrt(mean((observed - predicted)^2))
}

# Pembagian data untuk K-fold
kf <- kfold(nrow(d), k = 5)

# Vectors to store the RMSE values obtained with each method
rmse1 <- rep(NA, 5) # Closest observation
rmse2 <- rep(NA, 5) # IDW
rmse3 <- rep(NA, 5) # Nearest neighbors
rmse4 <- rep(NA, 5) # Ensemble
for(k in 1:5){
#Splitdataintest andtrain
test<-d[kf== k,]
train<-d[kf !=k,]


#Closestobservation
v<-voronoi(x=vect(train), bnd= map)
v<-st_as_sf(v)
p1<-st_intersection(v,test)$vble
rmse1[k] <-RMSE(test$vble,p1)

#IDW
gs<-gstat(formula=vble ~1,locations= train,
nmax=nrow(train), set= list(idp=1))
p2 <-predict(gs,test)$var1.pred
rmse2[k] <-RMSE(test$vble,p2)

#Nearestneighbors
nn<-gstat(formula=vble~1,locations= train,
nmax=5,set= list(idp=0))
p3 <-predict(nn,test)$var1.pred
rmse3[k] <-RMSE(test$vble,p3)

#Ensemble(weightsareinverseRMSEsolowerRMSEhigherweight)
w<-1/c(rmse1[k],rmse2[k],rmse3[k])
weights<-w/sum(w)
p4<-p1*weights[1]+p2*weights[2]+p3*weights[3]
rmse4[k] <-RMSE(test$vble,p4)
}

[inverse distance weighted interpolation]
[inverse distance weighted interpolation]

[inverse distance weighted interpolation]
[inverse distance weighted interpolation]

[inverse distance weighted interpolation]
[inverse distance weighted interpolation]

[inverse distance weighted interpolation]
[inverse distance weighted interpolation]

[inverse distance weighted interpolation]
[inverse distance weighted interpolation]

#RMSEobtainedforeachofthe5splits
data.frame(closest.obs= rmse1,IDW= rmse2,
nearest.neigh= rmse3,ensemble= rmse4)

  closest.obs      IDW nearest.neigh ensemble
1    857.0022 855.0812      793.3416 747.2656
2    857.5906 762.6661      736.2985 703.3591
3   1059.6057 962.7646      897.8722 868.1311
4    994.6159 921.5928      861.9015 819.1503
5    792.8463 844.1711      773.5493 714.6230

# Average RMSE over the 5 splits
data.frame(closest.obs = mean(rmse1), IDW = mean(rmse2), nearest.neigh = mean(rmse3), ensemble = mean(rmse4))

  closest.obs      IDW nearest.neigh ensemble
1    912.3322 869.2552      812.5926 770.5058

Metode Stokastik

Kriging

Kriging adalah metode interpolasi geostatistik yang memperhitungkan autokorelasi spasial untuk membuat prediksi. Ordinary Kriging (OK) adalah salah satu varian Kriging yang paling umum digunakan, di mana nilai yang akan diprediksi adalah variabel z yang tidak diketahui, yang dianggap sebagai variabel acak.

Variogram: Dalam Ordinary Kriging, penting untuk membangun variogram. Variogram adalah grafik yang menggambarkan bagaimana variabilitas antar titik pengamatan berubah seiring dengan jarak antar titik tersebut. Kriging menggunakan informasi variogram ini untuk menghitung covariance antar titik.
Model Semivariogram: Setelah variogram dihitung, model semivariogram (seperti spherical, exponential, atau Gaussian) di-fit ke data untuk menangkap hubungan spasial di antara titik-titik pengamatan.
Ordinary Kriging mengasumsikan bahwa nilai rata-rata populasi tidak diketahui, tetapi konstan di seluruh wilayah.

# Memuat paket yang dibutuhkan
library(sp)        
library(gstat)     
library(sf)        
library(mapview)   

# Memuat dataset meuse dan meuse.grid yang sudah tersedia di R
data(meuse)        
data(meuse.grid)   

# Mengonversi dataset meuse menjadi objek 'sf' (simple feature) dengan koordinat 'x' dan 'y', dan CRS 28992 (proyeksi UTM)
meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
meuse.grid <- st_as_sf(meuse.grid, coords = c("x", "y"), crs = 28992)

# Menampilkan peta interaktif 
mapview(meuse, zcol = "zinc", map.types = "CartoDB.Voyager")

mapview(meuse.grid, map.types = "CartoDB.Voyager")

Sumber : kodu.ut.ee (Lesson 6 : Spatial Interpolation - Geospatial Analysis with Python and R

Keterangan :

Nugget: Nilai semivarian ketika jarak antarobservasi mendekati nol, yang merepresentasikan variasi mikrospasial dan kesalahan pengukuran.
Sill: Nilai maksimum atau plateau semivarian yang merepresentasikan total varians spasial.
Range: Jarak ketika semivarian mencapai sill, yang menunjukkan batas jangkauan ketergantungan spasial.

Secara umum:

Semivarian rendah pada jarak dekat → indikasi bahwa lokasi berdekatan memiliki nilai yang mirip.
Semivarian meningkat seiring jarak → indikasi adanya struktur ketergantungan spasial.
Semivarian mencapai sill → variasi spasial relatif stabil setelah jarak tertentu.

# Variogram dengan cloud plot
vc <- variogram(log(zinc) ~ 1, meuse, cloud = TRUE)
plot(vc)

# Variogram dengan plot standar
v <- variogram(log(zinc) ~ 1, data = meuse)
plot(v)

vgm()

   short                                      long
1    Nug                              Nug (nugget)
2    Exp                         Exp (exponential)
3    Sph                           Sph (spherical)
4    Gau                            Gau (gaussian)
5    Exc        Exclass (Exponential class/stable)
6    Mat                              Mat (Matern)
7    Ste Mat (Matern, M. Stein's parameterization)
8    Cir                            Cir (circular)
9    Lin                              Lin (linear)
10   Bes                              Bes (bessel)
11   Pen                      Pen (pentaspherical)
12   Per                            Per (periodic)
13   Wav                                Wav (wave)
14   Hol                                Hol (hole)
15   Log                         Log (logarithmic)
16   Pow                               Pow (power)
17   Spl                              Spl (spline)
18   Leg                            Leg (Legendre)
19   Err                   Err (Measurement error)
20   Int                           Int (Intercept)

# Menampilkan daftar model variogram
show.vgms(par.strip.text = list(cex = 0.75))

# Menentukan model variogram awal dengan parameter psill, model, range, dan nugget
vinitial <- vgm(psill = 0.5, model = "Sph", range = 900, nugget = 0.1)

# Memvisualisasikan variogram dengan model awal (vinitial) dan batas cutoff 1000
plot(v, vinitial, cutoff = 1000, cex = 1.5)

# Menyesuaikan variogram dengan model "Sph" menggunakan data v dan parameter
fv<-fit.variogram(object=v,model= vgm(psill= 0.5,model= "Sph",range= 900,nugget= 0.1))
fv

  model      psill    range
1   Nug 0.05066017   0.0000
2   Sph 0.59060556 897.0044

plot(v,fv, cex= 1.5)

library(automap)
# Autofit variogram untuk data log-zinc
afvg <- autofitVariogram(log(zinc) ~ 1, meuse)

# Plot hasil
plot(afvg)

library(ggplot2)
library(viridis)
# Membuat objek gstat dengan model variogram hasil fitting (fv) untuk data zinc
k <- gstat(formula = log(zinc) ~ 1, data = meuse, model = fv)

# Melakukan prediksi variabel pada grid menggunakan model gstat yang sudah dibuat
kpred <- predict(k, meuse.grid)

[using ordinary kriging]

# Visualisasi prediksi log(zinc) dengan warna berdasarkan var1.pred
ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = meuse) +
  scale_color_viridis(name = "log(zinc)") + theme_bw()

# Visualisasi prediksi varians dengan warna berdasarkan var1.var
ggplot() + geom_sf(data = kpred, aes(color = var1.var)) +
  geom_sf(data = meuse) +
  scale_color_viridis(name = "variance") + theme_bw()

Tugas:

Buat interpolasi menggunakan data meuse, dengan menggunakan seluruh metode interpolasi, yaitu Closest Observation, NN, IDW, Ensembel (Bebas menggunakan metode interpolasi apa saja yang ingin di ensemble, minimal 3), dan Kriging. setelah itu bandingkan RMSE dengan menggunakan Cross Validation