Soal:

Chemical Oxygen Demand (COD), yaitu ukuran jumlah oksigen yang dibutuhkan untuk mengurai polutan organik dalam air secara kimiawi menggunakan oksidan kuat. Satuan yang digunakan untuk menyatakan nilai COD adalah miligram per liter (mg/L). Lakukan analisis dengan ordinary kriging untuk mengestimasi kandungan COD pada lokasi ke-7.

NO LOKASI GARIS LINTANG GARIS BUJUR COD
1 Kali Surabaya di Intake Jagir 07°03’762” 112°34’576” 14,058
2 Kali Surabaya di Dam Gunungsari 06°59’137” 112°22’434” 13,994
3 Kali Surabaya di intake Karangpilang 07°06’638” 112°09’247” 12,888
4 Kali Surabaya sesudah outlet PT. Suparma 07°03’300” 112°07’108” 12,107
5 Kali Surabaya sesudah pertemuan dengan Kali Tengah 07°09’666” 112°7’108” 12,086
6 Kali Surabaya sesudah outlet PT. Miwon 07°08’812” 111°35’930” 11,873
7 Kali Surabaya sebelum outlet PT. SAK 07°21’211” 111°08’932” ?
8 Kali Surabaya di DAM Mirip 07°23’457” 111°26’509” 12,201

Penyelesian:

Library

library(sp)
library(gstat)

Memanggil Data

data=read.table(file.choose(),header=TRUE)
data
##   Latitude Longitude    COD
## 1   7.2617  112.7267 14.058
## 2   7.0214  112.4872 13.994
## 3   7.2772  112.2186 12.888
## 4   7.1667  112.1467 12.107
## 5   7.3350  112.1467 12.086
## 6   7.3589  111.8417 11.873
## 7   7.5103  111.5747 12.201

1. Analisis Statistika Deskriptif

summary(data$COD)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   11.87   12.10   12.20   12.74   13.44   14.06
var(data$COD)
## [1] 0.8670091

Tabel 1. Statistika Deskriptif Data COD

Minimum Rata-rata Median Variansi Maksimum
11,87 12,74 12,20 0,8670091 14,06

Interpretasi:

Berdasarkan hasil Analisis Statistika Deskriptif Data Chemical Oxygen Demand (COD), diketahui bahwa nilai minimum sebesar 11,87 dan nilai maksimum sebesar 14,06 dengan nilai rata-rata sebesar 12,74 serta nilai variansi sebesar 0,8670091.

2. Jarak Euclid

data1=data.frame(data$Longitude,data$Latitude)
data1
##   data.Longitude data.Latitude
## 1       112.7267        7.2617
## 2       112.4872        7.0214
## 3       112.2186        7.2772
## 4       112.1467        7.1667
## 5       112.1467        7.3350
## 6       111.8417        7.3589
## 7       111.5747        7.5103
a=as.matrix(data1)
a
##      data.Longitude data.Latitude
## [1,]       112.7267        7.2617
## [2,]       112.4872        7.0214
## [3,]       112.2186        7.2772
## [4,]       112.1467        7.1667
## [5,]       112.1467        7.3350
## [6,]       111.8417        7.3589
## [7,]       111.5747        7.5103
b=as.matrix(dist(a))
b
##           1         2         3         4         5         6         7
## 1 0.0000000 0.3392703 0.5083364 0.5877287 0.5846135 0.8903218 1.1785185
## 2 0.3392703 0.0000000 0.3709172 0.3702058 0.4629095 0.7284068 1.0352195
## 3 0.5083364 0.3709172 0.0000000 0.1318327 0.0922521 0.3856533 0.6847940
## 4 0.5877287 0.3702058 0.1318327 0.0000000 0.1683000 0.3605078 0.6672668
## 5 0.5846135 0.4629095 0.0922521 0.1683000 0.0000000 0.3059350 0.5982592
## 6 0.8903218 0.7284068 0.3856533 0.3605078 0.3059350 0.0000000 0.3069380
## 7 1.1785185 1.0352195 0.6847940 0.6672668 0.5982592 0.3069380 0.0000000
n=length(a[,1])
n
## [1] 7

Tabel 2. Jarak Euclid

j 1 2 3 4 5 6 7
1 0 0,3392703 0,5083364 0,5877287 0,5846135 0,8903218 1,1785185
2 0,3392703 0 0,3709172 0,3702058 0,4629095 0,7284068 1,0352195
3 0,5083364 0,3709172 0 0,1318327 0,0922521 0,3856533 0,6847940
4 0,5877287 0,3702058 0,1318327 0 0,1683000 0,3605078 0,6672668
5 0,5846135 0,4629095 0,0922521 0,1683000 0 0,3059350 0,5982592
6 0,8903218 0,7284068 0,3856533 0,3605078 0,3059350 0 0,3069380
7 1,1785185 1,0352195 0,6847940 0,6672668 0,5982592 0,3069380 0

3. Kestasioneran Data

plot(data$COD)

Interpretasi:

Berdasarkan hasil plot di atas, diketahui bahwa data kandungan Chemical Oxygen Demand (COD) menunjukkan pada tiap titik lokasi tidak memiliki kecenderungan pola data tren tertentu, sehingga dapat dikatakan bahwa data kandungan COD adalah stasioner.

4. Menghitung nilai semivariogram eksperimental

coordinates(data1)=~data.Latitude+data.Longitude
coordinates(data1)
##   data.Latitude data.Longitude
## 1        7.2617       112.7267
## 2        7.0214       112.4872
## 3        7.2772       112.2186
## 4        7.1667       112.1467
## 5        7.3350       112.1467
## 6        7.3589       111.8417
## 7        7.5103       111.5747
maxdist=max(b)
maxdist
## [1] 1.178519
kelas=round(1+3.3*log10(n))
kelas
## [1] 4
interval=maxdist/kelas
interval
## [1] 0.2946296
varClassic=variogram(data$COD~1,data1,cutoff=maxdist,width=interval)
varClassic
##   np      dist     gamma dir.hor dir.ver   id
## 1  3 0.1307949 0.2089343       0       0 var1
## 2 11 0.4166378 0.8513902       0       0 var1
## 3  4 0.6696817 0.6240839       0       0 var1
## 4  3 1.0346866 1.9062538       0       0 var1

Tabel 3. Nilai Semivariogram Eksperimental

Kelas Interval Banyak
Pasangan 		Rata-rata
Jarak 		\(\gamma(h)\)
1 \(0 \le M < 0,2946\) 3 0,1307949 0,2089343
2 \(0,2946 \le M < 0,5892\) 11 0,4166378 0,8513902
3 \(0,5892 \le M < 0,8838\) 4 0,6696817 0,6240839
4 \(0,8838 \le M < 1,1784\) 3 1,0346866 1,9062538 
plot(varClassic)

r=0.4419

\[Sill(C) = 0,8670091 \text{ (Nilai Variansi Data)}\]

\[Range(a) = \frac{0,2946 + 0,5892}{2} = 0,4419\]

Nilai \(Range(a)\) diperoleh dari rata-rata interval yang nilai semivariogramnya paling mendekati nilai \(Sill(C)\), yaitu interval kelas 2.

s=var(data$COD)
dist=varClassic$dist
dist
## [1] 0.1307949 0.4166378 0.6696817 1.0346866
i=dist<=r
i
## [1]  TRUE  TRUE FALSE FALSE
j=dist>r
j
## [1] FALSE FALSE  TRUE  TRUE

5. SPHERICAL

sphe1=s*((1.5*(dist[i]/r)-(0.5*(dist[i]^3)/(r^3))))
sphe1
## [1] 0.3736894 0.8628399
sphe2=s+dist[j]*0
sphe2
## [1] 0.8670091 0.8670091
spherical=c(sphe1,sphe2)
spherical
## [1] 0.3736894 0.8628399 0.8670091 0.8670091

Model Spherical

\[\hat{\gamma}(h) = \begin{cases} 0,8670091 \left[ 1,5 \frac{f(h)}{0,4419} - 0,5 \left( \frac{f(h)}{0,4419} \right)^3 \right], \text{untuk } f(h) \le a \\ \quad\quad\quad 0,8670091, \text{untuk } f(h) > a \end{cases}\]

6. EXPONENSIAL

eksponensial=s*(1-(exp(-(3*dist)/r)))
eksponensial
## [1] 0.5102360 0.8157676 0.8578140 0.8662376

Model Eksponensial

\[\hat{\gamma}(h) = 0,8670091 \left[ 1 - e^{\frac{3 f(h)}{0,4419}} \right]\]

7. GAUSSIAN

gaussian=s*(1-exp(-3*dist^2/r^2))
gaussian
## [1] 0.2003816 0.8067744 0.8661265 0.8670091

Model Gaussian

\[\hat{\gamma}(h) = 0,8670091 \left[ 1 - e^{\frac{3 f(h)^2}{0,4419^2}}\right]\]

semivariogramteoritis=data.frame(varClassic$np,varClassic$dist,varClassic$gamma,spherical,eksponensial,gaussian)
semivariogramteoritis
##   varClassic.np varClassic.dist varClassic.gamma spherical eksponensial
## 1             3       0.1307949        0.2089343 0.3736894    0.5102360
## 2            11       0.4166378        0.8513902 0.8628399    0.8157676
## 3             4       0.6696817        0.6240839 0.8670091    0.8578140
## 4             3       1.0346866        1.9062538 0.8670091    0.8662376
##    gaussian
## 1 0.2003816
## 2 0.8067744
## 3 0.8661265
## 4 0.8670091

Tabel 5. Nilai Semivariogram Teoritis

Kelas Banyak
Pasangan 		Rata-rata
Jarak 		\(\mathbf{\gamma(h)}\) Spherical Eksponensial Gaussian
1 3 0,1307949 0,2089343 0,3736894 0,5102360 0,2003816
2 11 0.4166378 0.8513902 0.8628399 0.8157676 0.8067744
3 4 0.6696817 0.6240839 0.8670091 0.8578140 0.8661265
4 3 1.0346866 1.9062538 0.8670091 0.8662376 0.8670091

8. Menghitung MSE

n=length(varClassic$dist)
n
## [1] 4
k=c(1:n)
k
## [1] 1 2 3 4

a. Model Spherical

Selanjutnya, menghitung nilai MSE untuk setiap model dan diperoleh hasil sebagai berikut:

Model Spherical

\[MSE = \frac{\sum_{i=1}^n \left( \gamma(h_i) - \gamma^*(h_i) \right)^2}{n^*}\]

\[_{MSE} = \frac{(0,2089343 - 0,3736894)^2 + (0,8513902 - 0,8628399)^2 + (0,6240839 - 0,8670091)^2 + (1,9062538 - 0,8670091)^2}{4}\]

\[MSE = 0,2915794\]

MSE_Spherical=sum(((varClassic$gamma[k]-spherical[k])^2)/n)
MSE_Spherical
## [1] 0.2915794

b. Model Gaussian

Model Gaussian

\[MSE = \frac{\sum_{i=1}^n \left( \gamma(h_i) - \gamma^*(h_i) \right)^2}{n^*}\]

\[_{MSE} = \frac{(0,2089343 - 0,2003816)^2 + (0,8513902 - 0,8067744)^2 + (0,6240839 - 0,8661265)^2 + (1,9062538 - 0,8670091)^2}{4}\]

\[MSE = 0,2851695\]

MSE_Gaussian=sum(((varClassic$gamma[k]-gaussian[k])^2)/n)
MSE_Gaussian
## [1] 0.2851695

c. Model Eksponensial

Model Eksponensial

\[MSE = \frac{\sum_{i=1}^n \left( \gamma(h_i) - \gamma^*(h_i) \right)^2}{n^*}\]

\[_{MSE} = \frac{(0,2089343 - 0,5102360)^2 + (0,8513902 - 0,8157676)^2 + (0,6240839 - 0,8578140)^2 + (1,9062538 - 0,8662376)^2}{4}\]

\[MSE = 0,3070788\]

MSE_Eksponensial=sum(((varClassic$gamma[k]-eksponensial[k])^2)/n)
MSE_Eksponensial
## [1] 0.3070788
MSE2=data.frame(MSE_Spherical,MSE_Gaussian,MSE_Eksponensial)
MSE2
##   MSE_Spherical MSE_Gaussian MSE_Eksponensial
## 1     0.2915794    0.2851695        0.3070788

Tabel 6. Nilai MSE Semivariogram Teoritis

Spherical Eksponensial Gaussian
0,2915794 0,3070788 0,2851695

Interpretasi:

Berdasarkan hasil dari MSE diperoleh hasil perbandingan untuk setiap model semivariogram teoritis. Nilai MSE terkecil didapatkan pada model Gaussian sebesar 0,2851695, sehingga model Gaussian adalah model semivariogram teoritis terbaik.

9. Menghitung Nilai Estimasi

longitude = c(111.3922)
latitude = c(7.4086)
m=data.frame(longitude,latitude)
new_spatial_data <- SpatialPoints(m)
spatial_data=SpatialPointsDataFrame(data, data =data)
Koordinat Estimasi
(111,3922; 7,4086) 12,87317

Interpretasi:

Berdasarkan tabel di atas, hasil estimasi kandungan Chemical Oxygen Demand (COD), pada lokasi ke-7 dengan titik koordinat (111,3922; 7,4086) adalah sebesar 12,87317.

10. Model Terbaik

variogram=vgm(s,model="Gau",r)
variogram
##   model     psill  range
## 1   Gau 0.8670091 0.4419
OK=krige(data$COD~1,spatial_data,new_spatial_data,model=variogram)
## [using ordinary kriging]
print(OK)
##          coordinates var1.pred var1.var
## 1 (111.3922, 7.4086)  12.87317 1.074214