################################################### if (nchar(Sys.getenv("GISRC")) == 0) stop("Start R inside GRASS") library(spgrass6) G <- gmeta6() .LOC <- G$LOCATION_NAME if (length(grep("^spearfish", .LOC)) == 0) stop("not in spearfish") ################################################### soilsph <- readRAST6("soils.ph", ignore.stderr=TRUE) ################################################### summary(soilsph) image(soilsph, "soils.ph", col=rev(cm.colors(6))) legend("bottom", legend=0:5, fill=rev(cm.colors(6)), cex=0.8, bty="n", horiz=TRUE) ################################################### spear <- readRAST6(c("elevation.dem", "landcover.30m"), cat=c(FALSE, TRUE), ignore.stderr=TRUE) ################################################### summary(spear) ################################################### grd <- slot(spear, "grid") lcovareas <- table(spear$landcover.30m) * (prod(slot(grd, "cellsize"))/10000) oopar <- par(mar=c(3,16,2,2)+0.1) boxplot(elevation.dem ~ landcover.30m, data=spear, medlwd=1, horizontal=TRUE, las=1, width=lcovareas) par(oopar) ################################################### bugsDF <- readVECT6("bugsites", ignore.stderr=TRUE) vInfo("streams", ignore.stderr=TRUE) streams <- readVECT6("streams", type="line,boundary", remove.duplicates=FALSE, ignore.stderr=TRUE) ################################################### image(spear, "elevation.dem", col=terrain.colors(20), axes=TRUE) plot(streams, col="blue", add=TRUE) plot(bugsDF, pch=19, col="grey25", cex=0.5, add=TRUE) ################################################### ## system("r.in.gdal input=BMcD_fldf.txt output=BMcD_fldf location=BMcD") ################################################### system("g.gisenv set='LOCATION_NAME=BMcD'") library(gstat) ################################################### ## system("v.in.ogr dsn=. layer=BMcD output=BMcD") ################################################### BMcD <- readVECT6("BMcD", ignore.stderr=TRUE) BMcD$Fldf <- factor(BMcD$Fldf) names(BMcD) ################################################### bubble(BMcD, "Zn") ################################################### boxplot(Zn ~ Fldf, BMcD, width=table(BMcD$Fldf), col="grey") ################################################### BMcD_grid <- as(readRAST6("BMcD_fldf", ignore.stderr=TRUE), "SpatialPixelsDataFrame") names(BMcD_grid) <- "Fldf" BMcD_grid$Fldf <- as.factor(BMcD_grid$Fldf) proj4string(BMcD) <- CRS(proj4string(BMcD_grid)) ################################################### pts <- list("sp.points", BMcD, pch = 4, col = "white") spplot(BMcD_grid, "Fldf", col.regions=1:3, sp.layout=list(pts)) ################################################### bluepal <- colorRampPalette(c("azure1", "steelblue4")) brks <- c(0,130,155,195,250,330,450,630,890,1270,1850) cols <- bluepal(length(brks)-1) sepal <- colorRampPalette(c("peachpuff1", "tomato3")) brks.se <- c(0,240,250,260,270,280,290,300,350,400,1000) cols.se <- sepal(length(brks.se)-1) scols <- c("green", "red") ################################################### cvgm <- variogram(Zn~1, data=BMcD, width=100, cutoff=1000) efitted <- fit.variogram(cvgm, vgm(psill=1, model="Exp", range=100, nugget=1)) efitted ################################################### plot(cvgm, model=efitted, plot.numbers=TRUE, col="black") ################################################### OK_fit <- gstat(id="OK_fit", formula = Zn ~ 1, data = BMcD, model=efitted) pe <- gstat.cv(OK_fit, debug.level=0, random=FALSE)$residual round(sqrt(mean(pe^2)), 2) z <- predict(OK_fit, newdata=BMcD_grid, debug.level=0) BMcD_grid$OK_pred <- z$OK_fit.pred BMcD_grid$OK_se <- sqrt(z$OK_fit.var) ################################################### image(BMcD_grid, "OK_pred", breaks=brks, col=cols) title("Fitted exponential OK model") symbols(coordinates(BMcD), circles=sqrt(abs(pe)), fg="black", bg=scols[(pe < 0)+1], inches=FALSE, add=TRUE) legend("topleft", fill=cols, legend=leglabs(brks), bty="n", cex=0.8) ################################################### image(BMcD_grid, "OK_se", breaks=brks.se, col=cols.se) title("Fitted exponential OK standard errors") symbols(coordinates(BMcD), circles=sqrt(abs(pe)), fg="black", bg=scols[(pe < 0)+1], inches=FALSE, add=TRUE) legend("topleft", fill=cols.se, legend=leglabs(brks.se), bty="n", cex=0.8) ################################################### cvgm <- variogram(Zn~Fldf, data=BMcD, width=100, cutoff=1000) uefitted <- fit.variogram(cvgm, vgm(psill=1, model="Exp", range=100, nugget=1)) uefitted ################################################### plot(cvgm, model=uefitted, plot.numbers=TRUE, col="black") ################################################### UK_fit <- gstat(id="UK_fit", formula = Zn ~ Fldf, data = BMcD, model=uefitted) pe_UK <- gstat.cv(UK_fit, debug.level=0, random=FALSE)$residual round(sqrt(mean(pe_UK^2)), 2) z <- predict(UK_fit, newdata=BMcD_grid, debug.level=0) BMcD_grid$UK_pred <- z$UK_fit.pred BMcD_grid$UK_se <- sqrt(z$UK_fit.var) ################################################### image(BMcD_grid, "UK_pred", breaks=brks, col=cols) title("Flood frequency UK model") symbols(coordinates(BMcD), circles=sqrt(abs(pe_UK)), fg="black", bg=scols[(pe_UK < 0)+1], inches=FALSE, add=TRUE) legend("topleft", fill=cols, legend=leglabs(brks), bty="n", cex=0.8) ################################################### image(BMcD_grid, "UK_se", breaks=brks.se, col=cols.se) title("Flood frequency UK interpolation standard errors") symbols(coordinates(BMcD), circles=sqrt(abs(pe_UK)), fg="black", bg=scols[(pe_UK < 0)+1], inches=FALSE, add=TRUE) legend("topleft", fill=cols.se, legend=leglabs(brks.se), bty="n", cex=0.8) ################################################### pts <- list("sp.points", BMcD, pch = 4, col = "black", cex=0.5) spplot(BMcD_grid, c("OK_pred", "UK_pred"), at=brks, col.regions=cols, sp.layout=list(pts)) ################################################### writeRAST6(BMcD_grid, vname="UK_pred", zcol="UK_pred", ignore.stderr=TRUE) cat(system("r.info UK_pred", intern=TRUE, ignore.stderr=TRUE)) ################################################### system("g.remove rast=UK_pred", ignore.stderr=TRUE) system(paste("g.gisenv set=\'LOCATION_NAME=",.LOC, "\'", sep=""))