1### R code from vignette source 'csdacm.Rnw' 2 3################################################### 4### code chunk number 1: csdacm.Rnw:53-56 5################################################### 6owidth <- getOption("width") 7options("width"=90) 8.PngNo <- 0 9 10 11################################################### 12### code chunk number 2: figreset (eval = FALSE) 13################################################### 14## .iwidth <- 5 15## .iheight <- 6 16## .ipointsize <- 12 17 18 19################################################### 20### code chunk number 3: csdacm.Rnw:65-66 21################################################### 22.iwidth <- 5 23.iheight <- 6 24.ipointsize <- 12 25 26 27################################################### 28### code chunk number 4: afig (eval = FALSE) 29################################################### 30## .PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="") 31## pdf(file=file, width = .iwidth, height = .iheight, pointsize = .ipointsize) 32## opar <- par(mar=c(3,3,1,1)+0.1) 33 34 35################################################### 36### code chunk number 5: zfig (eval = FALSE) 37################################################### 38## dev.null <- dev.off() 39## cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="") 40 41 42################################################### 43### code chunk number 6: csdacm.Rnw:109-112 44################################################### 45myfun <- function(x) { 46 x + 2 47} 48 49 50################################################### 51### code chunk number 7: csdacm.Rnw:117-118 52################################################### 53myfun(1:3) 54 55 56################################################### 57### code chunk number 8: csdacm.Rnw:123-124 58################################################### 59myfun(x=1:3) 60 61 62################################################### 63### code chunk number 9: csdacm.Rnw:130-133 64################################################### 65plotXplus2Yminus3 <- function(x, y, ...) { 66 plot(x = x + 2, y = y - 3, ...) 67} 68 69 70################################################### 71### code chunk number 10: csdacm.Rnw:146-147 72################################################### 73methods("plot") 74 75 76################################################### 77### code chunk number 11: csdacm.Rnw:152-154 78################################################### 79library(sp) 80showMethods("plot") 81 82 83################################################### 84### code chunk number 12: csdacm.Rnw:165-168 85################################################### 86x <- rnorm(10) 87class(x) <- "foo" 88x 89 90 91################################################### 92### code chunk number 13: csdacm.Rnw:177-180 93################################################### 94plot.foo <- function(x, y, ...) { 95 plot.default(x, type = 'l', ...) 96} 97 98 99################################################### 100### code chunk number 14: csdacm.Rnw:189-190 101################################################### 102class(x) <- c("foo", "bar") 103 104 105################################################### 106### code chunk number 15: csdacm.Rnw:192-193 (eval = FALSE) 107################################################### 108## plot(x) 109 110 111################################################### 112### code chunk number 16: csdacm.Rnw:204-207 113################################################### 114data(meuse) 115class(meuse) 116class(lm(log(zinc)~sqrt(dist), meuse)) 117 118 119################################################### 120### code chunk number 17: csdacm.Rnw:226-227 121################################################### 122options("width"=60) 123 124 125################################################### 126### code chunk number 18: csdacm.Rnw:229-249 (eval = FALSE) 127################################################### 128## setClass("CRS", representation(projargs = "character")) 129## setClass("Spatial", 130## representation(bbox = "matrix", proj4string = "CRS"), 131## # NOT TOO WIDE 132## validity <- function(object) { 133## bb <- bbox(object) 134## if (!is.matrix(bb)) 135## return("bbox should be a matrix") 136## n <- dimensions(object) 137## if (n < 2) 138## return("spatial.dimension should be 2 or more") 139## if (any(is.na(bb))) 140## return("bbox should never contain NA values") 141## if (any(!is.finite(bb))) 142## return("bbox should never contain infinite values") 143## if (any(bb[,"max"] < bb[,"min"])) 144## return("invalid bbox: max < min") 145## TRUE 146## } 147## ) 148 149 150################################################### 151### code chunk number 19: csdacm.Rnw:251-252 152################################################### 153options("width"=70) 154 155 156################################################### 157### code chunk number 20: csdacm.Rnw:264-265 158################################################### 159isGeneric("show") 160 161 162################################################### 163### code chunk number 21: csdacm.Rnw:271-273 (eval = FALSE) 164################################################### 165## setGeneric("bbox", function(obj) standardGeneric("bbox")) 166## setMethod("bbox", signature = "Spatial", function(obj) obj@bbox) 167 168 169################################################### 170### code chunk number 22: csdacm.Rnw:304-315 171################################################### 172library(sp) 173setClass("trip", representation("SpatialPointsDataFrame", TOR.columns = "character"), 174 validity <- function(object) { 175 if (length(object@TOR.columns) != 2) 176 stop("Time/id column names must be of length 2") 177 if (!all(object@TOR.columns %in% names(object@data))) 178 stop("Time/id columns must be present in attribute table") 179 TRUE 180 } 181) 182showClass("trip") 183 184 185################################################### 186### code chunk number 23: csdacm.Rnw:328-341 187################################################### 188trip.default <- function(obj, TORnames) { 189 if (!is(obj, "SpatialPointsDataFrame")) 190 stop("trip only supports SpatialPointsDataFrame") 191 if (is.numeric(TORnames)) 192 TORnames <- names(obj)[TORnames] 193 new("trip", obj, TOR.columns = TORnames) 194} 195 196if (!isGeneric("trip")) 197 setGeneric("trip", function(obj, TORnames) 198 standardGeneric("trip")) 199 200setMethod("trip", signature(obj = "SpatialPointsDataFrame", TORnames = "ANY"), trip.default) 201 202 203################################################### 204### code chunk number 24: csdacm.Rnw:348-349 205################################################### 206turtle <- read.csv(system.file("external/seamap105_mod.csv", package="sp")) 207 208 209################################################### 210### code chunk number 25: csdacm.Rnw:351-360 211################################################### 212timestamp <- as.POSIXlt(strptime(as.character(turtle$obs_date), "%m/%d/%Y %H:%M:%S"), "GMT") 213turtle <- data.frame(turtle, timestamp = timestamp) 214turtle$lon <- ifelse(turtle$lon < 0, turtle$lon+360, turtle$lon) 215turtle <- turtle[order(turtle$timestamp),] 216coordinates(turtle) <- c("lon", "lat") 217proj4string(turtle) <- CRS("+proj=longlat +ellps=WGS84") 218turtle$id <- c(rep(1, 200), rep(2, nrow(coordinates(turtle)) - 200)) 219turtle_trip <- trip(turtle, c("timestamp", "id")) 220summary(turtle_trip) 221 222 223################################################### 224### code chunk number 26: csdacm.Rnw:372-382 225################################################### 226summary.trip <- function(object, ...) { 227 cat("Object of class \"trip\"\nTime column: ") 228 print(object@TOR.columns[1]) 229 cat("Identifier column: ") 230 print(object@TOR.columns[2]) 231 print(summary(as(object, "Spatial"))) 232 print(summary(object@data)) 233} 234setMethod("summary", "trip", summary.trip) 235summary(turtle_trip) 236 237 238################################################### 239### code chunk number 27: csdacm.Rnw:400-415 240################################################### 241setGeneric("lines", function(x, ...) standardGeneric("lines")) 242setMethod("lines", signature(x = "trip"), 243 function(x, ..., col = NULL) { 244# NOT TOO WIDE 245 tor <- x@TOR.columns 246 if (is.null(col)) { 247 l <- length(unique(x[[tor[2]]])) 248 col <- hsv(seq(0, 0.5, length = l)) 249 } 250 coords <- coordinates(x) 251 lx <- split(1:nrow(coords), x[[tor[2]]]) 252 for (i in 1:length(lx)) 253 lines(coords[lx[[i]], ], col = col[i], ...) 254 } 255) 256 257 258################################################### 259### code chunk number 28: csdacm.Rnw:436-437 260################################################### 261options("width"=50) 262 263 264################################################### 265### code chunk number 29: csdacm.Rnw:439-448 266################################################### 267setClass("SpatialMultiPoints", representation("SpatialLines"), 268 validity <- function(object) { 269 if (any(unlist(lapply(object@lines, function(x) length(x@Lines))) != 1)) 270# NOT TOO WIDE 271 stop("Only Lines objects with one Line element") 272 TRUE 273 } 274) 275SpatialMultiPoints <- function(object) new("SpatialMultiPoints", object) 276 277 278################################################### 279### code chunk number 30: csdacm.Rnw:450-451 280################################################### 281options("width"=70) 282 283 284################################################### 285### code chunk number 31: csdacm.Rnw:458-468 286################################################### 287n <- 5 288set.seed(1) 289x1 <- cbind(rnorm(n),rnorm(n, 0, 0.25)) 290x2 <- cbind(rnorm(n),rnorm(n, 0, 0.25)) 291x3 <- cbind(rnorm(n),rnorm(n, 0, 0.25)) 292L1 <- Lines(list(Line(x1)), ID="mp1") 293L2 <- Lines(list(Line(x2)), ID="mp2") 294L3 <- Lines(list(Line(x3)), ID="mp3") 295s <- SpatialLines(list(L1,L2,L3)) 296smp <- SpatialMultiPoints(s) 297 298 299################################################### 300### code chunk number 32: csdacm.Rnw:477-491 301################################################### 302plot.SpatialMultiPoints <- function(x, ..., pch = 1:length(x@lines), col = 1, cex = 1) { 303 n <- length(x@lines) 304 if (length(pch) < n) 305 pch <- rep(pch, length.out = n) 306 if (length(col) < n) 307 col <- rep(col, length.out = n) 308 if (length(cex) < n) 309 cex <- rep(cex, length.out = n) 310 plot(as(x, "Spatial"), ...) 311 for (i in 1:n) 312 points(x@lines[[i]]@Lines[[1]]@coords, pch = pch[i], col = col[i], cex = cex[i]) 313} 314setMethod("plot", signature(x = "SpatialMultiPoints", y = "missing"), 315 function(x, y, ...) plot.SpatialMultiPoints(x, ...)) 316 317 318################################################### 319### code chunk number 33: csdacm.Rnw:513-525 320################################################### 321cName <- "SpatialMultiPointsDataFrame" 322setClass(cName, representation("SpatialLinesDataFrame"), 323 validity <- function(object) { 324 lst <- lapply(object@lines, function(x) length(x@Lines)) 325 if (any(unlist(lst) != 1)) 326 stop("Only Lines objects with single Line") 327 TRUE 328 } 329) 330SpatialMultiPointsDataFrame <- function(object) { 331 new("SpatialMultiPointsDataFrame", object) 332} 333 334 335################################################### 336### code chunk number 34: csdacm.Rnw:531-536 337################################################### 338df <- data.frame(x1 = 1:3, x2 = c(1,4,2), row.names = c("mp1", "mp2", "mp3")) 339smp_df <- SpatialMultiPointsDataFrame(SpatialLinesDataFrame(smp, df)) 340setMethod("plot", signature(x = "SpatialMultiPointsDataFrame", y = "missing"), 341 function(x, y, ...) plot.SpatialMultiPoints(x, ...)) 342grys <- c("grey10", "grey40", "grey80") 343 344 345################################################### 346### code chunk number 35: csdacm.Rnw:538-539 (eval = FALSE) 347################################################### 348## plot(smp_df, col = grys[smp_df[["x1"]]], pch = smp_df[["x2"]], cex = 2, axes = TRUE) 349 350 351################################################### 352### code chunk number 36: csdacm.Rnw:546-553 353################################################### 354.iwidth <- 6 355.iheight <- 2.5 356.ipointsize <- 10 357.PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="") 358pdf(file=file, width = .iwidth, height = .iheight, pointsize = .ipointsize) 359opar <- par(mar=c(3,3,1,1)+0.1) 360plot(smp_df, col = grys[smp_df[["x1"]]], pch = smp_df[["x2"]], cex = 2, axes = TRUE) 361dev.null <- dev.off() 362cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="") 363.iwidth <- 5 364.iheight <- 6 365.ipointsize <- 12 366 367 368################################################### 369### code chunk number 37: csdacm.Rnw:569-573 370################################################### 371data(meuse.grid) 372gridded(meuse.grid)=~x+y 373xx <- spsample(meuse.grid, type="hexagonal", cellsize=200) 374class(xx) 375 376 377################################################### 378### code chunk number 38: csdacm.Rnw:581-582 379################################################### 380HexPts <- spsample(meuse.grid, type="hexagonal", cellsize=200) 381 382 383################################################### 384### code chunk number 39: csdacm.Rnw:584-585 (eval = FALSE) 385################################################### 386## spplot(meuse.grid["dist"], sp.layout = list("sp.points", HexPts, col = 1)) 387 388 389################################################### 390### code chunk number 40: csdacm.Rnw:587-590 391################################################### 392HexPols <- HexPoints2SpatialPolygons(HexPts) 393df <- over(HexPols, meuse.grid) 394HexPolsDf <- SpatialPolygonsDataFrame(HexPols, df, match.ID = FALSE) 395 396 397################################################### 398### code chunk number 41: csdacm.Rnw:592-593 (eval = FALSE) 399################################################### 400## spplot(HexPolsDf["dist"]) 401 402 403################################################### 404### code chunk number 42: csdacm.Rnw:599-611 405################################################### 406.iwidth <- 6 407.iheight <- 4 408.PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="") 409pdf(file=file, width = .iwidth, height = .iheight, pointsize = .ipointsize) 410opar <- par(mar=c(3,3,1,1)+0.1) 411library(lattice) 412# RSB quietening greys 413grys <- grey.colors(11, 0.95, 0.55, 2.2) 414print(spplot(meuse.grid["dist"], cuts=10, col.regions=grys, sp.layout = list("sp.points", HexPts, col = 1)), 415 split = c(1, 1, 2, 1), more = TRUE) 416print(spplot(HexPolsDf["dist"], cuts=10, col.regions=grys), 417 split = c(2, 1, 2, 1), more = FALSE) 418dev.null <- dev.off() 419cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="") 420.iwidth <- 5 421.iheight <- 6 422.ipointsize <- 12 423 424 425################################################### 426### code chunk number 43: csdacm.Rnw:624-631 427################################################### 428setClass("SpatialHexGrid", representation("SpatialPoints", dx = "numeric"), 429 validity <- function(object) { 430 if (object@dx <= 0) 431 stop("dx should be positive") 432 TRUE 433 } 434) 435 436 437################################################### 438### code chunk number 44: csdacm.Rnw:633-634 439################################################### 440options("width"=40) 441 442 443################################################### 444### code chunk number 45: csdacm.Rnw:636-644 445################################################### 446setClass("SpatialHexGridDataFrame", representation("SpatialPointsDataFrame", dx = "numeric"), 447# NOT TOO WIDE 448 validity <- function(object) { 449 if (object@dx <= 0) 450 stop("dx should be positive") 451 TRUE 452 } 453) 454 455 456################################################### 457### code chunk number 46: csdacm.Rnw:646-647 458################################################### 459options("width"=70) 460 461 462################################################### 463### code chunk number 47: csdacm.Rnw:664-669 464################################################### 465HexPts <- spsample(meuse.grid, type="hexagonal", cellsize=200) 466Hex <- new("SpatialHexGrid", HexPts, dx = 200) 467df <- over(Hex, meuse.grid) 468spdf <- SpatialPointsDataFrame(HexPts, df) 469HexDf <- new("SpatialHexGridDataFrame", spdf, dx = 200) 470 471 472################################################### 473### code chunk number 48: csdacm.Rnw:676-679 474################################################### 475is(HexDf, "SpatialHexGrid") 476setIs("SpatialHexGridDataFrame", "SpatialHexGrid") 477is(HexDf, "SpatialHexGrid") 478 479 480################################################### 481### code chunk number 49: csdacm.Rnw:689-690 482################################################### 483options("width"=50) 484 485 486################################################### 487### code chunk number 50: csdacm.Rnw:692-701 488################################################### 489# NOT TOO WIDE 490setAs("SpatialHexGrid", "SpatialPolygons", 491 function(from) 492 HexPoints2SpatialPolygons(from, from@dx) 493) 494setAs("SpatialHexGridDataFrame", "SpatialPolygonsDataFrame", 495 function(from) 496 SpatialPolygonsDataFrame(as(obj, "SpatialPolygons"), obj@data, match.ID = FALSE) 497) 498 499 500################################################### 501### code chunk number 51: csdacm.Rnw:703-704 502################################################### 503options("width"=70) 504 505 506################################################### 507### code chunk number 52: csdacm.Rnw:711-724 508################################################### 509setMethod("plot", signature(x = "SpatialHexGrid", y = "missing"), 510 function(x, y, ...) plot(as(x, "SpatialPolygons"), ...) 511) 512setMethod("spplot", signature(obj = "SpatialHexGridDataFrame"), 513 function(obj, ...) 514 spplot(SpatialPolygonsDataFrame( as(obj, "SpatialPolygons"), obj@data, match.ID = FALSE), ...) 515) 516setMethod("spsample", "SpatialHexGrid", function(x, n, type, ...) 517 spsample(as(x, "SpatialPolygons"), n = n, type = type, ...) 518) 519setMethod("over", c("SpatialHexGrid", "SpatialPoints"), function(x, y, ...) 520 over(as(x, "SpatialPolygons"), y) 521) 522 523 524################################################### 525### code chunk number 53: csdacm.Rnw:730-732 (eval = FALSE) 526################################################### 527## spplot(meuse.grid["dist"], sp.layout = list("sp.points", Hex, col = 1)) 528## spplot(HexDf["dist"]) 529 530 531################################################### 532### code chunk number 54: csdacm.Rnw:739-740 (eval = FALSE) 533################################################### 534## as(HexDf, "data.frame") 535 536 537################################################### 538### code chunk number 55: csdacm.Rnw:747-749 539################################################### 540bbox(Hex) 541bbox(as(Hex, "SpatialPolygons")) 542 543 544################################################### 545### code chunk number 56: csdacm.Rnw:770-776 546################################################### 547n <- 10 548x <- data.frame(expand.grid(x1 = 1:n, x2 = 1:n, x3 = 1:n), z = rnorm(n^3)) 549coordinates(x) <- ~x1+x2+x3 550gridded(x) <- TRUE 551fullgrid(x) <- TRUE 552summary(x) 553 554 555################################################### 556### code chunk number 57: csdacm.Rnw:791-792 557################################################### 558options("width"=50) 559 560 561################################################### 562### code chunk number 58: csdacm.Rnw:794-801 563################################################### 564# NOT TOO WIDE 565setClass("SpatialTimeGrid", "SpatialGrid", 566 validity <- function(object) { 567 stopifnot(dimensions(object) == 3) 568 TRUE 569 } 570) 571 572 573################################################### 574### code chunk number 59: csdacm.Rnw:803-804 575################################################### 576options("width"=70) 577 578 579################################################### 580### code chunk number 60: csdacm.Rnw:811-819 581################################################### 582setClass("SpatialTimeGridDataFrame", "SpatialGridDataFrame", 583 validity <- function(object) { 584 stopifnot(dimensions(object) == 3) 585 TRUE 586 } 587) 588setIs("SpatialTimeGridDataFrame", "SpatialTimeGrid") 589x <- new("SpatialTimeGridDataFrame", x) 590 591 592################################################### 593### code chunk number 61: csdacm.Rnw:824-835 594################################################### 595summary.SpatialTimeGridDataFrame <- function(object, ...) { 596 cat("Object of class SpatialTimeGridDataFrame\n") 597 x <- gridparameters(object) 598 t0 <- ISOdate(1970,1,1,0,0,0) 599 t1 <- t0 + x[3,1] 600 cat(paste("first time step:", t1, "\n")) 601 t2 <- t0 + x[3,1] + (x[3,3] - 1) * x[3,2] 602 cat(paste("last time step: ", t2, "\n")) 603 cat(paste("time step: ", x[3,2], "\n")) 604 summary(as(object, "SpatialGridDataFrame")) 605} 606 607 608################################################### 609### code chunk number 62: csdacm.Rnw:837-838 610################################################### 611options("width"=50) 612 613 614################################################### 615### code chunk number 63: csdacm.Rnw:840-843 616################################################### 617# NOT TOO WIDE 618setMethod("summary", "SpatialTimeGridDataFrame", summary.SpatialTimeGridDataFrame) 619summary(x) 620 621 622################################################### 623### code chunk number 64: csdacm.Rnw:845-846 624################################################### 625options("width"=70) 626 627 628################################################### 629### code chunk number 65: csdacm.Rnw:853-880 630################################################### 631subs.SpatialTimeGridDataFrame <- function(x, i, j, ..., drop=FALSE) { 632 t <- coordinates(x)[,3] + ISOdate(1970,1,1,0,0,0) 633 if (missing(j)) 634 j <- TRUE 635 sel <- t %in% i 636 if (! any(sel)) 637 stop("selection results in empty set") 638 fullgrid(x) <- FALSE 639 if (length(i) > 1) { 640 x <- x[i = sel, j = j,...] 641 fullgrid(x) <- TRUE 642 as(x, "SpatialTimeGridDataFrame") 643 } else { 644 gridded(x) <- FALSE 645 x <- x[i = sel, j = j,...] 646 cc <- coordinates(x)[,1:2] 647 p4s <- CRS(proj4string(x)) 648# NOT TOO WIDE 649 SpatialPixelsDataFrame(cc, x@data, proj4string = p4s) 650 } 651} 652setMethod("[", c("SpatialTimeGridDataFrame", "POSIXct", "ANY"), 653 subs.SpatialTimeGridDataFrame) 654t1 <- as.POSIXct("1970-01-01 0:00:03", tz = "GMT") 655t2 <- as.POSIXct("1970-01-01 0:00:05", tz = "GMT") 656summary(x[c(t1,t2)]) 657summary(x[t1]) 658 659 660################################################### 661### code chunk number 66: csdacm.Rnw:893-906 662################################################### 663spplot.stgdf <- function(obj, zcol = 1, ..., format = NULL) { 664# NOT TOO WIDE 665 if (length(zcol) != 1) 666 stop("can only plot a single attribute") 667 if (is.null(format)) format <- "%Y-%m-%d %H:%M:%S" 668 cc <- coordinates(obj) 669 df <- unstack(data.frame(obj[[zcol]], cc[,3])) 670 ns <- as.character(coordinatevalues(getGridTopology(obj))[[3]] + ISOdate(1970,1,1,0,0,0), format = format) 671 cc2d <- cc[cc[,3] == min(cc[,3]), 1:2] 672 obj <- SpatialPixelsDataFrame(cc2d, df) 673 spplot(obj, names.attr = ns,...) 674} 675setMethod("spplot", "SpatialTimeGridDataFrame", spplot.stgdf) 676 677 678################################################### 679### code chunk number 67: csdacm.Rnw:912-918 680################################################### 681.iwidth <- 6 682.iheight <- 4 683.PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="") 684pdf(file=file, width = .iwidth, height = .iheight, pointsize = .ipointsize) 685opar <- par(mar=c(3,3,1,1)+0.1) 686print(spplot(x, format = "%H:%M:%S", as.table=TRUE)) 687dev.null <- dev.off() 688cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="") 689.iwidth <- 5 690.iheight <- 6 691.ipointsize <- 12 692 693 694################################################### 695### code chunk number 68: csdacm.Rnw:927-932 (eval = FALSE) 696################################################### 697## library(lattice) 698## trellis.par.set(canonical.theme(color = FALSE)) 699## spplot(x, format = "%H:%M:%S", as.table=TRUE, cuts=6, 700## col.regions=grey.colors(7, 0.55, 0.95, 2.2)) 701## # RSB quietening greys 702 703 704################################################### 705### code chunk number 69: csdacm.Rnw:938-939 (eval = FALSE) 706################################################### 707## ?as.character.POSIXt 708 709 710################################################### 711### code chunk number 70: csdacm.Rnw:963-969 712################################################### 713library(gstat) 714data(meuse) 715coordinates(meuse) <- ~x+y 716v <- vgm(.5, "Sph", 800, .05) 717sim <- krige(log(zinc)~1, meuse, meuse.grid, v, nsim=100, nmax=30) 718sim@data <- exp(sim@data) 719 720 721################################################### 722### code chunk number 71: csdacm.Rnw:977-981 723################################################### 724quantile.Spatial <- function(x, ..., byLayer = FALSE) { 725 stopifnot("data" %in% slotNames(x)) 726 apply(x@data, ifelse(byLayer, 2, 1), quantile, ...) 727} 728 729 730################################################### 731### code chunk number 72: csdacm.Rnw:987-989 732################################################### 733sim$lower <- quantile.Spatial(sim[1:100], probs = 0.025) 734sim$upper <- quantile.Spatial(sim[1:100], probs = 0.975) 735 736 737################################################### 738### code chunk number 73: csdacm.Rnw:996-997 739################################################### 740medians <- quantile.Spatial(sim[1:100], probs = 0.5, byLayer = TRUE) 741 742 743################################################### 744### code chunk number 74: csdacm.Rnw:999-1000 (eval = FALSE) 745################################################### 746## hist(medians) 747 748 749################################################### 750### code chunk number 75: csdacm.Rnw:1014-1015 751################################################### 752options("width"=50) 753 754 755################################################### 756### code chunk number 76: csdacm.Rnw:1017-1022 757################################################### 758fractionBelow <- function(x, q, byLayer = FALSE) { 759 stopifnot(is(x, "Spatial") || !("data" %in% slotNames(x))) 760 apply(x@data < q, ifelse(byLayer, 2, 1), function(r) sum(r)/length(r)) 761# NOT TOO WIDE 762} 763 764 765################################################### 766### code chunk number 77: csdacm.Rnw:1024-1025 767################################################### 768options("width"=70) 769 770 771################################################### 772### code chunk number 78: csdacm.Rnw:1027-1030 773################################################### 774over500 <- 1 - fractionBelow(sim[1:100], 200, byLayer = TRUE) 775summary(over500) 776quantile(over500, c(0.025, 0.975)) 777 778 779################################################### 780### code chunk number 79: csdacm.Rnw:1046-1047 781################################################### 782run <- require(rgdal, quietly=TRUE) 783 784 785################################################### 786### code chunk number 80: csdacm.Rnw:1051-1054 787################################################### 788if (run) { 789 fn <- system.file("pictures/erdas_spnad83.tif", package = "rgdal")[1] 790} 791 792 793################################################### 794### code chunk number 81: csdacm.Rnw:1056-1059 (eval = FALSE) 795################################################### 796## x <- readGDAL(fn, output.dim = c(120, 132)) 797## x$band1[x$band1 <= 0] <- NA 798## spplot(x, col.regions=bpy.colors()) 799 800 801################################################### 802### code chunk number 82: csdacm.Rnw:1070-1082 803################################################### 804if (run) { 805library(rgdal) 806x <- GDAL.open(fn) 807class(x) 808} 809if (run) { 810x.subs <- x[1:100, 1:100, 1] 811class(x.subs) 812} 813if (run) { 814gridparameters(x.subs) 815} 816 817 818################################################### 819### code chunk number 83: csdacm.Rnw:1100-1101 820################################################### 821options("width"=50) 822 823 824################################################### 825### code chunk number 84: csdacm.Rnw:1103-1110 826################################################### 827if (run) { 828setClass("SpatialGDAL", 829 representation("Spatial", grid = "GridTopology", grod = "GDALReadOnlyDataset", 830# NOT TOO WIDE 831 name = "character")) 832setClass("SpatialGDALWrite", "SpatialGDAL") 833} 834 835 836################################################### 837### code chunk number 85: csdacm.Rnw:1112-1113 838################################################### 839options("width"=70) 840 841 842################################################### 843### code chunk number 86: csdacm.Rnw:1124-1140 (eval = FALSE) 844################################################### 845## x <- open.SpatialGDAL(fn) 846## nrows <- GDALinfo(fn)["rows"] 847## ncols <- GDALinfo(fn)["columns"] 848## xout <- copy.SpatialGDAL(x, "erdas_spnad83_out.tif") 849## bls <- 20 850## for (i in 1:(nrows/bls - 1)) { 851## r <- 1+(i-1)*bls 852## for (j in 1:(ncols/bls - 1)) { 853## c <- 1+(j-1)*bls 854## x.in <- x[r:(r+bls),c:(c+bls)] 855## xout[r:(r+bls),c:(c+bls)] <- x.in$band1 + 10 #$ 856## } 857## cat(paste("row-block", i, "\n")) 858## } 859## close(x) 860## close(xout) 861 862 863################################################### 864### code chunk number 87: csdacm.Rnw:1147-1154 (eval = FALSE) 865################################################### 866## setMethod("[", "SpatialGDAL", 867## function(x, i, j, ... , drop = FALSE) 868## x@grod[i = i, j = j, ...] 869## ) 870## setReplaceMethod("[", "SpatialGDALWrite", function(x, i, j, ..., value) { 871## ... 872## }) 873 874 875################################################### 876### code chunk number 88: csdacm.Rnw:1180-1181 877################################################### 878options("width"=owidth) 879 880 881