R Users: ...been working with the sp and gstat packages for the past couple of days in an effort to analyze a set of ~ 200 soil samples collected from various eastings, northings, and depths and containing a wide range of measured hydrocarbon concentrations.
Thus far, I've managed to import the data, log-transform the concentrations, assign coordinates, generate and fit a variogram model and krige the data to a 3-D grid, but I've come up hard against the limits of my R proficiency. I'm looking for helpful suggestions on: - how to visualize the 3-D grid of expected values, - how to subset the data (such as 2-d slices of the grid in the xy or xz planes), - and how to do some relatively simple grid math operations such as removing all coordinates that have higher elevations than ground surface (I've made a 2-d grid of ground surface values), reverse the log-transformation, and ultimately estimate the total mass of hydrocarbons in the soil. Here's what I've gathered from the supporting documentation: - need to cast the 3-d expected concentration values into sp's SpatialPixelsDataFrame class in order to make any irregular subsets of the rectangular grid. - subsetting the grid probably involves use of "[ ...]" Here's my code so far: library(RODBC) library(gstat) options(digits = 10) channel = odbcDriverConnect("") sqlTables(channel) PAHdata <- sqlFetch(channel, "All") PAHdata attach(PAHdata) logTotalPAH = log10(TotalPAH) PAHdata = cbind(PAHdata, logTotalPAH) class(PAHdata) names(PAHdata) coordinates(PAHdata) = ~Easting + Northing + SampleElevation class(PAHdata) summary(PAHdata) dimensions(PAHdata) hist(logTotalPAH) lpah.vgm = variogram(logTotalPAH ~1, PAHdata) lpah.vgm plot(lpah.vgm ) lpah.fit = fit.variogram(lpah.vgm, model = vgm( 6, "Sph", 80, 1.5)) lpah.fit plot(lpah.vgm, lpah.fit) x = GridTopology(c(534500,3531400, 20), c(2,2, 0.5), c(176,126, 31)) class(x) coordinatevalues(x) x = SpatialGrid(x) class(x) grd = as(x, "SpatialPixels") lpah.kriged = krige(logTotalPAH~1, PAHdata, x, model = lpah.fit) class(lpah.kriged) lpah.kriged = as(lpah.kriged, "SpatialPixelsDataFrame") Any help would be appreciated. ...using R 2.7.0 on a windows XP desktop computer. ... will happily share scrubbed and anonomized data if requested. Matt Findley CH2M HILL Corvallis, Oregon [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.