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 coyDataset 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 326Inisialisasi 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(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')
Perhatikan bahwa bentuk semi-variogram eksperimental di atas, memiliki bentuk yang tidak terlalu baik. Dengan demikian, model linear dapat dikatakan lebih cocok terhadap masalah semi-variogram di atas. Perhatikan pula pada model fitted semivariogram di atas, nilai psill, range, dan nugget effect diperoleh secara otomatis melalui fungsi fit.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('Lin'))## Warning in fit.variogram(variogram_eksperimental, model = vgm("Lin")):
## singular model in variogram fitplot(variogram_eksperimental, model = fitted_model, main = 'Variogram Eksperimental dan Fitted Model')
new <- data.frame(x = 0.93, y = -0.9)
coordinates(new) <- ~x+y
prediction_result <- krige(DZ ~ 1, new_coord_data, new, model = fitted_model)## [using ordinary kriging]print(prediction_result)## coordinates var1.pred var1.var
## 1 (0.93, -0.9) 297.7024 0.1783236Pertanyaan selanjutnya, apakah hasil prediksi di atas masuk akal? Untuk melihat apakah hasil prediksi di atas cukup masuk akal, akan dilakukan visualiasi scatter plot dengan pemetaan z pada besar point.
coord_data_df <- dataset %>% select(x,y,DZ)
hasil_prediksi <- data.frame(x = 0.93, y = -0.9, DZ = 297.7024, var = 0.17832)
ggplot(coord_data_df, aes(x = x, y = y, size = DZ)) + geom_point() +
geom_point(data = hasil_prediksi, aes(x = x, y = y, size = DZ, colour = 'red')) +
labs(title = 'Scatter Plot dengan Pemetaan DZ pada Size Titik',
x = 'X',
y = 'Y') +
t