Library
## Warning: package 'spData' was built under R version 4.5.3
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
## terra 1.8.93
## Warning: package 'tmap' was built under R version 4.5.3
## Warning: package 'viridis' was built under R version 4.5.3
## Loading required package: viridisLite
## Warning: package 'dismo' was built under R version 4.5.3
## Loading required package: raster
## Loading required package: sp
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:terra':
##
## intersect, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Warning: package 'spatstat.data' was built under R version 4.5.3
library(sp)
library(maptiles)
## Warning: package 'maptiles' was built under R version 4.5.3
## Warning: package 'gstat' was built under R version 4.5.3
library(dismo)
library(mapview)
## Warning: package 'mapview' was built under R version 4.5.3
## Warning: package 'automap' was built under R version 4.5.3
library(ggplot2)
library(rlang)
Data Diameter Pohon di Georgia, Amerika Serikat
data(longleaf)
longleaf_df <- data.frame(
x = longleaf$x,
y = longleaf$y,
vble = longleaf$marks
)
longleaf_sf <- st_as_sf(longleaf_df, coords = c("x", "y"), crs = 32617)
longleaf_sf
## Simple feature collection with 584 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 0 ymin: 0 xmax: 200 ymax: 200
## Projected CRS: WGS 84 / UTM zone 17N
## First 10 features:
## vble geometry
## 1 32.9 POINT (200 8.8)
## 2 53.5 POINT (199.3 10)
## 3 68.0 POINT (193.6 22.4)
## 4 17.7 POINT (167.7 35.6)
## 5 36.9 POINT (183.9 45.4)
## 6 51.6 POINT (182.5 47.2)
## 7 66.4 POINT (166.1 48.8)
## 8 17.7 POINT (160.7 42.1)
## 9 21.9 POINT (162.9 29)
## 10 25.7 POINT (166.4 33.6)
Peta Batas Administratif
map_local <- st_as_sf(st_as_sfc(st_bbox(c(xmin=0, ymin=0, xmax=200, ymax=200))))
st_crs(map_local) <- 32617
tmap_mode("plot")
## ℹ tmap modes "plot" - "view"
## ℹ toggle with `tmap::ttm()`
tm_shape(map_local) + tm_polygons(alpha = 0.3) +
tm_shape(longleaf_sf) +
tm_dots("vble", palette = "viridis")
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
grid <- rast(map_local, nrows = 100, ncols = 100)
xy <- xyFromCell(grid, 1:ncell(grid))
coop <- st_as_sf(as.data.frame(xy), coords = c("x","y"), crs = 32617)
qtm(coop)
Metode Deterministik
Closest Observation
v <- voronoi(x = vect(longleaf_sf), bnd = map_local)
plot(v)
points(vect(longleaf_sf), cex = 0.5)
v <- st_as_sf(v)
tm_shape(v) + tm_fill(col = "vble", palette = "viridis", alpha = 0.6)
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the map variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
resp <- st_intersection(v, coop)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
resp$pred <- resp$vble
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c1<-tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis") + tm_title("Closest Observation")
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
Inverse Distant Weight
res <- gstat(formula = vble ~ 1, locations = longleaf_sf, nmax = nrow(longleaf_sf), set = list(idp = 1))
resp <- predict(res, coop)
## [inverse distance weighted interpolation]
resp$x <- st_coordinates(resp)[,1]
resp$y <- st_coordinates(resp)[,2]
resp$pred <- resp$var1.pred
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c2<-tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis") + tm_title("IDW")
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
Nearest Neighbors
res <- gstat(formula = vble ~ 1, locations = longleaf_sf, nmax = 6, set = list(idp = 0))
resp <- predict(res, coop)
## [inverse distance weighted interpolation]
resp$x <- st_coordinates(resp)[,1]
resp$y <- st_coordinates(resp)[,2]
resp$pred <- resp$var1.pred
pred <- rasterize(resp, grid, field = "pred", fun = "mean")
c3<-tm_shape(pred) + tm_raster(alpha = 0.6, palette = "viridis") + tm_title("KNN")
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
Pendekatan Ensemble
Closest Observation
v <- voronoi(x = vect(longleaf_sf), bnd = vect(map_local))
v <- st_as_sf(v)
idx <- st_intersects(coop, v)
p1 <- vapply(idx,
function(i) if (length(i) == 0) NA_real_ else v$vble[i[1]],
numeric(1))
IDW
gs<-gstat(formula=vble~1,locations= longleaf_sf,nmax= nrow(longleaf_sf),
set= list(idp=1))
p2<-predict(gs,coop)$var1.pred
## [inverse distance weighted interpolation]
K-Nearest Neighbors
nn<-gstat(formula=vble~1,locations= longleaf_sf,nmax= 6,
set= list(idp=0))
p3<-predict(nn,coop)$var1.pred
## [inverse distance weighted interpolation]
Ensemble
weights<-c(1/3,1/3,1/3)
p4<-p1*weights[1]+p2*weights[2]+p3*weights[3]
resp <- data.frame(
x = st_coordinates(coop)[,1],
y = st_coordinates(coop)[,2],
pred = p4
)
resp <- st_as_sf(resp, coords = c("x","y"), crs = 32617)
pred4 <- rasterize(vect(resp), grid, field = "pred", fun = "mean")
c4 <- tm_shape(pred4) + tm_raster(alpha = 0.6, palette = "viridis") +
tm_title("Ensemble")
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## map variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_raster()`: use `col_alpha` instead of `alpha`.
tmap_arrange(c1, c2, c3, c4, ncol = 2)
Metode Stokastik (Kriging)
mapview(longleaf_sf, zcol = "vble", map.types = "CartoDB.Voyager")
mapview(coop, map.types = "CartoDB.Voyager")
## Warning in cbind(`Feature ID` = fid, mat): number of rows of result is not a
## multiple of vector length (arg 1)
vc <- variogram(log(vble) ~ 1, longleaf_sf, cloud = TRUE)
plot(vc)
v <- variogram(log(vble) ~ 1, data = longleaf_sf)
plot(v)
## 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)
show.vgms(par.strip.text = list(cex = 0.75))
vinitial <- vgm(psill = 0.76, model = "Sph", range = 30, nugget = 0.2)
plot(v, vinitial, cutoff = 100, cex = 1.5)
fv<-fit.variogram(object=v,model= vgm(psill= 0.76,model= "Sph",range= 30,nugget= 0.2))
fv
## model psill range
## 1 Nug 0.03703298 0.00000
## 2 Sph 0.90397406 28.24346
afvg <- autofitVariogram(log(vble) ~ 1, longleaf_sf)
plot(afvg)
k <- gstat(formula = log(vble) ~ 1, data = longleaf_sf, model = fv)
kpred <- predict(k, coop)
## [using ordinary kriging]
ggplot() + geom_sf(data = kpred, aes(color = var1.pred)) +
geom_sf(data = longleaf_sf) +
scale_color_viridis(name = "log(vble)") + theme_bw()
ggplot() + geom_sf(data = kpred, aes(color = var1.var)) +
geom_sf(data = longleaf_sf) +
scale_color_viridis(name = "variance") + theme_bw()
Evaluasi Model
Cross Validation
set.seed(048)
RMSE <- function(observed, predicted) {
sqrt(mean((observed - predicted)^2))
}
kf <- kfold(nrow(longleaf_sf), k = 5)
rmse1 <- rep(NA, 5)
rmse2 <- rep(NA, 5)
rmse3 <- rep(NA, 5)
rmse4 <- rep(NA, 5)
rmse5 <- rep(NA, 5)
for(k in 1:5){
test<-longleaf_sf[kf== k,]
train<-longleaf_sf[kf !=k,]
v <- voronoi(x = vect(train), bnd = vect(map_local))
v <- st_as_sf(v)
idx <- st_intersects(test, v)
p1 <- vapply(idx,
function(i) if (length(i) == 0) NA_real_ else v$vble[i[1]],
numeric(1))
rmse1[k] <- RMSE(test$vble, p1)
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)
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)
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)
afvg_k <- tryCatch(
autofitVariogram(log(vble) ~ 1, train, cutoff = 100),
error = function(e) NULL
)
if (!is.null(afvg_k)) {
k_k <- gstat(formula = log(vble) ~ 1, data = train,
model = afvg_k$var_model)
p5_log <- predict(k_k, test)$var1.pred
p5 <- exp(p5_log)
rmse5[k] <- RMSE(test$vble, p5)
} else {
p5 <- rep(NA_real_, nrow(test))
rmse5[k] <- NA
}
}
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [using ordinary kriging]
Tabel RMSE
data.frame(closest.obs= rmse1,IDW= rmse2,
nearest.neighbors= rmse3,ensemble= rmse4, kriging= rmse5)
## closest.obs IDW nearest.neighbors ensemble kriging
## 1 14.05815 16.94283 12.61120 12.78885 12.64418
## 2 12.20609 14.88017 11.22873 10.82253 12.15022
## 3 14.42222 14.17224 10.19146 11.00566 12.14936
## 4 16.89717 17.28433 15.10307 14.78798 14.87530
## 5 12.14880 15.32086 11.04565 10.94442 11.60855
Tabel Rata-Rata RMSE
data.frame(
closest.obs = mean(rmse1, na.rm = TRUE),
IDW = mean(rmse2, na.rm = TRUE),
nearest.neigh = mean(rmse3, na.rm = TRUE),
ensemble = mean(rmse4, na.rm = TRUE),
kriging = mean(rmse5, na.rm = TRUE)
)
## closest.obs IDW nearest.neigh ensemble kriging
## 1 13.94648 15.72009 12.03602 12.06989 12.68552