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