[R-sig-Geo] Indicator Simulation Help

Katherine Lockhart lockhart.katherine at gmail.com
Wed Nov 7 21:24:13 CET 2012


Hi There,
I am trying to do conditional indicator simulation on some TCE values
in groundwater.
I set my data up as a SpatialPointsDataFrame and expanded a grid to
predict to that is a SpatialPixelsDataFrame.

>From there, here is my code:
for(i in 1:length(data$TCE_.mg.L.)) (if(data$TCE_.mg.L.[i]<=0.005)
data$TCE_.mg.L.ind[i] = 1 else
data$TCE_.mg.L.ind[i] = 0)
#Compute and plot the sample variogram for the indicator transform
TCE.var<-variogram(TCE_.mg.L.ind~1, data)
plot(TCE.var)
#Set initial choice of variogram model
# calculate variance as estimate of the sill
var(data$TCE_.mg.L.ind)
# [1] 0.2491694
TCE.vgm<-vgm(.25, "Sph", 40)
TCE.model<-fit.variogram(TCE.var, TCE.vgm)
TCE.model
#  model     psill    range
#1   Sph 0.2770698 50.82271
# try gstat and predict.gstat for simulation
g <- gstat(id = "TCE_.mg.L.ind", formula = TCE_.mg.L.ind~1,data =
data,model=TCE.model)
class(g)
#[1] "gstat" "list"
TCE.sim<-predict.gstat(g, site1.grid,nsim=1,indicators=TRUE)
summary(TCE.sim)
#Object of class SpatialPixelsDataFrame
#Coordinates:
# min max
#x  90 110
#y  75  85
#Is projected: NA
#proj4string : [NA]
#Number of points: 231
#Grid attributes:
# cellcentre.offset cellsize cells.dim
#x                90        1        21
#y                75        1        11
#Data attributes:
#  Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#0.00000 0.00000 0.00000 0.09524 0.00000 1.00000

My issue is how to get the simulated values for each pixel. I tried
TCE.sim[1] and I just get a list of the coordinates in X and Y. I also
need to extract the value of the mean. Also, the output seems off to
me, does this look typical of indicator simulation? If I plot with
spplot(TCE.sim[1]) I get a plot with just zeros and 1s. From my
understanding, indicator simulation gives values between 0-1. Any
advice is appreciated!
Thank you!
Katie

--
Katherine Lockhart
Graduate Student Researcher
UC Davis
Department of Land, Air, and Water Resources



More information about the R-sig-Geo mailing list