Analisis Semivariogram

Praktikum Analisis Data, Kamis 09.00 - 11.00

Load package yang akan digunakan pada praktikum kali ini. Instalasi package dapat dilakukan melalui install.packages(‘nama_package’). Dengan demikian, untuk menginstall package di bawah ini dapat menggunakan code berikut:
install.packages(‘ggplot2’)
install.packages(‘akima’)
install.packages(‘dplyr’)
install.packages(‘tidyr’)
install.packages(‘gstat’)
install.packages(‘sp’)
Untuk praktikum kali ini, tidak akan dijelaskan secara mendalam mengenai kegunaan dari setiap package. Peserta praktikum dapat langsung menggunakan code di bawah ini, termasuk untuk visualisasi data.

library(ggplot2) #Package untuk melakukan visualisasi dengan grammer of graphics
library(akima) #Manipulasi data
library(tidyr) #Manipulasi data
library(dplyr) #Manipulasi data
library(gstat) #Package utama untuk melakukan analisis spasial
library(sp) #Spasial coy

Dataset Minyak Jatibarang
Perhatikan bahwa terapat 4 buah kolom/features pada data set di bawah ini yakni (1) koordinat x, (2) koordinat y, (3) level (K-FRACTURE), (4) DZ atau kedalaman sumur.
Dataset di bawah ini akan digunakan untuk kepentingan simulasi awal analisis semivariogram.

dataset <- data.frame(x = c(0.6309, 0.9241, 0.5677, 1.0561, 0.1290, 0.9848, 1.2489, 1.2200, 0.5042, 0.5174, 0.4646, 0.8448), 
                      y = c(-1.3109, -1.2761, -1.0127, -0.8133, -0.5951, -0.3901, -0.4034, -0.6375, -1.5291, -1.2708, -0.6614, -0.9009), 
                      level = c(35.445, 26.645, 51.72, 14.331, 13.333, 4.683, 49.754, 35.832, 30.467, 26.741, 16.776, 40.154), 
                      DZ = c(213, 291, 291, 388, 138, 108, 168, 256,238, 326, 181, 359))
head(dataset, n = 10)

##         x       y  level  DZ
## 1  0.6309 -1.3109 35.445 213
## 2  0.9241 -1.2761 26.645 291
## 3  0.5677 -1.0127 51.720 291
## 4  1.0561 -0.8133 14.331 388
## 5  0.1290 -0.5951 13.333 138
## 6  0.9848 -0.3901  4.683 108
## 7  1.2489 -0.4034 49.754 168
## 8  1.2200 -0.6375 35.832 256
## 9  0.5042 -1.5291 30.467 238
## 10 0.5174 -1.2708 26.741 326

Inisialisasi Fungsi
Berikut ini adalah fungsi yang digunakan untuk kepentingan manipulasi data sebelum dilakukannya visualisasi. grid_make digunakan untuk menghasilkan grid matriks lokasi K-Fracture. interpolation_func digunakan untuk menghasilkan grid matriks lokasi K-Fracture dengan interpolasi data antar titiknya. wo_interpolation pada intinya sama dengan grid_make namun, bentuk grid dibuat berbeda karena menggunakan fungsi visualisasi yang berbeda. t adalah format theme untuk kepentingan visualiasi dengan ggplot2.

grid_make <- function(df){
  if(sum((c('x','y','level') %in% names(df)) == TRUE) < 3){
    stop("Data frame harus mengandung kolom bernama 'x', 'y', dan 'level'", call. = FALSE)
  }
  x_sorted <- sort(df$x)
  y_sorted <- sort(df$y)
  value_mat <- matrix(0, ncol = length(x_sorted), nrow = length(y_sorted))
  for (a in df$x){
    for(b in df$y){
      val <- df %>% filter(x == a & y == b)
      if(nrow(val) == 0){
        value_mat[which(a == x_sorted), which(b == y_sorted)] <- 0}
      else{
        value_mat[which(a == x_sorted),which(b == y_sorted)] <- val$level
      }
    }
  }
  return(list(grid = value_mat, x_sorted = x_sorted, y_sorted = y_sorted))
}

interpolation_func <- function(df){
  if(sum((c('x','y','level') %in% names(df)) == TRUE) < 3){
    stop("Data frame harus mengandung kolom bernama 'x', 'y', dan 'level'", call. = FALSE)
  }
  the_grid <- interp(df$x, df$y, df$level)
  the_grid1 <- expand.grid(x = the_grid$x, y = the_grid$y)
  the_grid1$z <- as.vector(the_grid$z)
  return(as.data.frame(the_grid1) %>% mutate(lvl = ifelse(is.na(z), 0, z)))
}

wo_interpolation <- function(df){
  x_sorted <- sort(df$x)
  y_sorted <- sort(df$y)
  value_mat <- matrix(0, ncol = length(x_sorted), nrow = length(y_sorted))
  for (a in df$x){
    for(b in df$y){
      val <- df %>% filter(x == a & y == b)
      if(nrow(val) == 0){
        value_mat[which(a == x_sorted), which(b == y_sorted)] <- 0}
      else{
        value_mat[which(a == x_sorted),which(b == y_sorted)] <- val$level
      }
    }
  }
  thegrid <- expand.grid(x = x_sorted, y = y_sorted)
  thegrid$z <- as.vector(value_mat)
  return(as.data.frame(thegrid))
}
t <- theme(panel.background = element_blank(),
           panel.grid.major = element_line(colour = 'grey', linetype = 2),
           plot.title = element_text(hjust = 0.5),
           plot.subtitle = element_text(hjust = 0.5))

Implementasi Fungsi

raw_grid <- grid_make(dataset)
int_grid <- interpolation_func(dataset)
wo_int_grid <- wo_interpolation(dataset)

Visualisasi 3D Plot

persp(raw_grid$x_sorted, raw_grid$y_sorted, raw_grid$grid, ticktype = 'detailed', 
      shade = 0.2, main = '3D Plot of Raw Data', theta = 30, xlab = 'X', ylab = 'Y', zlab = 'K-FRACTURE')

Visualisasi Peta Kontur dengan Interpolasi

ggplot(int_grid, aes(x = x, y = y, z = lvl)) + 
  stat_contour(bins = 20, geom = 'polygon', aes(fill = ..level..)) +
  scale_fill_gradient(low = 'yellow', high = 'red') +
  labs(title = 'Contour Map Data K-FRACTURE Berdasarkan Koordinat X dan Y',
       subtitle = 'Note: dengan interpolasi',
       fill = 'Level',
       x = 'X',
       y = 'Y') +
  t

Visualisasi Peta Kontur Tanpa Interpolasi

ggplot(wo_int_grid, aes(x = x, y = y, z = z)) + 
  stat_contour(bins = 20, geom = 'polygon', aes(fill = ..level..)) +
  scale_fill_gradient(low = 'yellow', high = 'red') +
  labs(title = 'Contour Map Data K-FRACTURE Berdasarkan Koordinat X dan Y',
       subtitle = 'Note: tanpa interpolasi',
       fill = 'Level',
       x = 'X',
       y = 'Y') +
  t

Perhatikan bahwa software Surfer dapat menghasilkan peta kontur yang lebih baik. Hal ini dikarenakan software Surfer memilili fungsi interpolasi yang lebih baik (berbayar gitu lhoo…).

Analisis Semivariogram

Kali ini, akan dibahas mengenai proses pembuatan model Semivariogram menggunakan package gstat. Akan dilakukan analisis semivariogram terhadap kedalaman sumur sehingga kali ini akan digunakan kolom DZ. Pada kode di bawah ini, perhatikan bahwa variogram adalah fungsi yang digunaan untuk menghitung nilai variogram eksperimental. vgm adalah fungsi untuk mendeklarasikan model. fit.variogram adalah fungsi untuk fitting model variogram.

new_coord_data <- dataset %>% select(x,y,DZ) #dapat juga menggunakan dataset[,c('x','y','DZ')]
coordinates(new_coord_data) <- ~x+y
variogram_eksperimental <- variogram(log(DZ)~1, data = new_coord_data)
fitted_model <- fit.variogram(variogram_eksperimental, model = vgm('Exp'))

## Warning in fit.variogram(variogram_eksperimental, model = vgm("Exp")): No
## convergence after 200 iterations: try different initial values?

plot(variogram_eksperimental, model = fitted_model, main = 'Variogram Eksperimental dan Fitted Model')