inst/doc/GRTS.R

### R code from vignette source 'GRTS.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: foo
###################################################
n <- 19
reps <- 1000

library(GRTS)
library(xtable)
library(plyr)
library(ggplot2)
library(reshape)
library(sp)
theme_set(theme_grey(8))
options(width = 60)
options(str = strOptions(strict.width = "cut"))
foo <- packageDescription("GRTS")


###################################################
### code chunk number 2: GRTS.Rnw:56-58
###################################################
tmp <- matrix(paste("(", rep(7:0, 8), ",", rep(0:7, each = 8), ")", sep = ""), ncol = 8, nrow = 8)
print(xtable(tmp, align = "r|rrrrrrrr|", label = "tab:2D", caption = "2D coordinates of a 8 x 8 matrix"), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 8))


###################################################
### code chunk number 3: GRTS.Rnw:65-67
###################################################
tmp <- matrix(paste("(", rep(rep(1:0, each = 4), 8), ",", rep(rep(0:1, each = 4), each = 8), ")", sep = ""), ncol = 8, nrow = 8)
print(xtable(tmp, align = "r|rrrr|rrrr|", label = "tab:Level1_2D", caption = "Binary 2D codes for the first level of submatrices."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 4, 8))


###################################################
### code chunk number 4: GRTS.Rnw:72-74
###################################################
tmp <- matrix(c(rep(rep(c(1, 0), each = 4), 4), rep(rep(c(3, 2), each = 4), 4)), ncol = 8, nrow = 8)
print(xtable(tmp, align = "r|rrrr|rrrr|", digits = 0, label = "tab:Level1_1D", caption = "Base 4 1D codes for the first level of submatrices."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 4, 8))


###################################################
### code chunk number 5: GRTS.Rnw:79-82
###################################################
tmp2 <- matrix(c(rep(rep(c(1, 0), each = 2), 4), rep(rep(c(3, 2), each = 2), 4)), ncol = 8, nrow = 8)
print(xtable(tmp2, align = "r||rr|rr||rr|rr||", digits = 0, label = "tab:Level2_1D", caption = "Base 4 1D codes for the second level of submatrices."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))
print(xtable(matrix(sprintf("%02i", 10 * tmp2 + tmp), ncol = 8), align = "r||rr|rr||rr|rr||", label = "tab:Level2_1Dc", caption = "Prepending the base 4 1D codes for the second level of submatrices to the base 4 1D codes for the first level of submatrices."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))


###################################################
### code chunk number 6: GRTS.Rnw:87-90
###################################################
tmp3 <- matrix(c(rep(1:0, 4), rep(3:2, 4)), ncol = 8, nrow = 8)
print(xtable(tmp3, align = "r||rr|rr||rr|rr||", digits = 0, label = "tab:Level3_1D", caption = "Base 4 1D codes for the third level of submatrices."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))
print(xtable(matrix(sprintf("%03i", 100 * tmp3 + 10 * tmp2 + tmp), ncol = 8), label = "tab:Level3_1Dc", caption = "Combining the base 4 1D codes from all levels.", align = "r||rr|rr||rr|rr||"), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))


###################################################
### code chunk number 7: GRTS.Rnw:98-105
###################################################
tmp <- GRTS(4)
tmp2 <- tmp %/% 4
tmp1 <- tmp %% 4
print(xtable(tmp1, digits = 0, align = "r|rr|rr|", label = "tab:Ex4L1", caption = "Level 1 submatrices with randomised base 4 1D code."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 2, 4))
print(xtable(tmp2, digits = 0, align = "r|rr|rr|", label = "tab:Ex4L2", caption = "Level 2 submatrices with randomised base 4 1D code."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 2, 4))
print(xtable(matrix(sprintf("%02i", 10 * tmp2 + tmp1), ncol = 4), align = "r|rr|rr|", label = "tab:Ex4LC", caption = "Combined base 4 1D codes."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 2, 4))
print(xtable(tmp, digits = 0, align = "r|rr|rr|", label = "tab:Ex4L", caption = "Order of randomised base 4 1D code."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 2, 4))


###################################################
### code chunk number 8: GRTS.Rnw:110-119
###################################################
tmp <- GRTS(8)
tmp3 <- tmp %/% 16
tmp2 <- tmp %/% 4 %% 4
tmp1 <- tmp %% 4
print(xtable(tmp1, digits = 0, align = "r||rr|rr||rr|rr||", label = "tab:Ex8L1", caption = "Level 1 submatrices with randomised base 4 1D code."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))
print(xtable(tmp2, digits = 0, align = "r||rr|rr||rr|rr||", label = "tab:Ex8L2", caption = "Level 2 submatrices with randomised base 4 1D code."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))
print(xtable(tmp3, digits = 0, align = "r||rr|rr||rr|rr||", label = "tab:Ex8L3", caption = "Level 3 submatrices with randomised base 4 1D code."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))
print(xtable(matrix(sprintf("%03i", 100 * tmp3 + 10 * tmp2 + tmp1), ncol = 8), align = "r||rr|rr||rr|rr||", label = "tab:Ex8LC", caption = "Combined base 4 1D codes."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))
print(xtable(tmp, digits = 0, align = "r||rr|rr||rr|rr||", label = "tab:Ex8L", caption = "Order of randomised base 4 1D code."), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))


###################################################
### code chunk number 9: GRTS.Rnw:128-129
###################################################
print(xtable(ifelse(tmp < n, "X", ""), align = "c||cc|cc||cc|cc||", label = "tab:Ex8Sample", caption = paste("A sample of", n , "points for table \\ref{tab:Ex8L}")), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))


###################################################
### code chunk number 10: GRTS.Rnw:136-141
###################################################
samples <- replicate(reps, GRTS(8) < n)
p3 <- apply(samples, 1:2, sum)
print(xtable(matrix(sprintf("%2.1f%%", 100 * p3 / reps), ncol = 8), align = "c||cc|cc||cc|cc||", label = "tab:PropL3", caption = paste("Proportion of", reps , "replications in which the grid cell is selected when sampling", n, "grid cells using GRTS. The expected proportion is", 
sprintf("$\\frac{%i}{8^2}=%2.1f\\%%$", n, 100 * n/(8^2))
)), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(0, 0, 2, 4, 4, 6, 8, 8))


###################################################
### code chunk number 11: GRTS.Rnw:148-156
###################################################
dataset <- expand.grid(X = 0:7, Y = 0:7)
dataset$X1 <- factor(dataset$X %/% 4, levels = 1:0)
dataset$Y1 <- factor(dataset$Y %/% 4)
dataset$L2 <- with(dataset, factor(X %% 4 %/% 2 + 2 * Y %% 4 %/% 2))
tmp <- ddply(dataset, .(X1, Y1), function(x){
  data.frame(Samples = apply(samples[unique(x$X + 1), unique(x$Y + 1), ], 3, sum))
})
print(ggplot(tmp, aes(x = Samples)) + geom_histogram(binwidth = 1) + facet_grid(X1 ~ Y1))


###################################################
### code chunk number 12: GRTS.Rnw:164-170
###################################################
dataset$X2 <- dataset$X1:factor(dataset$X %% 4 %/% 2, levels = 1:0)
dataset$Y2 <- dataset$Y1:factor(dataset$Y %% 4 %/% 2)
tmp <- ddply(dataset, .(X2, Y2), function(x){
  data.frame(Samples = apply(samples[unique(x$X + 1), unique(x$Y + 1), ], 3, sum))
})
print(ggplot(tmp, aes(x = Samples)) + geom_histogram(binwidth = 1) + facet_grid(X2 ~ Y2))


###################################################
### code chunk number 13: GRTS.Rnw:182-184
###################################################
GRTS(8)
GRTS(7)


###################################################
### code chunk number 14: GRTS.Rnw:191-202
###################################################
#define a SpatialPolygons object
Sr1 <- Polygon(cbind(c(2, 4, 4, 1,2), c(2, 3, 5, 4, 2)))
Sr2 <- Polygon(cbind(c(5, 4, 2, 5), c(1.5, 2.5, 1.5, 1.5)))
Sr3 <- Polygon(cbind(c(4, 4, 5, 10.1, 4), c(5, 3, 2, 5.1, 5)))
Sr4 <- Polygon(cbind(c(4.5, 5.5, 6, 5.5, 4.5), c(4, 3, 3, 4, 4)), hole = TRUE)

Srs1 <- Polygons(list(Sr1), "s1")
Srs2 <- Polygons(list(Sr2), "s2")
Srs3 <- Polygons(list(Sr3, Sr4), "s3/4")
SpP <- SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3)
plot(SpP, col = 1:3, pbg="white", axes = TRUE)


###################################################
### code chunk number 15: GRTS.Rnw:207-214
###################################################
pls <- list("sp.polygons", SpP, col = "black", first = FALSE)
output <- GRTS(SpP, cellsize = 0.1)
limits <- apply(cbind(bbox(output), bbox(SpP)), 1, function(x){
  range(pretty(x))
})
spplot(output, sp.layout = list(pls), col.regions = terrain.colors(100), 
       scales = list(draw = TRUE), xlim = limits[, 1], ylim = limits[, 2])


###################################################
### code chunk number 16: GRTS.Rnw:219-225
###################################################
output <- GRTS(SpP, cellsize = 0.5)
limits <- apply(cbind(bbox(output), bbox(SpP)), 1, function(x){
  range(pretty(x))
})
spplot(output, sp.layout = list(pls), col.regions = terrain.colors(100), 
       scales = list(draw = TRUE), xlim = limits[, 1], ylim = limits[, 2])


###################################################
### code chunk number 17: GRTS.Rnw:230-236
###################################################
output <- GRTS(SpP, cellsize = 0.5, RandomStart = TRUE)
limits <- apply(cbind(bbox(output), bbox(SpP)), 1, function(x){
  range(pretty(x))
})
spplot(output, sp.layout = list(pls), col.regions = terrain.colors(100), 
       scales = list(draw = TRUE), xlim = limits[, 1], ylim = limits[, 2])


###################################################
### code chunk number 18: GRTS.Rnw:241-248
###################################################
output <- GRTS(SpP, cellsize = 0.1, Subset = TRUE)
limits <- apply(cbind(bbox(output), bbox(SpP)), 1, function(x){
  range(pretty(x))
})
spplot(output, sp.layout = list(pls), scales = list(draw = TRUE), 
       col.regions = terrain.colors(100),
       xlim = limits[, 1], ylim = limits[, 2])


###################################################
### code chunk number 19: GRTS.Rnw:253-261
###################################################
n <- 19
#calculate the treshold value
MaxRanking <- max(head(sort(output$Ranking), n))
#do the selection
Selection <- subset(output, Ranking <= MaxRanking)
spplot(Selection, sp.layout = list(pls), scales = list(draw = TRUE), 
       col.regions = rainbow(n),
       xlim = limits[, 1], ylim = limits[, 2])


###################################################
### code chunk number 20: GRTS.Rnw:268-305 (eval = FALSE)
###################################################
## testGRTS <- t(replicate(reps, {
##   #do the randomisation
##   output <- GRTS(SpP, cellsize = 0.1, Subset = TRUE, RandomStart = TRUE)
##   #calculate the treshold value
##   MaxRanking <- max(head(sort(output$Ranking), n))
##   #do the selection
##   Selection <- subset(output, Ranking <= MaxRanking)
##   #do the overlay
##   table(Polygon = factor(over(Selection, SpP), levels = 1:3, 
##     labels = c("A", "B", "C")))
## }))
## testGRTS <- melt(data = testGRTS)
## testGRTS$Type <- "GRTS"
## 
## testSRS <- t(replicate(reps, {
##   Selection <- spsample(SpP, n = n, type = "random")
##   table(Polygon = factor(over(Selection, SpP), levels = 1:3, 
##     labels = c("A", "B", "C")))
## }))
## testSRS <- melt(data = testSRS)
## testSRS$Type <- "SRS"
## 
## test <- rbind(testGRTS, testSRS)
## 
## areas <- sapply(SpP@polygons, function(x){
##   tmp <- sapply(x@Polygons, function(y){
##     c(ifelse(y@hole, -1, 1), y@area)
##   })
##   sum(tmp[1, ] * tmp[2, ])
## })
## reference <- data.frame(Polygon = factor(c("A", "B", "C")), 
##   Expected = n * areas / sum(areas))
## ggplot(test) + 
##   geom_density(aes(x = value, colour = Type), adjust = 2) + 
##   geom_vline(data = reference, aes(xintercept = Expected), colour = "blue") +
##   facet_grid(Polygon ~ .) + 
##   xlab("Number of samples per polygon")


###################################################
### code chunk number 21: GRTS.Rnw:308-315 (eval = FALSE)
###################################################
## p <- 
##   ggplot(test) + 
##   geom_density(aes(x = value, colour = Type), adjust = 3) + 
##   geom_vline(data = reference, aes(xintercept = Expected), colour = "blue") +
##   facet_grid(Polygon ~ .) + 
##   xlab("Number of samples per polygon")
## ggsave(p, filename = "GRTS_SRS_distribution.pdf", width = 4, height = 4, path = 'pkg/GRTS/vignettes')


###################################################
### code chunk number 22: GRTS.Rnw:322-362 (eval = FALSE)
###################################################
## Weights <- data.frame(
##   ID = c("s1", "s2", "s3/4"), 
##   Weight = c(1, 5, 2)
## )
## rownames(Weights) <- Weights$ID
## Weights$Expected <- n * areas * Weights$Weight / sum(areas * Weights$Weight)
## SpP <- SpatialPolygonsDataFrame(SpP, data = Weights)
## SpP$ID <- factor(SpP$ID)
## test <- replicate(reps, {
##   GRTSorder <- GRTS(SpP, cellsize = 0.1, Subset = TRUE)
##   GRTSorder$Weight <- over(GRTSorder, SpP[, "Weight"])$Weight / (max(SpP$Weight))
##   GRTSorder$Weight <- 
##     GRTSorder$Weight * 
##     pmin(
##       1, 
##       length(GRTSorder) / sum(GRTSorder$Weight)
##     )
##   GRTSups <- 
##     GRTSorder[
##       rbinom(
##         length(GRTSorder), 
##         size = 1, 
##         prob = GRTSorder$Weight
##       ) == 1, 
##       "Ranking"]
##   table(
##     over(
##       GRTSups[order(GRTSups$Ranking) <= n, ], 
##       SpP[, "ID"]
##     )$ID,
##     useNA = 'ifany')
## })
## test <- melt(t(test))
## colnames(test) <- c("Run", "ID", "Estimate")
## test <- merge(test, Weights)
## ggplot(test, aes(x = Estimate)) + 
##   geom_histogram(binwidth = 1) + 
##   geom_vline(aes(xintercept = Expected), colour = "red") + 
##   xlab("Number of samples per polygon") + 
##   facet_wrap(~ID)


###################################################
### code chunk number 23: GRTS.Rnw:365-368 (eval = FALSE)
###################################################
## p <- 
##   ggplot(test, aes(x = Estimate)) + geom_histogram(binwidth = 1) + geom_vline(aes(xintercept = Expected), colour = "red") + xlab("Number of samples per polygon") + facet_wrap(~ID)
## ggsave(p, filename = "GRTS_UPS_distribution.pdf", width = 4, height = 4, path = 'pkg/GRTS/vignettes')

Try the GRTS package in your browser

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

GRTS documentation built on May 2, 2019, 6:39 p.m.