[R-sig-Geo] Indicator Simulation, Simulated Values Help

[Katherine Lockhart]

Hi Everyone,
I am using Indicator Simulation to interpolate indicator transformed
TCE sample values from 40 monitoring wells to a
SpatialPixelsDataFrame. The grid that covers all of my sample points
in approximately 55,000 pixels (1m resolution) and I was not able to
get a single simulation to finish running. However, I am only
concerned with a 10X20m section within my larger grid. If I
interpolate only to this smaller grid (231 pixels), it runs instantly
and I can compute 1000 simulations in about 2 hours. The trouble I am
having is that each simulation gives highly variable results, either
all 0s or all 1s for example. The histogram of average field values
for each simulation is flat to bimodal. Also, if I simulate to another
smaller grid of interest, I get almost the same distribution of
average field values, though I have reason to believe that they should
not be the same. When I average the averages, I get nearly the name
number, 0.29.
Could this be cause by simulating only to a smaller grid? Or perhaps
my variogram model is incorrect?

TCE.vgm<-vgm(.25, "Sph", 40)
TCE.model<-fit.variogram(TCE.var, TCE.vgm)

g <- gstat(id = "TCE_.mg.L.ind", formula = TCE_.mg.L.ind~1,data =
data,model=TCE.model)

#Now perform the simulations for Site2
# set a new projection grid for the site squares to speed up predictions
#  set up a regular grid data.frame for my predictions to site 2
xs2<-seq(from=190,to=210,by=1)
ys2<-seq(from=55,to=55,by=1)
xysite2<-expand.grid(x=xs1,y=ys1)

# change the data frame class to a SpatialPixelsDataFrame
site2.grid<-SpatialPixelsDataFrame(points=xysite2[c('x','y')],
data=xysite2,proj4string=CRS(as.character(NA)))

# write the loop for site 2
vec2<-numeric(1000)
for (i in 1:1000) {
TCE.sim.i<-predict.gstat(g, site2.grid,nsim=i,indicators=TRUE)
average.i<-mean(TCE.sim.i[[i]])
vec2[i]<-average.i
}


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