Nothing
## ----setup, include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
options(width = 999)
knitr::opts_chunk$set(echo = TRUE, fig.width=6, fig.height=4, out.width = 6, out.heigth = 4)
## ----load, eval=TRUE, echo=FALSE, include=FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
load("spatial_vignette.RData")
## ---- eval = FALSE,echo=FALSE, message=FALSE-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# library(sp)
# # locations (155 observed points)
# data("meuse")
# # grid of points (3103)
# data("meuse.grid")
# meuse.grid$id <- c(1:nrow(meuse.grid))
# coordinates(meuse)<-c("x","y")
# coordinates(meuse.grid)<-c("x","y")
## ---- eval = FALSE,echo=FALSE,include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# par(mfrow=c(1,2))
# plot(meuse.grid)
# title("Meuse territory (3103 points)",sub="with reported distance from the river",cex.main=0.8,cex.sub=0.8)
# plot(meuse)
# title("Subset of 155 points",sub="with also observed metals concentration",cex.main=0.8,cex.sub=0.8)
# par(mfrow=c(1,1))
## ---- eval = F,echo=TRUE,message=FALSE,fig.width=6,fig.height=8----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# library(automap)
# kriging_lead = autoKrige(log(lead) ~ dist, meuse, meuse.grid)
# plot(kriging_lead,sp.layout = NULL, justPosition = TRUE)
## ---- eval = F,echo=TRUE,message=FALSE,fig.width=6,fig.height=8----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# kriging_zinc = autoKrige(log(zinc) ~ dist, meuse, meuse.grid)
# plot(kriging_zinc, sp.layout = list(pts = list("sp.points", meuse)))
## ---- eval = F,echo=TRUE,message=FALSE,fig.width=6,fig.height=8----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# r2_lead <- 1 - kriging_lead$sserr/sum((log(meuse$lead)-mean(log(meuse$lead)))^2)
# r2_lead
# ## [1] 0.9999997
# r2_zinc <- 1 - kriging_zinc$sserr/sum((log(meuse$zinc)-mean(log(meuse$zinc)))^2)
# r2_zinc
# ## [1] 0.9999999
## ---- eval = F,echo=TRUE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# df <- NULL
# df$id <- meuse.grid$id
# df$lead.pred <- kriging_lead$krige_output@data$var1.pred
# df$lead.var <- kriging_lead$krige_output@data$var1.var
# df$zinc.pred <- kriging_zinc$krige_output@data$var1.pred
# df$zinc.var <- kriging_zinc$krige_output@data$var1.var
# df$lon <- meuse.grid$x
# df$lat <- meuse.grid$y
# df$dom1 <- 1
# df <- as.data.frame(df)
# head(df)
# ## id lead.pred lead.var zinc.pred zinc.var lon lat dom1
# ## 1 1 5.509360 0.1954937 6.736502 0.2007150 181180 333740 1
# ## 2 2 5.546006 0.1716895 6.785460 0.1749260 181140 333700 1
# ## 3 3 5.488913 0.1784052 6.698883 0.1826314 181180 333700 1
# ## 4 4 5.388320 0.1855561 6.558216 0.1906426 181220 333700 1
# ## 5 5 5.584415 0.1463018 6.841612 0.1465346 181100 333660 1
# ## 6 6 5.525538 0.1533757 6.749216 0.1549663 181140 333660 1
## ---- eval = F,echo=TRUE, message=FALSE, warning=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# library(SamplingStrata)
# frame <- buildFrameSpatial(df=df,
# id="id",
# X=c("lead.pred","zinc.pred"),
# Y=c("lead.pred","zinc.pred"),
# variance=c("lead.var","zinc.var"),
# lon="lon",
# lat="lat",
# domainvalue = "dom1")
# cv <- as.data.frame(list(DOM=rep("DOM1",1),
# CV1=rep(0.01,1),
# CV2=rep(0.01,1),
# domainvalue=c(1:1) ))
# cv
# # DOM CV1 CV2 domainvalue
# # 1 DOM1 0.01 0.01 1
## ---- eval = F,echo=TRUE, message=FALSE, warning=FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# set.seed(1234)
# solution <- optimStrata (
# method = "spatial",
# errors=cv,
# framesamp=frame,
# iter = 15,
# pops = 10,
# nStrata = 5,
# fitting = c(r2_lead,r2_zinc),
# range = c(kriging_lead$var_model$range[2],kriging_zinc$var_model$range[2]),
# kappa=1,
# writeFiles = FALSE,
# showPlot = FALSE,
# parallel = FALSE)
# framenew <- solution$framenew
# outstrata <- solution$aggr_strata
# ##
# ## Input data have been checked and are compliant with requirements
# ## Sequential optimization as parallel = FALSE, defaulting number of cores = 1
# ## *** Domain : 1 1
# ## Number of strata : 3103
# ## *** Sample cost: 61.92901
# ## *** Number of strata: 4
## ---- eval = F,echo=TRUE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# plotStrata2d(framenew,outstrata,domain=1,vars=c("X1","X2"),
# labels=c("Lead","Zinc"))
## ---- eval = F,echo=TRUE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# frameres <- SpatialPointsDataFrame(data=framenew, coords=cbind(framenew$LON,framenew$LAT) )
# frameres2 <- SpatialPixelsDataFrame(points=frameres[c("LON","LAT")], data=framenew)
# frameres2$LABEL <- as.factor(frameres2$LABEL)
# spplot(frameres2,c("LABEL"), col.regions=bpy.colors(5))
## ---- eval = F,echo=TRUE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# s <- selectSampleSpatial(framenew,outstrata,coord_names=c("LON","LAT"))
# ##
# ## *** Sample has been drawn successfully ***
# ## 62 units have been selected from 4 strata
## ---- eval = F,echo=TRUE,message=FALSE,fig.width=6, fig.height=8---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# s <- selectSampleSpatial(framenew,outstrata,coord_names = c("LON","LAT"))
# coordinates(s) <- ~LON+LAT
# proj4string(s) <- CRS("+init=epsg:28992")
# s$LABEL <- as.factor(s$LABEL)
# library(mapview)
# mapview(s,zcol="LABEL", map.types = c("OpenStreetMap"))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.