inst/doc/spatial.R

## ----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"))

Try the SamplingStrata package in your browser

Any scripts or data that you put into this service are public.

SamplingStrata documentation built on Nov. 16, 2022, 1:08 a.m.