Pertemuan 9: Interpolasi Spasial

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

Pada bagian ini digunakan dataset nz_height yang tersedia dalam paket spData pada R. Dataset ini berisi informasi ketinggian permukaan tanah pada sejumlah titik pengamatan di wilayah Selandia Baru (New Zealand), yang dilengkapi dengan informasi koordinat spasial pada setiap lokasi pengamatan.

Variabel yang digunakan dalam proses interpolasi adalah ketinggian (elevation). Variabel ini merupakan data numerik kontinu yang menunjukkan ketinggian suatu lokasi dalam satuan meter. Variabel tersebut digunakan sebagai target interpolasi untuk memperkirakan ketinggian pada lokasi yang tidak memiliki pengamatan.

Langkah-langkah: 1. Memuat dataset nz_height dari paket spData. 2. Memilih variabel elevation sebagai variabel yang akan diinterpolasi. 3. Membentuk variabel baru vble yang berisi nilai ketinggian. 4. Mengubah data menjadi objek spasial menggunakan koordinat lokasi pengamatan. 5. Membuat batas wilayah pengamatan berdasarkan convex hull dari sebaran titik observasi.

library(sp)
library(spData)
library(sf)
library(terra)
library(tmap)
library(viridis)
library(dismo)
library(dplyr)
library(gstat)

# Memuat dataset
data("nz_height", package = "spData")

head(nz_height)
## Simple feature collection with 6 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1204143 ymin: 5048309 xmax: 1389460 ymax: 5168749
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
##   t50_fid elevation                geometry
## 1 2353944      2723 POINT (1204143 5049971)
## 2 2354404      2820 POINT (1234725 5048309)
## 3 2354405      2830 POINT (1235915 5048745)
## 4 2369113      3033 POINT (1259702 5076570)
## 5 2362630      2749 POINT (1378170 5158491)
## 6 2362814      2822 POINT (1389460 5168749)
str(nz_height)
## Classes 'sf' and 'data.frame':   101 obs. of  3 variables:
##  $ t50_fid  : int  2353944 2354404 2354405 2369113 2362630 2362814 2362817 2363991 2363993 2363994 ...
##  $ elevation: int  2723 2820 2830 3033 2749 2822 2778 3004 3114 2882 ...
##  $ geometry :sfc_POINT of length 101; first list element:  'XY' num  1204143 5049971
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
##   ..- attr(*, "names")= chr [1:2] "t50_fid" "elevation"
# Menyalin dataset
d <- nz_height

# Membentuk variabel target interpolasi
d$vble <- d$elevation

# Mengubah menjadi objek sf
d <- st_as_sf(d)

# Memeriksa CRS
st_crs(d)
## Coordinate Reference System:
##   User input: EPSG:2193 
##   wkt:
## PROJCRS["NZGD2000 / New Zealand Transverse Mercator 2000",
##     BASEGEOGCRS["NZGD2000",
##         DATUM["New Zealand Geodetic Datum 2000",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4167]],
##     CONVERSION["New Zealand Transverse Mercator 2000",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",173,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",1600000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",10000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["New Zealand - North Island, South Island, Stewart Island - onshore."],
##         BBOX[-47.33,166.37,-34.1,178.63]],
##     ID["EPSG",2193]]
# Membuat batas wilayah pengamatan (Convex Hull)
map <- st_sf(
  geometry = st_convex_hull(
    st_union(st_geometry(d))
  )
)

# Visualisasi
tmap_mode("plot")

tm_shape(map) +
  tm_borders(col = "black") +
  tm_shape(d) +
  tm_dots(
    col = "vble",
    palette = "viridis",
    size = 0.01,
    title = "Elevation (m)"
  ) +
  tm_layout(
    title = "Sebaran Ketinggian pada Dataset nz_height",
    legend.outside = TRUE
  )

Pada bagian ini, digunakan dataset nz_height untuk melakukan interpolasi ketinggian permukaan tanah di wilayah Selandia Baru (New Zealand). Area interpolasi ditentukan berdasarkan cakupan wilayah pengamatan yang dibentuk menggunakan convex hull dari seluruh titik lokasi pengukuran ketinggian.

Grid prediksi dibuat dalam bentuk raster dengan resolusi 100 × 100 sel. Setiap titik pusat sel (cell center) digunakan sebagai lokasi untuk memperkirakan ketinggian berdasarkan hasil pengamatan pada titik sampel.

Langkah-langkah: 1. Membuat grid raster berukuran 100 × 100 sel yang mencakup area pengamatan. 2. Mengambil koordinat X dan Y dari setiap titik pusat sel pada grid. 3. Mengubah koordinat grid menjadi objek sf agar dapat digunakan dalam proses prediksi interpolasi spasial. 4. Menyaring titik-titik prediksi agar hanya mencakup lokasi yang berada di dalam batas wilayah pengamatan. 5. Menampilkan peta sebaran lokasi prediksi.

# Membuat grid raster 100 x 100 sel
grid <- terra::rast(
  terra::vect(map),
  nrows = 100,
  ncols = 100
)

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

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

# Menyaring lokasi prediksi di dalam batas wilayah
coop <- sf::st_filter(coop, map)

# Menampilkan peta lokasi prediksi
tmap::qtm(coop)

Metode Deterministik

Closest Observation

Metode ini menggunakan diagram Voronoi untuk membagi ruang menjadi poligon-poligon yang masing-masing memiliki satu titik pembangkit (generator point). Setiap lokasi di dalam suatu poligon memiliki jarak lebih dekat ke titik pembangkitnya dibandingkan dengan titik pengamatan lainnya.

Pada dataset nz_height, prediksi ketinggian dilakukan berdasarkan nilai elevasi dari titik pengamatan terdekat. Dengan demikian, setiap poligon Voronoi akan memiliki nilai prediksi elevasi yang sama dengan titik pembangkitnya.

Langkah-langkah: 1. Membuat diagram Voronoi berdasarkan lokasi titik pengamatan ketinggian. 2. Membatasi diagram Voronoi menggunakan batas wilayah pengamatan (convex hull). 3. Menggunakan st_intersection() untuk menentukan poligon Voronoi yang mencakup setiap titik lokasi prediksi. 4. Memberikan nilai elevasi berdasarkan titik pengamatan terdekat.

# Membuat diagram Voronoi dari lokasi pengamatan
v <- dismo::voronoi(
  x = terra::vect(d),
  bnd = terra::vect(map)
)

# Menampilkan diagram Voronoi
plot(v)

# Menampilkan titik pengamatan di atas diagram Voronoi
points(
  terra::vect(d),
  cex = 0.5,
  pch = 16
)

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

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

Inverse Distance Weighting (IDW)

Metode IDW menggunakan bobot berdasarkan jarak antara lokasi yang belum teramati dan lokasi yang sudah teramati. Semakin dekat jaraknya, semakin besar bobot yang diberikan pada nilai lokasi teramati.

Langkah-langkah:

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

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

# Melakukan prediksi pada seluruh lokasi dalam grid
resp <- predict(res, coop)
## [inverse distance weighted interpolation]
# Menambahkan koordinat X dan Y
resp$x <- sf::st_coordinates(resp)[, 1]
resp$y <- sf::st_coordinates(resp)[, 2]

# Menyimpan nilai hasil prediksi
resp$pred <- resp$var1.pred

# Membuat peta hasil interpolasi IDW
c2 <- tmap::tm_shape(resp) +
  tmap::tm_dots(
    col = "pred",
    palette = "viridis",
    size = 0.08,
    title = "Elevation (m)"
  ) +
  tmap::tm_layout(
    title = "Hasil Interpolasi IDW - nz_height",
    legend.outside = TRUE
  )

# Menampilkan peta
c2

Nearest Neighbors

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

Langkah-langkah:

  1. Menggunakan gstat() dengan parameter nmax yang menentukan jumlah tetangga dan idp = 0 (bobot sama).
  2. Menggunakan predict() untuk memprediksi harga di lokasi yang tidak teramati.
# Mengubah hasil prediksi menjadi SpatVector
resp_vect <- terra::vect(resp)

# Menggunakan grid awal 100 x 100
# untuk membentuk raster prediksi
pred <- terra::rasterize(
  resp_vect,
  grid,
  field = "pred"
)

# Membatasi raster berdasarkan convex hull
pred <- terra::mask(
  pred,
  terra::vect(map)
)

# Memeriksa hasil raster
print(pred)
## class       : SpatRaster
## size        : 100, 100, 1  (nrow, ncol, nlyr)
## resolution  : 6183.496, 6021.83  (x, y)
## extent      : 1204143, 1822492, 5048309, 5650492  (xmin, xmax, ymin, ymax)
## coord. ref. : NZGD2000 / New Zealand Transverse Mercator 2000 (EPSG:2193)
## source(s)   : memory
## name        :        last
## min value   : 2772.852158
## max value   :   3072.3105
terra::global(pred, "range", na.rm = TRUE)
##           min      max
## last 2772.852 3072.311
# Memeriksa nama layer raster
names(pred)
## [1] "last"
# Mengubah nama layer menjadi pred
names(pred) <- "pred"

# Membuat peta interpolasi IDW
c3 <- tmap::tm_shape(pred) +
  tmap::tm_raster(
    col = "pred",
    palette = "viridis",
    title = "Elevation (m)"
  ) +
  tmap::tm_layout(
    title = "Interpolasi IDW dengan 5 Tetangga Terdekat",
    legend.outside = TRUE
  )

# Menampilkan peta
c3

Pendekatan Ensembel

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

Langkah-langkah:

  1. Menghitung prediksi harga dengan masing-masing metode.
  2. Menggabungkan prediksi setiap metode menggunakan bobot tertentu.
# Membuat diagram Voronoi
v <- terra::voronoi(
  terra::vect(d),
  bnd = terra::vect(map)
)

v <- sf::st_as_sf(v)

# Mengambil nilai elevation dari poligon Voronoi
# untuk setiap titik grid prediksi
p1 <- sf::st_join(
  coop,
  v["vble"],
  join = sf::st_intersects
)$vble

gs <- gstat::gstat(
  formula = vble ~ 1,
  locations = d,
  nmax = nrow(d),
  set = list(idp = 1)
)

p2 <- predict(gs, coop)$var1.pred
## [inverse distance weighted interpolation]
nn <- gstat::gstat(
  formula = vble ~ 1,
  locations = d,
  nmax = 5,
  set = list(idp = 0)
)

p3 <- predict(nn, coop)$var1.pred
## [inverse distance weighted interpolation]
weights <- c(1/3, 1/3, 1/3)

p4 <- p1 * weights[1] +
      p2 * weights[2] +
      p3 * weights[3]

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

resp_sf <- sf::st_as_sf(
  resp,
  coords = c("x", "y"),
  crs = sf::st_crs(coop)
)

resp_vect <- terra::vect(resp_sf)

pred <- terra::rasterize(
  resp_vect,
  grid,
  field = "pred"
)

# Membatasi raster berdasarkan convex hull
pred <- terra::mask(
  pred,
  terra::vect(map)
)

# Memberikan nama layer
names(pred) <- "pred"

c4 <- tmap::tm_shape(pred) +
  tmap::tm_raster(
    col = "pred",
    alpha = 0.6,
    palette = "viridis",
    title = "Elevation (m)"
  ) +
  tmap::tm_layout(
    title = "Interpolasi Kombinasi Voronoi dan IDW",
    legend.outside = TRUE
  )

c4

library(sf)
library(sp)
library(terra)
library(gstat)
library(tmap)
# Konversi raster terra menjadi stars
pred_stars <- stars::st_as_stars(pred)

# Buat ulang peta Ensemble
c4 <- tmap::tm_shape(pred_stars) +
  tmap::tm_raster(
    alpha = 0.6,
    palette = "viridis",
    title = "Log Zinc"
  )

# Gabungkan seluruh peta
tmap::tmap_arrange(
  c1, c2, c3, c4,
  ncol = 2
)

Evaluasi Model

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

Langkah-langkah:

  1. Membagi data menjadi \(K\) bagian dengan dismo::kfold().
  2. Menghitung RMSE untuk setiap metode dengan membandingkan nilai prediksi dan nilai sebenarnya.
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)

# Vektor untuk menyimpan RMSE tiap metode
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) {
  # Pisahkan data menjadi data uji dan data latih
  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 (bobot = invers RMSE, sehingga RMSE lebih kecil -> bobot lebih besar)
  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]
# RMSE pada masing-masing dari 5 fold
data.frame(closest.obs = rmse1, IDW = rmse2,
           nearest.neigh = rmse3, ensemble = rmse4)
##   closest.obs      IDW nearest.neigh  ensemble
## 1    132.9607 125.1960      106.1677  97.69485
## 2    185.6355 192.5830      199.2079 170.72730
## 3    125.5661 165.9054      160.3733 130.45587
## 4    233.7472 171.5258      180.4246 169.90364
## 5    177.9350 204.8474      185.1128 171.64743
# Rata-rata RMSE dari 5 fold
data.frame(closest.obs = mean(rmse1), IDW = mean(rmse2),
           nearest.neigh = mean(rmse3), ensemble = mean(rmse4))
##   closest.obs      IDW nearest.neigh ensemble
## 1    171.1689 172.0115      166.2573 148.0858

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, di mana nilai yang diprediksi adalah variabel \(z\) yang tidak diketahui dan dianggap sebagai variabel acak.

Variogram. Dalam Ordinary Kriging, penting untuk membangun variogram. Variogram adalah grafik yang menggambarkan bagaimana variabilitas antar titik pengamatan berubah seiring jarak antar titik tersebut. Kriging menggunakan informasi variogram ini untuk menghitung kovarians 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.

library(sp)
library(spData)
library(gstat)
library(sf)
library(mapview)

# Memuat dataset nz_height
data("nz_height", package = "spData")

# Melihat struktur data
head(nz_height)
## Simple feature collection with 6 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1204143 ymin: 5048309 xmax: 1389460 ymax: 5168749
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
##   t50_fid elevation                geometry
## 1 2353944      2723 POINT (1204143 5049971)
## 2 2354404      2820 POINT (1234725 5048309)
## 3 2354405      2830 POINT (1235915 5048745)
## 4 2369113      3033 POINT (1259702 5076570)
## 5 2362630      2749 POINT (1378170 5158491)
## 6 2362814      2822 POINT (1389460 5168749)
str(nz_height)
## Classes 'sf' and 'data.frame':   101 obs. of  3 variables:
##  $ t50_fid  : int  2353944 2354404 2354405 2369113 2362630 2362814 2362817 2363991 2363993 2363994 ...
##  $ elevation: int  2723 2820 2830 3033 2749 2822 2778 3004 3114 2882 ...
##  $ geometry :sfc_POINT of length 101; first list element:  'XY' num  1204143 5049971
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
##   ..- attr(*, "names")= chr [1:2] "t50_fid" "elevation"
# Mengubah menjadi objek sf
d <- sf::st_as_sf(nz_height)

# Memeriksa CRS
sf::st_crs(d)
## Coordinate Reference System:
##   User input: EPSG:2193 
##   wkt:
## PROJCRS["NZGD2000 / New Zealand Transverse Mercator 2000",
##     BASEGEOGCRS["NZGD2000",
##         DATUM["New Zealand Geodetic Datum 2000",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4167]],
##     CONVERSION["New Zealand Transverse Mercator 2000",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",173,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",1600000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",10000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["New Zealand - North Island, South Island, Stewart Island - onshore."],
##         BBOX[-47.33,166.37,-34.1,178.63]],
##     ID["EPSG",2193]]
# Membentuk variabel target interpolasi
d$vble <- d$elevation

# Peta sebaran ketinggian
mapview(
  d,
  zcol = "elevation",
  map.types = "CartoDB.Voyager"
)
# Peta sebaran variabel interpolasi
mapview(
  d,
  zcol = "vble",
  map.types = "CartoDB.Voyager"
)

Keterangan komponen semivariogram:

  • Nugget: nilai semivarian ketika jarak antarobservasi mendekati nol, yang merepresentasikan variasi mikrospasial dan kesalahan pengukuran.
  • Sill: nilai maksimum (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: lokasi berdekatan memiliki nilai yang mirip.
  • Semivarian meningkat seiring jarak: ada struktur ketergantungan spasial.
  • Semivarian mencapai sill: variasi spasial relatif stabil setelah jarak tertentu.
# Menghitung variogram empiris
v <- gstat::variogram(
  vble ~ 1,
  data = d,
  cutoff = 1000,
  width = 50
)

# Visualisasi variogram empiris
plot(
  v,
  main = "Experimental Variogram - nz_height Dataset",
  cex = 1.2
)

# Menghitung variogram eksperimental
v <- gstat::variogram(
  vble ~ 1,
  data = d
)

# Menampilkan variogram standar
plot(
  v,
  main = "Experimental Variogram - nz_height Dataset"
)

# Daftar model variogram yang tersedia
gstat::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 grafik model variogram yang tersedia
gstat::show.vgms(
  par.strip.text = list(cex = 0.75)
)

# Model variogram awal
vinitial <- gstat::vgm(
  psill = 0.5,
  model = "Sph",
  range = 900,
  nugget = 0.1
)

# Visualisasi variogram empiris dengan model awal
plot(
  v,
  vinitial,
  cutoff = 1000,
  cex = 1.5,
  main = "Empirical Variogram dan Initial Spherical Model"
)

# Fitting model Spherical (Sph) ke variogram empiris
fv <- gstat::fit.variogram(
  object = v,
  model = gstat::vgm(
    psill = 0.5,
    model = "Sph",
    range = 900,
    nugget = 0.1
  )
)

# Menampilkan parameter hasil fitting
fv
##   model psill range
## 1   Nug   0.1     0
## 2   Sph   0.5   900
# Visualisasi variogram empiris dan model hasil fitting
plot(
  v,
  fv,
  cex = 1.5,
  main = "Fitting Spherical Variogram - nz_height Dataset"
)

# Memanggil package
library(automap)

# Autofit variogram untuk data elevation
afvg <- automap::autofitVariogram(
  vble ~ 1,
  input_data = d
)

# Menampilkan hasil fitting variogram
afvg
## $exp_var
##      np       dist     gamma dir.hor dir.ver   id
## 1   960   3230.345 50696.157       0       0 var1
## 2  1000   9112.754 47691.628       0       0 var1
## 3   835  15089.592 61582.966       0       0 var1
## 4   520  21243.253 64476.302       0       0 var1
## 5   196  31173.856 34053.000       0       0 var1
## 6    81  38839.896 64411.877       0       0 var1
## 7   222  64225.297 52425.939       0       0 var1
## 8    36  80820.969 88380.153       0       0 var1
## 9    44 143699.900 37604.193       0       0 var1
## 10  181 179121.151 40408.367       0       0 var1
## 11  118 209530.666 45071.953       0       0 var1
## 12   29 272846.318  4788.345       0       0 var1
## 
## $var_model
##   model    psill    range
## 1   Nug 50224.21     0.00
## 2   Gau 18055.79 29268.79
## 
## $sserr
## [1] 647.6276
## 
## attr(,"class")
## [1] "autofitVariogram" "list"
# Plot hasil autofit variogram
plot(afvg)

library(gstat)
library(ggplot2)
library(viridis)
library(sf)

# 1. Membuat model Ordinary Kriging
k <- gstat::gstat(
  formula = vble ~ 1,
  locations = d,
  model = fv,
  nmax = nrow(d)
)

# 2. Prediksi pada seluruh grid
kpred <- predict(k, newdata = coop)
## [using ordinary kriging]
# 3. Visualisasi prediksi elevation
p1 <- ggplot() +
  geom_sf(
    data = kpred,
    aes(color = var1.pred),
    size = 1
  ) +
  geom_sf(
    data = d,
    color = "black",
    alpha = 0.15,
    size = 0.1
  ) +
  scale_color_viridis(name = "Elevation (m)") +
  labs(
    title = "Ordinary Kriging Prediction",
    subtitle = "Prediksi Ketinggian Permukaan Tanah"
  ) +
  theme_bw()

print(p1)

# 4. Visualisasi varians prediksi
p2 <- ggplot() +
  geom_sf(
    data = kpred,
    aes(color = var1.var),
    size = 1
  ) +
  scale_color_viridis(name = "Prediction Variance") +
  labs(
    title = "Ordinary Kriging Prediction Variance",
    subtitle = "Varians Ketidakpastian Prediksi"
  ) +
  theme_bw()

print(p2)

library(sf)
library(gstat)

# Menggunakan seluruh observasi
d_sample <- d

# Leave-One-Out Cross Validation
cv_kriging <- gstat::krige.cv(
  formula = vble ~ 1,
  locations = d_sample,
  model = fv,
  nmax = nrow(d)
)

# Menghitung RMSE
RMSE_Kriging <- sqrt(
  mean(cv_kriging$residual^2, na.rm = TRUE)
)

# Menghitung MAE
MAE_Kriging <- mean(
  abs(cv_kriging$residual),
  na.rm = TRUE
)

# Menampilkan hasil evaluasi
cat("Jumlah observasi:", nrow(d_sample), "\n")
## Jumlah observasi: 101
cat("RMSE Ordinary Kriging:", RMSE_Kriging, "\n")
## RMSE Ordinary Kriging: 180.226
cat("MAE Ordinary Kriging :", MAE_Kriging, "\n")
## MAE Ordinary Kriging : 133.119

Kesimpulan

Metode dengan rata-rata RMSE terkecil pada tabel di atas merupakan metode dengan kinerja prediksi terbaik untuk data nz_height. Berdasarkan perbandingan RMSE, metode Ensemble menghasilkan kesalahan prediksi paling kecil sebesar 148,0858, diikuti oleh IDW sebesar 166,2573, Closest Observation sebesar 171,1689, dan Nearest Neighbor sebesar 172,0115. Sementara itu, Ordinary Kriging menghasilkan RMSE sebesar 180,226. Dengan demikian, metode Ensemble memiliki kinerja prediksi terbaik dibandingkan metode lainnya karena menghasilkan nilai RMSE paling rendah.