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]]

