R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("GRASS-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('GRASS') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "contour.G" > > ### * contour.G > > flush(stderr()); flush(stdout()) > > ### Name: contour.G > ### Title: Equal scale contour plots for GRASS raster and site data > ### Aliases: contour.G > ### Keywords: spatial IO > > ### ** Examples > > data(utm.maas) > Zn.o <- as.ordered(cut(utm.maas$Zn, labels=c("insignificant", "low", + "medium", "high", "crisis"), breaks=c(100, 200, 400, 700, 1000, 2000), + include.lowest=TRUE)) > G <- maas.metadata > contour.G(G) > points(utm.maas$east, utm.maas$north, pch=unclass(Zn.o)) > legend(x=c(269800, 270500), y=c(5652300, 5653000), pch=c(1:5), + legend=levels(Zn.o)) > title("Note equal east and north scales") > contour.G(G, kde2d.G(utm.maas$east, utm.maas$north, h=c(300,300), G=G, Z=utm.maas$Zn)*maasmask) > points(utm.maas$east, utm.maas$north, pch=unclass(Zn.o)) > legend(x=c(269800, 270500), y=c(5652300, 5653000), pch=c(1:5), + legend=levels(Zn.o)) > title("Kernel density representation of soil Zn") > > > > cleanEx(); ..nameEx <- "gmeta" > > ### * gmeta > > flush(stderr()); flush(stdout()) > > ### Name: gmeta > ### Title: Reads GRASS metadata from the current LOCATION > ### Aliases: gmeta east north east.default north.default east.grassmeta > ### north.grassmeta obsno reverse make.maas.location > ### Keywords: spatial IO > > ### ** Examples > > if(!get("maas.loc", env = .GRASS.meta)) make.maas.location() [1] TRUE > G <- gmeta() > summary(G) Data from GRASS 5.0 LOCATION maas with 259 columns and 232 rows; UTM, zone: 32 The west-east range is: 269870, 272460, and the south-north: 5650610, 5652930; West-east cell sizes are 10 units, and south-north 10 units. > > > > cleanEx(); ..nameEx <- "interp.new.G" > > ### * interp.new.G > > flush(stderr()); flush(stdout()) > > ### Name: interp.new.G > ### Title: Interpolation from sites to raster, Akima's new method > ### Aliases: interp.new.G > ### Keywords: spatial IO > > ### ** Examples > > data(utm.maas) > Zn.o <- as.ordered(cut(utm.maas$Zn, labels=c("insignificant", "low", + "medium", "high", "crisis"), breaks=c(100, 200, 400, 700, 1000, 2000), + include.lowest=TRUE)) > G <- maas.metadata > pdata <- cbind(utm.maas$elev, utm.maas$d.river, log(utm.maas$d.river), + utm.maas$Zn, log(utm.maas$Zn)) > colnames(pdata) <- c("Elevation", "Distance", "Ln(Distance)", "Zinc", + "Ln(Zinc)") > pairs(pdata) > mod1 <- lm(log(Zn) ~ elev + log(d.river), data=utm.maas) > summary(mod1) Call: lm(formula = log(Zn) ~ elev + log(d.river), data = utm.maas) Residuals: Min 1Q Median 3Q Max -0.87849 -0.26956 0.01677 0.26452 0.85981 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 10.00167 0.27394 36.510 < 2e-16 *** elev -0.29266 0.03718 -7.872 5.58e-12 *** log(d.river) -0.33314 0.03542 -9.407 3.07e-15 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3944 on 95 degrees of freedom Multiple R-Squared: 0.7486, Adjusted R-squared: 0.7433 F-statistic: 141.4 on 2 and 95 DF, p-value: < 2.2e-16 > anova(mod1) Analysis of Variance Table Response: log(Zn) Df Sum Sq Mean Sq F value Pr(>F) elev 1 30.2237 30.2237 194.339 < 2.2e-16 *** log(d.river) 1 13.7610 13.7610 88.483 3.067e-15 *** Residuals 95 14.7745 0.1555 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > inregion <- (utm.maas$east >= G$w & utm.maas$east <= G$e) & + (utm.maas$north >= G$s & utm.maas$north <= G$n) > if(all(!inregion)) stop("None of the site locations are inside the current GRASS region") > if(any(!inregion)) + warning("Some site locations are outside the current GRASS region") > require(akima) Loading required package: akima [1] TRUE > elev.grid <- interp.new.G(G, utm.maas$east, utm.maas$north, utm.maas$elev, + extrap=TRUE) Warning in if (extrap & linear) warning("Cannot extrapolate with linear option") : the condition has length > 1 and only the first element will be used > brks <- c(-Inf, seq(5.4, 10.6, 0.52), +Inf) > plot(G, elev.grid*maasmask, breaks=brks, + col=c("yellow", grey(11:2/11), "red")) > plot(G, ifelse(is.na(maasmask), 1, NA), breaks=c(0,2), col="wheat", add=TRUE) > title("Bicubic spline interpolation: local elevation") > legend(c(269900, 270600), c(5652000, 5652900), bty="n", bg="wheat", + legend=c(rev(leglabs(brks)), "mask=NA"), + fill=c("red", rev(grey(11:2/11)), "yellow", "wheat")) > ldist <- loess(d.river ~ east * north, data=utm.maas, span=0.2, + control = loess.control(surface = "direct")) > bank <- predict(ldist, newdata=data.frame(east=east(maas.metadata), + north=north(maas.metadata))) > b1 <- bank*maasmask > b1[b1 < 5] <- 5 > brks <- c(seq(0, 1000, 200), +Inf) > cols <- grey(6:2/6) > plot(maas.metadata, b1, breaks=brks, col=c(cols, "red")) > plot(maas.metadata, ifelse(is.na(maasmask), 1, NA), breaks=c(0,2), col="wheat", add=TRUE) > title("Loess predictions of distance from river bank") > legend(c(269900, 270600), c(5652000, 5652900), bty="n", bg="wheat", + legend=c(rev(leglabs(brks)), "mask=NA"), + fill=c("red", rev(cols), "wheat")) > new <- data.frame(elev=elev.grid*maasmask, d.river=b1) > pr.mod1 <- predict(mod1, new, se.fit=TRUE) > rm(new, elev.grid, bank) > v.pr <- rep(NA, G$Ncells) > v.pr[as.integer(names(pr.mod1$fit))] <- pr.mod1$fit > summary(v.pr) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 3.724 4.980 5.386 5.419 5.819 7.575 17424.000 > summary(exp(v.pr)) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 41.42 145.40 218.40 286.70 336.80 1948.00 17424.00 > plot(G, exp(v.pr), breaks=c(-200, round(seq(142, 1154, length=9)), 5000), + col=c("yellow", grey(9:2/9), "red")) > plot(G, ifelse(is.na(maasmask), 1, NA), col="wheat", add=TRUE) > title("Regression predictions of Zn levels, (B&McD p. 113)") > legend(c(269900, 270600), c(5652000, 5652900), bty="n", bg="wheat", + legend=c("> 1154", rev(legtext(round(seq(142, 1154, length=9)))), "< 142", + "mask=NA"), fill=c("red", rev(grey(9:2/9)), "yellow", "wheat")) > v.pr[as.integer(names(pr.mod1$se.fit))] <- pr.mod1$se.fit > plot(G, exp(v.pr), col=grey(16:2/16)) > plot(G, ifelse(is.na(maasmask), 1, NA), col="wheat", add=TRUE) > points(utm.maas$east, utm.maas$north, pch=18) > title("Standard error of predictions") > Zn.grid <- interp.new.G(G, utm.maas$east, utm.maas$north, utm.maas$Zn, + extrap=TRUE) Warning in if (extrap & linear) warning("Cannot extrapolate with linear option") : the condition has length > 1 and only the first element will be used > plot(G, Zn.grid*maasmask, breaks=c(-500, round(seq(15,1994,length=9)), 3000), + col=c("yellow", grey(9:2/9), "red")) > plot(G, ifelse(is.na(maasmask), 1, NA), col="wheat", add=TRUE) > title(xlab="B&McD p. 118", main="Bicubic spline interpolation: Zn ppm") > legend(c(269900, 270600), c(5652000, 5652900), bty="n", bg="wheat", + legend=c("> 1994", rev(legtext(round(seq(15,1994,length=9)))), "< 15", + "mask=NA"), fill=c("red", rev(grey(9:2/9)), "yellow", "wheat")) > ## Not run: > ##D sites.put2(G, utm.maas, dims=c("east", "north", "elev"), > ##D lname="ex.utm.maas1", check=FALSE) > ##D system("s.surf.idw input=ex.utm.maas1 output=Zn.idw npoints=98 field=7") > ##D idw <- rast.get(G, "Zn.idw") > ##D plot(G, idw$Zn.idw*maasmask, breaks=c(-100, round(seq(15,1994,length=9)), 3000), > ##D col=c("yellow", grey(9:2/9), "red")) > ##D plot(G, ifelse(is.na(maasmask), 1, NA), col="wheat", add=TRUE) > ##D title(xlab="B&McD p. 118", > ##D main="Inverse squared distance interpolation: Zn ppm") > ##D legend(c(269900, 270600), c(5652000, 5652900), bty="n", bg="wheat", > ##D legend=c("> 1994", rev(legtext(round(seq(15,1994,length=9)))), "< 15", > ##D "mask=NA"), fill=c("red", rev(grey(9:2/9)), "yellow", "wheat")) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "kde2d.G" > > ### * kde2d.G > > flush(stderr()); flush(stdout()) > > ### Name: kde2d.G > ### Title: Two-Dimensional Kernel Density Estimation on a GRASS Grid > ### Aliases: kde2d.G > ### Keywords: spatial IO > > ### ** Examples > > data(utm.maas) > G <- maas.metadata > inregion <- (utm.maas$east >= G$w & utm.maas$east <= G$e) & + (utm.maas$north >= G$s & utm.maas$north <= G$n) > if(all(!inregion)) stop("None of the site locations are inside the current GRASS region") > if(any(!inregion)) warning("Some site locations are outside the current GRASS region") > plot(G, kde2d.G(G=G, utm.maas$east, utm.maas$north, h=c(300,300))*maasmask) > points(utm.maas$east, utm.maas$north) > rug(utm.maas$east, side=1, ticksize=0.02) > rug(utm.maas$north, side=2, ticksize=0.02) > title(main="Kernel estimates of density of soil sample placing", + xlab="(Bailey & Gatrell, pp. 84-88") > plot(G, kde2d.G(G=G, utm.maas$east, utm.maas$north, h=c(300,300), Z=utm.maas$Zn)*maasmask) > points(utm.maas$east, utm.maas$north) > rug(utm.maas$east, side=1, ticksize=0.02) > rug(utm.maas$north, side=2, ticksize=0.02) > title(main="Kernel density weighted average, Zn ppm", + xlab="(Bailey & Gatrell, pp. 159-161") > > > > cleanEx(); ..nameEx <- "krige.G" > > ### * krige.G > > flush(stderr()); flush(stdout()) > > ### Name: krige.G > ### Title: Kriging prediction for GRASS maps > ### Aliases: krige.G prmat2.G semat2.G > ### Keywords: spatial IO > > ### ** Examples > > data(utm.maas) > G <- maas.metadata > inregion <- (utm.maas$east >= G$w & utm.maas$east <= G$e) & + (utm.maas$north >= G$s & utm.maas$north <= G$n) > if(all(!inregion)) + stop("None of the site locations are inside the current GRASS region") > if(any(!inregion)) + warning("Some site locations are outside the current GRASS region") > require(sgeostat) Loading required package: sgeostat [1] TRUE > maas.pts <- point(utm.maas, x="east", y="north") > maas.pairs <- pair(maas.pts, num.lags=10, maxdist=1000) ................................................................................................. > maas.evg <- est.variogram(maas.pts, maas.pairs, a1="Zn") > plot(maas.evg) > text(maas.evg$bins, maas.evg$classic/2, labels=maas.evg$n, pos=1) > abline(h=25000) > abline(h=170000) > maas.fit.exp <- fit.exponential(maas.evg, c0=24803.4, ce=158232, + ae=361.604, plot.it=TRUE, iterations=0) Convergence not achieved! > cat("Using parameters from B&McD, p. 142 without fitting.\n") Using parameters from B&McD, p. 142 without fitting. > res.G1 <- krige.G(maas.pts, "Zn", maas.fit.exp, G, maasmask) Loading required package: spatial > summary(res.G1$zhat) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 127.6 226.0 352.4 403.9 524.0 1655.0 17424.0 > summary(res.G1$sehat) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 188.2 243.4 275.1 292.9 341.5 416.0 17424.0 > plot(G, res.G1$zhat, breaks=round(seq(120,1660,length=9)), col=grey(9:2/9)) > plot(G, ifelse(is.na(maasmask), 1, NA), breaks=c(0,2), col="wheat", add=TRUE) > legend(c(269900, 270600), c(5652000, 5652900), bty="n", bg="wheat", + legend=c(rev(legtext(round(seq(120,1660,length=9)))), "mask=NA"), + fill=c(rev(grey(9:2/9)), "wheat")) > title("Ordinary kriging predictions") > > > > cleanEx(); ..nameEx <- "pcbs" > > ### * pcbs > > flush(stderr()); flush(stdout()) > > ### Name: pcbs > ### Title: PCBs in an area of South Wales > ### Aliases: pcbs > ### Keywords: datasets > > ### ** Examples > > data(pcbs) > pcbs.o <- as.ordered(cut(pcbs$pcbs, labels=c("insignificant", "low", + "medium", "high", "crisis"), breaks=c(1,20,100,500,1000,5000), + include.lowest=TRUE)) > table(pcbs.o) pcbs.o insignificant low medium crisis 53 14 1 2 > plot(pcbs$east, pcbs$north, pch=unclass(pcbs.o), xlab="", ylab="", asp=1) > legend(x=c(67980, 68480), y=c(74710, 75180), pch=c(1:5), legend=levels(pcbs.o)) > > > > > cleanEx(); ..nameEx <- "plot.grassmeta" > > ### * plot.grassmeta > > flush(stderr()); flush(stdout()) > > ### Name: plot.grassmeta > ### Title: Equal scale plots for GRASS raster and site data > ### Aliases: plot.grassmeta legtext leglabs > ### Keywords: spatial IO hplot > > ### ** Examples > > data(utm.maas) > Zn.o <- as.ordered(cut(utm.maas$Zn, labels=c("insignificant", "low", + "medium", "high", "crisis"), breaks=c(100, 200, 400, 700, 1000, 2000), + include.lowest=TRUE)) > G <- maas.metadata > plot(G) > points(utm.maas$east, utm.maas$north, pch=unclass(Zn.o)) > legend(x=c(269800, 270500), y=c(5652300, 5653000), pch=c(1:5), + legend=levels(Zn.o)) > title("Note equal east and north scales") > example(kde2d.G) kd2d.G> data(utm.maas) kd2d.G> G <- maas.metadata kd2d.G> inregion <- (utm.maas$east >= G$w & utm.maas$east <= G$e) & (utm.maas$north >= G$s & utm.maas$north <= G$n) kd2d.G> if (all(!inregion)) stop("None of the site locations are inside the current GRASS region") kd2d.G> if (any(!inregion)) warning("Some site locations are outside the current GRASS region") kd2d.G> plot(G, kde2d.G(G = G, utm.maas$east, utm.maas$north, h = c(300, 300)) * maasmask) kd2d.G> points(utm.maas$east, utm.maas$north) kd2d.G> rug(utm.maas$east, side = 1, ticksize = 0.02) kd2d.G> rug(utm.maas$north, side = 2, ticksize = 0.02) kd2d.G> title(main = "Kernel estimates of density of soil sample placing", xlab = "(Bailey & Gatrell, pp. 84-88") kd2d.G> plot(G, kde2d.G(G = G, utm.maas$east, utm.maas$north, h = c(300, 300), Z = utm.maas$Zn) * maasmask) kd2d.G> points(utm.maas$east, utm.maas$north) kd2d.G> rug(utm.maas$east, side = 1, ticksize = 0.02) kd2d.G> rug(utm.maas$north, side = 2, ticksize = 0.02) kd2d.G> title(main = "Kernel density weighted average, Zn ppm", xlab = "(Bailey & Gatrell, pp. 159-161") > > > > cleanEx(); ..nameEx <- "rast.get" > > ### * rast.get > > flush(stderr()); flush(stdout()) > > ### Name: rast.get > ### Title: Import GRASS raster files > ### Aliases: rast.get > ### Keywords: spatial IO > > ### ** Examples > > if(!get("maas.loc", env = .GRASS.meta)) make.maas.location() > data(utm.maas) > Zn.o <- as.ordered(cut(utm.maas$Zn, labels=c("insignificant", "low", + "medium", "high", "crisis"), breaks=c(100, 200, 400, 700, 1000, 2000), + include.lowest=TRUE)) > G <- gmeta() > if(length(ls(pat="nameR"))==0){example(rast.put)} rst.pt> if (!get("maas.loc", env = .GRASS.meta)) make.maas.location() rst.pt> data(utm.maas) rst.pt> G <- gmeta() rst.pt> inregion <- (utm.maas$east >= G$w & utm.maas$east <= G$e) & (utm.maas$north >= G$s & utm.maas$north <= G$n) rst.pt> if (all(!inregion)) stop("None of the site locations are inside the current GRASS region") rst.pt> if (any(!inregion)) warning("Some site locations are outside the current GRASS region") rst.pt> require(akima) Loading required package: akima [1] TRUE rst.pt> Zn.grid <- interp.new.G(G, utm.maas$east, utm.maas$north, utm.maas$Zn, extrap = TRUE) * maasmask Warning in if (extrap & linear) warning("Cannot extrapolate with linear option") : the condition has length > 1 and only the first element will be used rst.pt> require(spatial) Loading required package: spatial [1] TRUE rst.pt> s3 <- trmat.G(G, surf.ls(3, utm.maas$east, utm.maas$north, utm.maas$Zn)) * maasmask rst.pt> Zn.grid.o <- as.ordered(cut(Zn.grid, labels = c("insignificant", "low", "medium", "high", "crisis"), breaks = c(-300, 200, 400, 700, 1000, 5000), include.lowest = TRUE)) rst.pt> nameR <- c("ex.Zn.grid.in", "ex.tr3.in", "ex.Zn.grid.o.in") rst.pt> rast.put(G, lname = nameR[1], layer = Zn.grid, title = "Akima spline interpolation", check = FALSE) rst.pt> rast.put(G, lname = nameR[2], layer = s3, title = "Cubic trend surface", check = FALSE) rst.pt> rast.put(G, lname = nameR[3], layer = Zn.grid.o, cat = TRUE, title = "Interpolated surface categories", check = FALSE) > exget <- rast.get(G, rlist=c("ex.Zn.grid.in", "ex.tr3.in", + "ex.Zn.grid.o.in"), catlabels=as.logical(c(FALSE, FALSE, TRUE))) > oldpar <- par(mfrow=c(1,2)) > plot(G, exget$ex.Zn.grid.in) > points(utm.maas$east, utm.maas$north, pch=18) > title("Bicubic spline interpolation") > plot(G, exget$ex.tr3.in) > points(utm.maas$east, utm.maas$north, pch=18) > title("Cubic trend surface") > par(oldpar) > table(exget$ex.Zn.grid.o.in) insignificant low medium high crisis 18740 12694 6431 2649 2150 > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "rast.put" > > ### * rast.put > > flush(stderr()); flush(stdout()) > > ### Name: rast.put > ### Title: Exports to GRASS raster data files > ### Aliases: rast.put > ### Keywords: spatial IO > > ### ** Examples > > if(!get("maas.loc", env = .GRASS.meta)) make.maas.location() > data(utm.maas) > G <- gmeta() > inregion <- (utm.maas$east >= G$w & utm.maas$east <= G$e) & (utm.maas$north >= G$s & utm.maas$north <= G$n) > if(all(!inregion)) stop("None of the site locations are inside the current GRASS region") > if(any(!inregion)) warning("Some site locations are outside the current GRASS region") > require(akima) Loading required package: akima [1] TRUE > Zn.grid <- interp.new.G(G, utm.maas$east, utm.maas$north, utm.maas$Zn, + extrap=TRUE)*maasmask Warning in if (extrap & linear) warning("Cannot extrapolate with linear option") : the condition has length > 1 and only the first element will be used > require(spatial) Loading required package: spatial [1] TRUE > s3 <- trmat.G(G, surf.ls(3, utm.maas$east, utm.maas$north, + utm.maas$Zn))*maasmask > Zn.grid.o <- as.ordered(cut(Zn.grid, labels=c("insignificant", "low", + "medium", "high", "crisis"), breaks=c(-300, 200, 400, 700, 1000, 5000), include.lowest=TRUE)) > nameR <- c("ex.Zn.grid.in", "ex.tr3.in", "ex.Zn.grid.o.in") > rast.put(G, lname=nameR[1], layer=Zn.grid, title="Akima spline interpolation", + check=FALSE) > rast.put(G, lname=nameR[2], layer=s3, title="Cubic trend surface", + check=FALSE) > rast.put(G, lname=nameR[3], layer=Zn.grid.o, cat=TRUE, + title="Interpolated surface categories", check=FALSE) > > > > cleanEx(); ..nameEx <- "sites.get" > > ### * sites.get > > flush(stderr()); flush(stdout()) > > ### Name: sites.get > ### Title: Import GRASS sites file > ### Aliases: sites.get > ### Keywords: spatial IO > > ### ** Examples > > if(!get("maas.loc", env = .GRASS.meta)) make.maas.location() > G <- gmeta() > example(sites.put) sts.pt> if (!get("maas.loc", env = .GRASS.meta)) make.maas.location() sts.pt> data(utm.maas) sts.pt> Zn.o <- as.ordered(cut(utm.maas$Zn, labels = c("insignificant", "low", "medium", "high", "crisis"), breaks = c(100, 200, 400, 700, 1000, 2000), include.lowest = TRUE)) sts.pt> G <- gmeta() sts.pt> nameQ <- c("ex.Zn.in", "ex.Znlog.in", "ex.Zncat.in") sts.pt> Zn <- data.frame(east = utm.maas$east, north = utm.maas$north, Zn = utm.maas$Zn) sts.pt> sites.put2(G, Zn, dims = c("east", "north"), lname = nameQ[1], check = FALSE) sts.pt> log.Zn <- data.frame(east = utm.maas$east, north = utm.maas$north, log.Zn = log(utm.maas$Zn)) sts.pt> sites.put2(G, log.Zn, dims = c("east", "north"), lname = nameQ[2], check = FALSE) sts.pt> Zn.o.df <- data.frame(east = utm.maas$east, north = utm.maas$north, Zn.o = Zn.o) sts.pt> sites.put2(G, Zn.o.df, dims = c("east", "north"), lname = nameQ[3], check = FALSE) sts.pt> newdf <- data.frame(utm.maas, Zn.o) sts.pt> sites.put2(G, newdf, dims = c("east", "north", "elev"), lname = "ex.utm.maas", check = FALSE) > ex.Zn.log <- sites.get(G, slist="ex.Znlog.in") > mean(ex.Zn.log$log.Zn - log(utm.maas$Zn)) [1] -4.80406e-10 > ex.Zn.cat <- sites.get(G, slist="ex.Zncat.in") > table(Zn.o, ex.Zn.cat$Zn.o) Zn.o crisis high insignificant low medium insignificant 0 0 31 0 0 low 0 0 0 24 0 medium 0 0 0 0 18 high 0 15 0 0 0 crisis 10 0 0 0 0 > utm.maas.new <- sites.get(G, slist="ex.utm.maas") > names(utm.maas.new) [1] "east" "north" "elev" "id" "x" "y" "d.river" [8] "Cd" "Cu" "Pb" "Zn" "LOI" "Fldf" "Soil" [15] "Zn.o" > table(Zn.o, utm.maas.new$Zn.o) Zn.o crisis high insignificant low medium insignificant 0 0 31 0 0 low 0 0 0 24 0 medium 0 0 0 0 18 high 0 15 0 0 0 crisis 10 0 0 0 0 > > > > cleanEx(); ..nameEx <- "sites.put" > > ### * sites.put > > flush(stderr()); flush(stdout()) > > ### Name: sites.put > ### Title: Export GRASS sites file > ### Aliases: sites.put sites.put2 > ### Keywords: spatial IO > > ### ** Examples > > if(!get("maas.loc", env = .GRASS.meta)) make.maas.location() > data(utm.maas) > Zn.o <- as.ordered(cut(utm.maas$Zn, labels=c("insignificant", "low", + "medium", "high", "crisis"), breaks=c(100, 200, 400, 700, 1000, 2000), + include.lowest=TRUE)) > G <- gmeta() > nameQ <- c("ex.Zn.in", "ex.Znlog.in", "ex.Zncat.in") > Zn <- data.frame(east=utm.maas$east, north=utm.maas$north, Zn=utm.maas$Zn) > sites.put2(G, Zn, dims=c("east", "north"), lname=nameQ[1], check=FALSE) > log.Zn <- data.frame(east=utm.maas$east, north=utm.maas$north, + log.Zn=log(utm.maas$Zn)) > sites.put2(G, log.Zn, dims=c("east", "north"), lname=nameQ[2], check=FALSE) > Zn.o.df <- data.frame(east=utm.maas$east, north=utm.maas$north, Zn.o=Zn.o) > sites.put2(G, Zn.o.df, dims=c("east", "north"), lname=nameQ[3], check=FALSE) > newdf <- data.frame(utm.maas, Zn.o) > sites.put2(G, newdf, dims=c("east", "north", "elev"), lname="ex.utm.maas", + check=FALSE) > > > > cleanEx(); ..nameEx <- "summary.grassmeta" > > ### * summary.grassmeta > > flush(stderr()); flush(stdout()) > > ### Name: summary.grassmeta > ### Title: Display GRASS metadata > ### Aliases: summary.grassmeta > ### Keywords: spatial IO > > ### ** Examples > > if(!get("maas.loc", env = .GRASS.meta)) make.maas.location() > G <- gmeta() > summary(G) Data from GRASS 5.0 LOCATION maas with 259 columns and 232 rows; UTM, zone: 32 The west-east range is: 269870, 272460, and the south-north: 5650610, 5652930; West-east cell sizes are 10 units, and south-north 10 units. > > > > cleanEx(); ..nameEx <- "trmat.G" > > ### * trmat.G > > flush(stderr()); flush(stdout()) > > ### Name: trmat.G > ### Title: Evaluate Trend Surface over a GRASS Grid > ### Aliases: trmat.G > ### Keywords: spatial IO > > ### ** Examples > > data(utm.maas) > G <- maas.metadata > if(G$Ncells != length(maasmask)) + stop("example not running in Maas GRASS location") > require(spatial) Loading required package: spatial [1] TRUE > s1 <- surf.ls(1, utm.maas$east, utm.maas$north, utm.maas$Zn) > s2 <- surf.ls(2, utm.maas$east, utm.maas$north, utm.maas$Zn) > s3 <- surf.ls(3, utm.maas$east, utm.maas$north, utm.maas$Zn) > s4 <- surf.ls(4, utm.maas$east, utm.maas$north, utm.maas$Zn) > summary(s2) Analysis of Variance Table Model: surf.ls(np = 2, x = utm.maas$east, y = utm.maas$north, z = utm.maas$Zn) Sum Sq Df Mean Sq F value Pr(>F) Regression 7320803 5 1464160.60 16.616 1.1560e-11 Deviation 8106812 92 88117.52 Total 15427615 97 Multiple R-Squared: 0.4745, Adjusted R-squared: 0.446 AIC: (df = 92) 1293.678 Fitted: Min 1Q Median 3Q Max 110.2 276.5 433.3 644.6 1231.7 Residuals: Min 1Q Median 3Q Max -706.60 -169.80 -51.62 118.34 881.91 > summary(s3) Analysis of Variance Table Model: surf.ls(np = 3, x = utm.maas$east, y = utm.maas$north, z = utm.maas$Zn) Sum Sq Df Mean Sq F value Pr(>F) Regression 8639815 9 959979.48 12.44559 1.6536e-12 Deviation 6787800 88 77134.09 Total 15427615 97 Multiple R-Squared: 0.56, Adjusted R-squared: 0.515 AIC: (df = 88) 1268.276 Fitted: Min 1Q Median 3Q Max -59.95 245.35 434.17 713.54 1136.49 Residuals: Min 1Q Median 3Q Max -555.52 -199.00 -24.40 136.78 798.99 > summary(s4) Analysis of Variance Table Model: surf.ls(np = 4, x = utm.maas$east, y = utm.maas$north, z = utm.maas$Zn) Sum Sq Df Mean Sq F value Pr(>F) Regression 10596769 14 756912.10 13.00470 1.5543e-15 Deviation 4830846 83 58202.96 Total 15427615 97 Multiple R-Squared: 0.6869, Adjusted R-squared: 0.634 AIC: (df = 83) 1224.945 Fitted: Min 1Q Median 3Q Max 29.07 201.36 408.27 710.75 1396.72 Residuals: Min 1Q Median 3Q Max -421.16 -135.54 -24.69 104.58 742.15 > anova(s2,s3,s4) Analysis of Variance Table Model 1: surf.ls(np = 2, x = utm.maas$east, y = utm.maas$north, z = utm.maas$Zn) Model 2: surf.ls(np = 3, x = utm.maas$east, y = utm.maas$north, z = utm.maas$Zn) Model 3: surf.ls(np = 4, x = utm.maas$east, y = utm.maas$north, z = utm.maas$Zn) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 92 8106812 2 88 6787800 4 1319012 4.2751 0.003298 3 83 4830846 5 1956954 6.7246 2.63e-05 > plot(G, trmat.G(G, s3)*maasmask, breaks=c(-10000, seq(-100, 1500, 200), + 10000), col=c("yellow", grey(8:1/8), "red")) > plot(G, ifelse(is.na(maasmask), 1, NA), breaks=c(0,2), col="wheat", add=TRUE) > plot(s3, add=TRUE) > legend(c(269900, 270600), c(5652000, 5652900), bty="n", bg="wheat", + legend=c("> 1500", rev(legtext(seq(-100,1500,200))), "< -100", "mask=NA"), + fill=c("red", rev(grey(8:1/8)), "yellow", "wheat")) > oldpar <- par(mfrow=c(2,2)) > plot(G, trmat.G(G, s1)*maasmask, breaks=c(-10000, seq(-100, 1500, 200), + 10000), col=c("yellow", grey(8:1/8), "red")) > plot(G, ifelse(is.na(maasmask), 1, NA), breaks=c(0,2), col="wheat", add=TRUE) > plot(s1, add=TRUE) > plot(G, trmat.G(G, s2)*maasmask, breaks=c(-10000, seq(-100, 1500, 200), + 10000), col=c("yellow", grey(8:1/8), "red")) > plot(G, ifelse(is.na(maasmask), 1, NA), breaks=c(0,2), col="wheat", add=TRUE) > plot(s2, add=TRUE) > plot(G, trmat.G(G, s3)*maasmask, breaks=c(-10000, seq(-100, 1500, 200), + 10000), col=c("yellow", grey(8:1/8), "red")) > plot(G, ifelse(is.na(maasmask), 1, NA), breaks=c(0,2), col="wheat", add=TRUE) > plot(s3, add=TRUE) > plot(G, trmat.G(G, s4)*maasmask, breaks=c(-10000, seq(-100, 1500, 200), + 10000), col=c("yellow", grey(8:1/8), "red")) > plot(G, ifelse(is.na(maasmask), 1, NA), breaks=c(0,2), col="wheat", add=TRUE) > plot(s4, add=TRUE) > par(oldpar) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "utm.maas" > > ### * utm.maas > > flush(stderr()); flush(stdout()) > > ### Name: utm.maas > ### Title: Soil pollution data set > ### Aliases: utm.maas maasmask maas.metadata > ### Keywords: datasets > > ### ** Examples > > data(utm.maas) > Zn.o <- as.ordered(cut(utm.maas$Zn, labels=c("insignificant", "low", + "medium", "high", "crisis"), breaks=c(100, 200, 400, 700, 1000, 2000), + include.lowest=TRUE)) > table(Zn.o) Zn.o insignificant low medium high crisis 31 24 18 15 10 > plot(utm.maas$east, utm.maas$north, pch=unclass(Zn.o), xlab="", ylab="", asp=1) > legend(x=c(269900, 270500), y=c(5652000, 5652700), pch=c(1:5), + legend=levels(Zn.o)) > cat("Burrough & McDonnell p. 107\n") Burrough & McDonnell p. 107 > round(tapply(utm.maas$Zn, as.factor(utm.maas$Fldf), mean), 2) 1 2 3 769.77 264.98 205.78 > round(tapply(utm.maas$Zn, as.factor(utm.maas$Fldf), sd), 2) 1 2 3 423.17 176.62 105.33 > round(tapply(log(utm.maas$Zn), as.factor(utm.maas$Fldf), mean), 3) 1 2 3 6.484 5.421 5.239 > round(tapply(log(utm.maas$Zn), as.factor(utm.maas$Fldf), sd), 3) 1 2 3 0.609 0.531 0.416 > round(exp(round(tapply(log(utm.maas$Zn), as.factor(utm.maas$Fldf), mean), 3)), 2) 1 2 3 654.58 226.11 188.48 > hist(utm.maas$Zn, breaks=seq(0,2000,100), col="grey") > hist(log(utm.maas$Zn), breaks=seq(3.5,8.5,0.25), col="grey") > t.test(utm.maas$Zn[utm.maas$Fldf == 2], utm.maas$Zn[utm.maas$Fldf == 3]) Welch Two Sample t-test data: utm.maas$Zn[utm.maas$Fldf == 2] and utm.maas$Zn[utm.maas$Fldf == 3] t = 1.3543, df = 18.243, p-value = 0.1922 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -32.54811 150.94907 sample estimates: mean of x mean of y 264.9783 205.7778 > cat("NB: B&McD p. 108, their relative variance is (1-(RSquared))\n") NB: B&McD p. 108, their relative variance is (1-(RSquared)) > anova(lm(Zn ~ as.factor(Fldf), data=utm.maas)) Analysis of Variance Table Response: Zn Df Sum Sq Mean Sq F value Pr(>F) as.factor(Fldf) 2 6413959 3206979 33.8 8.196e-12 *** Residuals 95 9013656 94881 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > summary(lm(Zn ~ as.factor(Fldf), data=utm.maas)) Call: lm(formula = Zn ~ as.factor(Fldf), data = utm.maas) Residuals: Min 1Q Median 3Q Max -650.77 -125.73 -63.38 49.57 1069.23 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 769.77 46.97 16.387 < 2e-16 *** as.factor(Fldf)2 -504.79 65.34 -7.726 1.13e-11 *** as.factor(Fldf)3 -563.99 112.91 -4.995 2.67e-06 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 308 on 95 degrees of freedom Multiple R-Squared: 0.4157, Adjusted R-squared: 0.4034 F-statistic: 33.8 on 2 and 95 DF, p-value: 8.196e-12 > anova(lm(log(Zn) ~ as.factor(Fldf), data=utm.maas)) Analysis of Variance Table Response: log(Zn) Df Sum Sq Mean Sq F value Pr(>F) as.factor(Fldf) 2 29.0984 14.5492 46.599 7.899e-15 *** Residuals 95 29.6608 0.3122 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > summary(lm(log(Zn) ~ as.factor(Fldf), data=utm.maas)) Call: lm(formula = log(Zn) ~ as.factor(Fldf), data = utm.maas) Residuals: Min 1Q Median 3Q Max -1.70515 -0.34744 -0.07094 0.29820 1.30316 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.48428 0.08521 76.097 < 2e-16 *** as.factor(Fldf)2 -1.06361 0.11853 -8.974 2.58e-14 *** as.factor(Fldf)3 -1.24547 0.20482 -6.081 2.49e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.5588 on 95 degrees of freedom Multiple R-Squared: 0.4952, Adjusted R-squared: 0.4846 F-statistic: 46.6 on 2 and 95 DF, p-value: 7.899e-15 > > > > ### *