inst/doc/spcosa.R

In spcosa: Spatial Coverage Sampling and Random Sampling from Compact Geographical Strata

## ---- echo=FALSE--------------------------------------------------------------
knitr::opts_chunk$set(comment = NA)

## -----------------------------------------------------------------------------
library(spcosa)

## -----------------------------------------------------------------------------
library(ggplot2)

## -----------------------------------------------------------------------------
set.seed(314)

## ---- message=FALSE-----------------------------------------------------------
library(sp)
grd <- expand.grid(s1 = 1:100, s2 = 1:50)
gridded(grd) <- ~ s1 * s2

## -----------------------------------------------------------------------------
stratification <- stratify(grd, nStrata = 75, nTry = 10)

## ---- fig.width=7, fig.height=4, out.width=500--------------------------------
plot(stratification)

## ---- fig.width=7, fig.height=4, out.width=500, message=FALSE-----------------
plot(stratification) +
    scale_x_continuous(name = "Easting (km)") +
    scale_y_continuous(name = "Northing (km)")

## -----------------------------------------------------------------------------
sampling_pattern <- spsample(stratification)

## ---- fig.width=7, fig.height=4, out.width=500--------------------------------
plot(sampling_pattern)

## ---- fig.width=7, fig.height=4, out.width=500--------------------------------
plot(stratification, sampling_pattern)

## -----------------------------------------------------------------------------
sampling_points <- as(sampling_pattern, "data.frame")
head(sampling_points)

## ---- echo=FALSE--------------------------------------------------------------
prior_points <- spsample(grd, n = 50, type = "random")
prior_points <- as(prior_points, "data.frame")
names(prior_points) <- c("s1", "s2")

## -----------------------------------------------------------------------------
head(prior_points)

## -----------------------------------------------------------------------------
coordinates(prior_points) <- ~ s1 * s2
stratification <- stratify(grd, priorPoints = prior_points, nStrata = 75, nTry = 100)
sampling_pattern <- spsample(stratification)

## ---- fig.width=7, fig.height=4, out.width=500--------------------------------
plot(stratification, sampling_pattern)

## -----------------------------------------------------------------------------
stratification <- stratify(grd, nStrata = 25, nTry = 10)

## ---- fig.width=7, fig.height=4, out.width=500--------------------------------
plot(stratification)

## -----------------------------------------------------------------------------
sampling_pattern <- spsample(stratification, n = 2)

## ---- fig.width=7, fig.height=4, out.width=500--------------------------------
plot(stratification, sampling_pattern)

## -----------------------------------------------------------------------------
sampling_points <- as(sampling_pattern, "data.frame")
head(sampling_points)

## -----------------------------------------------------------------------------
my_data <- data.frame(clay = rbeta(50, 2, 25), SOM = rbeta(50, 1, 50))
head(my_data)

## -----------------------------------------------------------------------------
estimate("spatial mean", stratification, sampling_pattern, my_data)

## -----------------------------------------------------------------------------
estimate("standard error", stratification, sampling_pattern, my_data)

## -----------------------------------------------------------------------------
scdf <- estimate("scdf", stratification, sampling_pattern, my_data)

## -----------------------------------------------------------------------------
head(scdf$clay)
head(scdf$SOM)

## ---- fig.width=6, fig.height=5, out.width=400, echo=FALSE--------------------
tmp <- as.data.frame(scdf$clay)
tmp$value <- tmp$value * 0.01
dx <- mean(diff(tmp$value))
dat <- data.frame(
    x = c(tmp$value[1] - dx, tmp$value),
    xend = c(tmp$value, tmp$value[nrow(tmp)]+ dx),
    yend =  c(tmp$cumFreq, 1)
)
ggplot() +
    geom_segment(data = dat, mapping = aes(x = x, y = yend, xend = xend, yend = yend)) +
    geom_point(data = dat[-1, ], mapping = aes(x = x, y = yend), shape = 16) +
        scale_x_continuous(name = "x : clay content (g/dag)") +
        scale_y_continuous(name = "cumulative frequency")
tmp <-  unique(as.data.frame(scdf$SOM))
tmp$value <- tmp$value * 0.01
dx <- mean(diff(tmp$value))
dat <- data.frame(
    x = c(tmp$value[1] - dx, tmp$value),
    xend = c(tmp$value, tmp$value[nrow(tmp)]+ dx),
    yend =  c(tmp$cumFreq, 1)
)
ggplot() +
    geom_segment(data = dat, mapping = aes(x = x, y = yend, xend = xend, yend = yend)) +
    geom_point(data = dat[-1, ], mapping = aes(x = x, y = yend), shape = 16) +
    scale_x_continuous(name = "x : soil organic matter content (g/dag)") +
    scale_y_continuous(name = "cumulative frequency")

## -----------------------------------------------------------------------------
library(sf)
directory <- system.file("maps", package = "spcosa")
shp_farmsum <- as(st_read(dsn = directory, layer = "farmsum"), "Spatial")

## ---- fig.width=6, fig.height=6, out.width=400--------------------------------
stratification <- stratify(shp_farmsum, nStrata = 20, equalArea = TRUE, nTry = 10)
plot(stratification)

## ---- fig.width=6, fig.height=6, out.width=400--------------------------------
sampling_pattern <- spsample(stratification, n = 2, type = "composite")
plot(stratification, sampling_pattern)

## -----------------------------------------------------------------------------
sampling_points <- as(sampling_pattern, "data.frame")
head(sampling_points)

## ---- echo=FALSE--------------------------------------------------------------
my_data <- data.frame(
  clay = c(9.7, 10.4),
  SOM = c(4.9, 5.2)
)

## -----------------------------------------------------------------------------
estimate("spatial mean", stratification, sampling_pattern, my_data)
estimate("standard error", stratification, sampling_pattern, my_data)

## -----------------------------------------------------------------------------
sampling_pattern <- spsample(stratification, n = 2)
sampling_points <- as(sampling_pattern, "data.frame")
head(sampling_points)

## ---- fig.width=6, fig.height=6, out.width=400--------------------------------
plot(stratification, sampling_pattern)

## -----------------------------------------------------------------------------
grd <- expand.grid(
    longitude = seq(from = -176, to = 180, by = 4),
    latitude  = seq(from =  -86, to =  86, by = 4)
)
gridded(grd) <- ~ longitude * latitude

grd_crs <- grd
slot(grd_crs, "proj4string") <- CRS("EPSG:4326")


## -----------------------------------------------------------------------------
strata     <- stratify(grd,     nStrata = 50)
strata_crs <- stratify(grd_crs, nStrata = 50)

## ---- fig.width=8, fig.height=5, out.width=500--------------------------------
plot(strata)
plot(strata_crs)

## -----------------------------------------------------------------------------
shp_mijdrecht <- as(st_read(
    dsn = system.file("maps", package = "spcosa"), 
    layer = "mijdrecht"), "Spatial")
stratification <- stratify(shp_mijdrecht, nStrata = 1, nGridCells = 5000)
sampling_pattern <- spsample(stratification, n = 30)

## ---- fig.width=5, fig.height=7, out.width=300--------------------------------
plot(stratification, sampling_pattern)

## -----------------------------------------------------------------------------
doughnut <- expand.grid(s1 = -25:25, s2 = -25:25)
d <- with(doughnut, sqrt(s1^2 + s2^2))
doughnut <- doughnut[(d < 25) & (d > 15), ]
coordinates(doughnut) <- ~ s1 * s2
gridded(doughnut) <- TRUE

## ---- fig.width=5, fig.height=5, out.width=350--------------------------------
stratification <- stratify(doughnut, nStrata = 2, nTry = 100)
sampling_pattern <- spsample(stratification)
plot(stratification, sampling_pattern)

## ---- echo=FALSE--------------------------------------------------------------
sessionInfo()