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 gempa quakes (bawaan R, paket datasets) yang berisi informasi 1000 kejadian gempa di sekitar Kepulauan Fiji sejak tahun 1964. Data ini diubah menjadi objek sf yang mencakup kolom kedalaman gempa (depth, dalam km) dan geometri yang memuat koordinat dari masing-masing kejadian gempa. Koordinat diproyeksikan ke sistem koordinat dengan satuan kilometer.

Langkah-langkah:

Memuat data quakes dan mengubahnya menjadi objek sf.
Membuat variabel baru vble untuk menyimpan kedalaman gempa.

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

# Memuat data gempa
data(quakes)
q <- quakes

# Bujur > 180 diubah ke rentang -180..180
q$long <- ifelse(q$long > 180, q$long - 360, q$long)

# Mengubah menjadi objek sf, lalu diproyeksikan (satuan km) di sekitar Fiji
crs_km <- "+proj=laea +lat_0=-21 +lon_0=179 +datum=WGS84 +units=km +no_defs"
d <- st_as_sf(q, coords = c("long", "lat"), crs = 4326) %>% st_transform(crs_km)

# Menghapus titik dengan koordinat kembar (agar kriging tidak singular)
d <- d[!duplicated(st_coordinates(d)), ]

# Menambahkan kolom kedalaman gempa ke dalam dataset
d$vble <- d$depth

# Membuat peta batas wilayah dari convex hull seluruh titik gempa
map <- st_union(d) %>% st_convex_hull() %>% st_sf()

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

Interpretasi: Data terdiri dari 998 titik gempa dengan koordinat unik. Kedalaman gempa berkisar dari 40 km sampai 680 km, dengan rata-rata 310.8 km dan simpangan baku 215.4 km, sehingga variasinya cukup besar. Titik-titik gempa tidak tersebar merata: titik mengelompok membentuk sabuk memanjang di sisi timur dan lengkung di sisi barat laut, sedangkan bagian tengah dan barat daya relatif kosong. Dari warna titik terlihat bahwa gempa dalam (hijau sampai kuning, sekitar 400 sampai 680 km) terkonsentrasi di sabuk bagian tengah dan pada satu kelompok kecil di barat laut, sedangkan gempa dangkal (ungu) berada di tepi barat dan tepi timur. Hal ini penting karena prediksi di area yang jauh dari titik pengamatan akan kurang andal.

Di sini, map adalah objek yang memuat batas wilayah (convex hull dari seluruh titik gempa), yang nantinya akan digunakan untuk membatasi area perhitungan.

Untuk memprediksi kedalaman, 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 wilayah gempa.
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 wilayah gempa
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 wilayah
coop <- st_filter(coop, map)

# Menampilkan peta lokasi prediksi
qtm(coop)

Interpretasi: Terdapat 5563 lokasi prediksi yang berada di dalam batas wilayah. Karena batas wilayah dibuat dari convex hull seluruh titik gempa, sebagian lokasi prediksi berada di area yang tidak memiliki titik pengamatan di sekitarnya, sehingga prediksi di area tersebut bersifat ekstrapolasi.

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 kedalaman dilakukan berdasarkan nilai titik terdekat.

Langkah-langkah:

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

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

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

Interpretasi: Poligon Voronoi berukuran kecil di area yang padat titik gempa dan berukuran besar di area yang jarang titik (bagian tengah). Seluruh lokasi di dalam satu poligon diberi nilai yang sama, yaitu nilai titik pembentuknya.

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

Interpretasi: Peta Closest Observation berbentuk petak-petak dengan batas yang tajam antarwilayah, karena setiap lokasi hanya mengambil nilai dari satu titik terdekat. Metode ini mempertahankan rentang nilai asli (dari sekitar 0 sampai 700 km), tetapi tidak menghaluskan permukaan sehingga sensitif terhadap pencilan dan tidak memanfaatkan informasi dari tetangga lain. Poligon yang luas di bagian tengah menunjukkan bahwa prediksi di area tersebut hanya bergantung pada satu titik.

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 pada lokasi yang tidak teramati.

library(gstat)

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

# Mendapatkan prediksi 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 dengan metode IDW
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c2 <- tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis")
c2

Interpretasi: Peta IDW paling halus dan sebagian besar nilainya berada pada kisaran 200 sampai 400 km (rentang legend 100 sampai 600 km), sehingga nilai yang sangat dangkal maupun sangat dalam tidak terlihat. Hal ini terjadi karena IDW dengan idp = 1 memakai seluruh titik pengamatan, sehingga titik yang jauh masih memiliki bobot cukup besar dan prediksi cenderung tertarik ke rata-rata global. Akibatnya, IDW kurang mampu menangkap pola lokal (lihat RMSE pada bagian evaluasi).

Nearest Neighbors

Metode Nearest Neighbors menghitung rata-rata nilai 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 nilai 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

Interpretasi: Peta Nearest Neighbors lebih detail daripada IDW dan polanya mirip dengan Closest Observation, tetapi lebih mulus karena prediksi merupakan rata-rata dari 5 titik terdekat. Nilai ekstrem sedikit berkurang dibandingkan Closest Observation, namun pola lokal masih tertangkap dengan baik.

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 dengan masing-masing metode.
Menghitung prediksi ensembel dengan menggabungkan hasil prediksi setiap metode menggunakan bobot tertentu.

# Closest observation (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]

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

## [inverse distance weighted interpolation]

# Ensemble (equal weights)
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

Interpretasi: Peta ensembel (bobot sama, masing-masing 1/3) merupakan kompromi dari tiga metode sebelumnya. Pola utama dari Closest Observation dan Nearest Neighbors masih terlihat, tetapi diperhalus oleh IDW sehingga rentang nilainya menyempit. Perlu diperhatikan bahwa peta ini memakai bobot sama, sedangkan pada evaluasi cross validation bobot ensembel ditentukan dari invers RMSE.

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

Interpretasi: Urutan peta: kiri atas Closest Observation, kanan atas IDW, kiri bawah Nearest Neighbors, dan kanan bawah Ensembel. Keempatnya memperlihatkan pola spasial yang secara umum konsisten: kedalaman besar berada di bagian barat laut dan di sabuk bagian tengah yang memanjang dari utara ke selatan, sedangkan kedalaman dangkal berada di tepi barat dan tepi timur. Perbedaan utama ada pada tingkat kehalusan permukaan: Closest Observation paling kasar (berpetak), sedangkan IDW paling halus dan paling mendekati rata-rata.

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)

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

mapview(coop, map.types = "CartoDB.Voyager")

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(vble ~ 1, d, cloud = TRUE)
plot(vc)

# Variogram dengan plot standar
v <- variogram(vble ~ 1, data = d)
plot(v)

Interpretasi: Semivarian rendah pada jarak dekat, meningkat seiring bertambahnya jarak, lalu mendatar pada jarak sekitar 1000 km. Pola ini menunjukkan adanya ketergantungan spasial yang kuat: gempa yang lokasinya berdekatan cenderung memiliki kedalaman yang mirip, sedangkan gempa yang berjauhan semakin berbeda. Kondisi ini mendukung penggunaan Kriging.

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
# (nilai awal diambil dari variogram empiris v karena skala data berbeda)
psill0 <- max(v$gamma) - min(v$gamma)
range0 <- max(v$dist) / 2
nugget0 <- min(v$gamma)
vinitial <- vgm(psill = psill0, model = "Sph", range = range0, nugget = nugget0)

# Memvisualisasikan variogram dengan model awal (vinitial)
plot(v, vinitial, cex = 1.5)

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

##   model    psill    range
## 1   Nug     0.00    0.000
## 2   Sph 69471.15 1068.278

plot(v, fv, cex = 1.5)

Interpretasi: Model Spherical yang di-fit memiliki nugget sebesar 0, partial sill sebesar 69471, dan range sebesar 1068 km. Nugget yang sangat kecil berarti hampir tidak ada variasi mikrospasial atau galat pengukuran, sehingga kedalaman gempa pada lokasi yang sangat berdekatan hampir identik. Sill sebesar 69471 menunjukkan total varians spasial, dan range sebesar 1068 km menunjukkan bahwa gempa yang berjarak kurang dari itu masih saling berkorelasi, sedangkan di luar jarak tersebut dianggap tidak saling bergantung. Kurva model mengikuti titik-titik variogram empiris dengan cukup baik.

library(automap)
# Autofit variogram untuk data kedalaman
afvg <- autofitVariogram(vble ~ 1, d)

# Plot hasil
plot(afvg)

Interpretasi: Fungsi autofitVariogram dipakai sebagai pembanding dan secara otomatis memilih model Gaussian dengan nugget 1355, sill 64762, dan range 389. Nilai sill-nya mirip dengan model Spherical di atas (sekitar 65.000 sampai 70.000), sehingga kedua model memberikan gambaran struktur spasial yang serupa. Pada tugas ini prediksi Kriging tetap memakai model Spherical hasil fit.variogram seperti pada materi.

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

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

## [using ordinary kriging]

# Visualisasi prediksi kedalaman dengan warna berdasarkan var1.pred
ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
  geom_sf(data = d) +
  scale_color_viridis(name = "Kedalaman (km)") + theme_bw()

Interpretasi (peta prediksi): Peta prediksi Kriging berupa permukaan yang mulus dan kontinu. Kedalaman besar (sekitar 500 sampai 680 km, warna kuning) terdapat di bagian barat laut, di bagian tengah, dan di sabuk memanjang ke arah selatan, sedangkan kedalaman dangkal (warna ungu) berada di tepi barat dan tepi timur. Pola ini sejalan dengan peta Closest Observation dan Nearest Neighbors, tetapi tanpa batas petak yang tajam.

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

Interpretasi (peta varians): Varians Kriging rendah (warna ungu) di sekitar titik pengamatan dan meningkat pada area yang jauh dari titik. Varians terbesar berada di bagian barat daya yang hampir tidak memiliki titik gempa, dengan nilai mendekati sill (sekitar 69.000). Artinya prediksi di area tersebut paling tidak pasti. Ini adalah keunggulan Kriging dibandingkan metode deterministik, karena selain memberikan prediksi juga memberikan ukuran ketidakpastian.

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
rmse5 <- rep(NA, 5) # Kriging

for(k in 1:5){
  # Split data into test and train
  test <- d[kf == k, ]
  train <- d[kf != k, ]

  # Closest observation
  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)

  # Nearest neighbors
  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 (weights are inverse RMSE so lower RMSE higher weight)
  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)

  # Kriging (variogram di-fit ulang pada data latih)
  vk <- variogram(vble ~ 1, data = train)
  fvk <- suppressWarnings(fit.variogram(object = vk,
           model = vgm(psill = max(vk$gamma) - min(vk$gamma), model = "Sph",
                       range = max(vk$dist) / 2, nugget = min(vk$gamma))))
  kr <- gstat(formula = vble ~ 1, data = train, model = fvk)
  p5 <- predict(kr, test)$var1.pred
  rmse5[k] <- RMSE(test$vble, p5)
}

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]

## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]

# RMSE obtained for each of the 5 splits
data.frame(closest.obs = rmse1, IDW = rmse2,
           nearest.neigh = rmse3, ensemble = rmse4, kriging = rmse5)

##   closest.obs      IDW nearest.neigh ensemble  kriging
## 1    78.05844 130.8342      56.60994 63.44850 61.10442
## 2    58.11179 133.4721      61.27105 55.01635 51.31094
## 3    83.55971 147.5291      76.25051 78.87651 76.25196
## 4    80.15431 133.5280      74.89965 72.02851 69.83357
## 5    84.55318 140.2018      64.94348 70.70188 74.90035

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

##   closest.obs     IDW nearest.neigh ensemble  kriging
## 1    76.88749 137.113      66.79493 68.01435 66.68025

Interpretasi Hasil Evaluasi

Interpretasi: RMSE rata-rata dari 5-fold cross validation (dalam km) adalah: Kriging 66.68, Nearest Neighbors 66.79, Ensembel 68.01, Closest Observation 76.89, dan IDW 137.11. Sebagai pembanding, simpangan baku kedalaman gempa adalah 215.4 km, sehingga metode terbaik (RMSE 66.68 km) menghasilkan galat sekitar 69% lebih kecil dibandingkan hanya menebak dengan rata-rata.

Kriging memiliki RMSE terkecil karena bobot tetangga ditentukan dari struktur ketergantungan spasial (variogram), bukan hanya dari jarak.
Nearest Neighbors hampir setara dengan Kriging (selisih sekitar 0.11 km). Pada beberapa fold NN bahkan lebih baik, sehingga keduanya dapat dianggap setara.
Ensembel lebih baik daripada Closest Observation dan IDW, tetapi sedikit di bawah NN dan Kriging. Hal ini karena ensembel masih memuat IDW yang galatnya besar, dan Kriging tidak termasuk dalam ensembel.
Closest Observation cukup baik, tetapi lebih buruk daripada NN karena hanya bergantung pada satu titik.
IDW memiliki RMSE terbesar karena pengaturan idp = 1 dengan seluruh titik membuat prediksi terlalu halus dan mendekati rata-rata global.

Kesimpulan

Berdasarkan hasil 5-fold cross validation, metode dengan RMSE rata-rata terkecil adalah kriging (RMSE = 66.68 km), sedangkan RMSE terbesar dihasilkan oleh IDW (RMSE = 137.11 km). Semakin kecil RMSE, semakin baik kemampuan metode dalam memprediksi kedalaman pada lokasi yang tidak teramati.

Kedalaman gempa di sekitar Fiji memiliki ketergantungan spasial yang kuat. Hal ini terlihat dari variogram yang naik terhadap jarak hingga mencapai sill pada range sekitar 1068 km.
Metode yang memanfaatkan informasi lokal (Kriging dan Nearest Neighbors) memberikan prediksi paling akurat, sedangkan IDW dengan idp = 1 yang memakai seluruh titik paling kurang akurat.
Kriging memiliki keunggulan tambahan berupa peta varians yang menunjukkan tingkat ketidakpastian prediksi, terutama di area yang jauh dari titik pengamatan.
Keterbatasan: titik gempa mengelompok dan terdapat area yang kosong, sehingga prediksi di area tanpa data bersifat ekstrapolasi dan kurang andal. Selain itu, pembagian fold secara acak dapat membuat titik uji berdekatan dengan titik latih, sehingga RMSE cenderung sedikit optimistis.

Pertemuan 9: Interpolasi Spasial

Aurelia Fury - M0401241029

Interpolasi Spasial

Konsep Dasar

Data yang digunakan

Metode Deteministik

Closest Observation

Invers Distance Weight

Nearest Neighbors

Pendekatan Ensembel

Metode Stokastik

Kriging

Evaluasi Model

Interpretasi Hasil Evaluasi

Kesimpulan