1 Trend surface interpolation
Data yang digunakan adalah data topo yang mencatat spatial topographic. Namun sebelumnya, kita perlu memanggil beberapa package terlbih dulu.
library(MASS)
library(raster)
library(gstat)
library(sp)
data(topo)
str(topo)## 'data.frame': 52 obs. of 3 variables:
## $ x: num 0.3 1.4 2.4 3.6 5.7 1.6 2.9 3.4 3.4 4.8 ...
## $ y: num 6.1 6.2 6.1 6.2 6.2 5.2 5.1 5.3 5.7 5.6 ...
## $ z: int 870 793 755 690 800 800 730 728 710 780 ...head(topo)help(topo)Selanjutnya kita perlu mendefinisikan koordinat pada data topo.
coordinates(topo)<-~x+y
par(mfrow=c(1,2))
plot(topo)
plot(topo, pch=1, cex=topo$z/600)
Tahapan selanjutnya adalah mempersiapkan grid untuk interpolasi.
x=seq(0,6.3, by=0.1);y=seq(0,6.3, by=0.1)
topo.grid <- expand.grid(x=x,y=y)
gridded(topo.grid) <- ~x+y
plot(topo.grid)
1.1 1st-Order Trend Surface
Selanjutnya akan dilakukan 1st-order trend surface.
x1 <- gstat(formula=topo$z ~ 1, data=topo, degree=1)
kt <- predict(x1, newdata=topo.grid)## [ordinary or weighted least squares prediction]spplot(kt[1], contour=TRUE,main="topo trend surface (1)-gstat")
1.2 2nd-Order Trend Surface
Berikut ini adalah contoh syntax untuk melakukan 2nd-order trend surface.
x <- gstat(formula=topo$z ~ 1, data=topo, degree=2)
kt2 <- predict(x, newdata=topo.grid)## [ordinary or weighted least squares prediction]spplot( kt2[1], contour=TRUE, main="topo trend surface (2)-gstat")
2 Inverse Distance Weighting (IDW)
Beberapa fitur dari pendekatan ini adalah:
titik yang lebih dekat kepada titik yang akan diinterpolasi akan lebih berpengaruh terhadap hasil interpolasi di titik tersebut
teknik yang digunakan untuk memprediksi suatu titik adalah dengan memberikan bobot dari pengaruh titik di sekitarnya berdasarkan jarak terhadap titik yang akan diinterpolasi
2.1 IDW on topo data set
kid <- krige(topo$z ~ 1, topo, topo.grid, NULL)## [inverse distance weighted interpolation]spplot(kid["var1.pred"], main = "Topo IDW",contour=TRUE)
3 Thiessen Polygons
Data meuse akan digunakan sebagai ilustrasi pendekatan ini.
library(sp)
data(meuse)
names(meuse)## [1] "x" "y" "cadmium" "copper" "lead" "zinc" "elev"
## [8] "dist" "om" "ffreq" "soil" "lime" "landuse" "dist.m"hist(meuse$zinc)
plot(meuse$x,meuse$y,asp=1)
Perhatikan bahwa data meuse masih berupa data.frame.
class(meuse)## [1] "data.frame"Oleh karenanya, kita akan mendefinisikan koordinat lokasi pada data meuse, agar data tersebut menjadi class data spatial. Terlihat pula pada gambar bahwa kandungan zinc cenderung lebih tinggi pada daerah di pinggir sungai (daerah paling kiri) yang ditunjukkan dengan lingkaran yang lebih besar.
coordinates(meuse) = ~x + y # convert data to spatial points data frame
class(meuse)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"plot(meuse,pch=1, cex=meuse$zinc/500)
Selanjutnya, kita perlu mempersiapkan grid untuk interpolasi. Kebetulan, R sudah memepersiapkan data grid tersebut yang berpadanan dengan data meuse.
data(meuse.grid)
str(meuse.grid)## 'data.frame': 3103 obs. of 7 variables:
## $ x : num 181180 181140 181180 181220 181100 ...
## $ y : num 333740 333700 333700 333700 333660 ...
## $ part.a: num 1 1 1 1 1 1 1 1 1 1 ...
## $ part.b: num 0 0 0 0 0 0 0 0 0 0 ...
## $ dist : num 0 0 0.0122 0.0435 0 ...
## $ soil : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ ffreq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...Sama seperti sebelumnya, kita perlu mendefinisikan koordinat lokasi untuk dataa meuse.grid.
coordinates(meuse.grid) = ~x + y # converts to spatial class
gridded(meuse.grid) <- TRUEKemudian setelah tersedia data titik dan data grid, kita dapat melakukan interpolasi pada data meuse.
zn.tp = krige(log(zinc) ~ 1, meuse, meuse.grid, nmax = 1)## [inverse distance weighted interpolation]#for this search, the neighborhood is set to nmax=1
image(zn.tp["var1.pred"])
points(meuse, pch = "+", cex = 0.65)
cc = coordinates(meuse)Plot interpolasi dapat pula ditampilkan dengan package tripack.
library(tripack)
image(zn.tp["var1.pred"])
points(meuse, pch = "+", cex = 0.65)
cc = coordinates(meuse)
plot(voronoi.mosaic(cc[, 1], cc[, 2]), do.points = FALSE, add = TRUE)
4 Spline Interpolation
Pendekatan ini menggunakan fungsi interpolasi pada beberapa bagian yang dipisahkan oleh beberapa titik (knot).
Kita dapat melakukan interpolasi ini dengan menggunakan fungsi pada package akima.
library(akima)
data(meuse)
meuse.zinc <- data.frame(cbind(meuse$x,
meuse$y,meuse$zinc))
meuse.akim <- interp(meuse.zinc$X1, meuse.zinc$X2, meuse.zinc$X3)
image(meuse.akim,,xlab="", ylab="",asp=1)
contour(meuse.akim,add=TRUE, levels = seq(100,2000,200))
5 Choosing the Right Interpolation
The quality of sample point set can affect choice of interpolation method as well. If the sample points are poorly distributed or there are few of them, the surface might not represent the actual terrain very well. If we have too few sample points, add more sample points in areas where the terrain changes abruptly or frequently, then try using Kriging.
The real-world knowledge of the subject matter will initially affect which interpolation method to use. If it is know that some of the features in surface exceed the z value, for example, and that IDW will result in a surface that does not exceed the highest or lowest z value in the sample point set, choose the Spline method.
If we know that the splined surface might end up with features that we know don’t exist because Spline interpolation doesn’t work well with sample points that are close together and have extreme differences in value, you might decide to try IDW.
6 References
Anonymous. 2017. Interpolation in R. https://mgimond.github.io/Spatial/interpolation-in-r.html [17 Okt 2017]
[gisresource]. 2013. Choosing the Right Interpolation Method. https://www.gisresources.com/choosing-the-right-interpolation-method_2/ [17 Okt 2017)
Lewis, J. Spatial Interpolation: Methods for Extending Information on Spatial Data Fields. Lewis, J. 2013. GEOSTAT13 Tutorial: Spatial Interpolation with R.
Sang. 2006. Lecture Notes: Spatial Interpolation. http://sbs.mnsu.edu/geography/people/sang/LECTURE8_SPATIAL_INTERPOLATION.pdf [17 Okt 2017]
Pustaka lain yang relevan.