if (FALSE) {
#############################
# Example of "spatial" method
#############################
library(sp)
data("meuse")
data("meuse.grid")
meuse.grid$id <- c(1:nrow(meuse.grid))
coordinates(meuse) <- c('x','y')
coordinates(meuse.grid) <- c('x','y')
library(gstat)
library(automap)
v <- variogram(lead ~ dist + soil, data = meuse)
fit.vgm.lead <- autofitVariogram(lead ~ dist + soil,
meuse,
model = "Exp")
plot(v, fit.vgm.lead$var_model)
lead.kr <- krige(lead ~ dist + soil, meuse, meuse.grid,
model = fit.vgm.lead$var_model)
lead.pred <- ifelse(lead.kr[1]$var1.pred < 0,0, lead.kr[1]$var1.pred)
lead.var <- ifelse(lead.kr[2]$var1.var < 0, 0, lead.kr[2]$var1.var)
df <- as.data.frame(list(dom = rep(1,nrow(meuse.grid)),
lead.pred = lead.pred,
lead.var = lead.var,
lon = meuse.grid$x,
lat = meuse.grid$y,
id = c(1:nrow(meuse.grid))))
frame <-buildFrameSpatial(df = df,
id = "id",
X = c("lead.pred"),
Y = c("lead.pred"),
variance = c("lead.var"),
lon = "lon",
lat = "lat",
domainvalue = "dom")
cv <- as.data.frame(list(DOM = rep("DOM1",1),
CV1 = rep(0.05,1),
domainvalue = c(1:1) ))
solution <- optimizeStrataSpatial(errors = cv,
framesamp = frame,
iter = 25,
pops = 10,
nStrata = 5,
fitting = 1,
kappa = 1,
range = fit.vgm.lead$var_model$range[2])
framenew <- solution$framenew
outstrata <- solution$aggr_strata
frameres <- SpatialPixelsDataFrame(points = framenew[c("LON","LAT")],
data = framenew)
frameres$LABEL <- as.factor(frameres$LABEL)
spplot(frameres,c("LABEL"), col.regions=bpy.colors(5))
s <- selectSample(framenew,outstrata)
}
Run the code above in your browser using DataLab