[R-sig-Geo] Fwd: help - MC simulations at universal co-kriging

[Mercedes Román]

Dear r-sig-geo members,

I am trying to perform Monte-Carlo simulations for a universal co-kriging
system. The script I am using would work with a small dataset, as it is
with the meuse dataset. However, my calibration dataset has ~ 37000
observations and four variables (two providing the trend, two to be
modelled). The dataset where I am trying to do the simulation consists
aproximatelly of 2000 locations. Despite using the closest 20-30
observations for local co-kriging, and testing first with just 10 of the
2000 locations, and even for 2-5 simulations, R seems to not finish the
task.

Should I change anything else to speed up the calculations, or complement
gstat with another package?

Thanks in advance for your help,

Mercedes



data(meuse)
data(meuse.grid)

### check correlation with covariables
cor.test(log10(meuse$lead), log10(meuse$zinc))
cor.test(log10(meuse$lead), (meuse$dist))
cor.test(log10(meuse$zinc), (meuse$dist))

meuse$ltpb <- log10(meuse$lead)
meuse$ltzn <- log10(meuse$zinc)

### Transofrm to spatial
coordinates(meuse) <- ~ x + y
coordinates(meuse.grid) <- ~ x + y


# experimental variogam
v.ltpb <- variogram(ltpb ~ dist, data=meuse, cutoff=1800, width=200)
plot(v.ltpb, pl=T)
# estimate variogram model form and parameters by eye
m.ltpb <- vgm(0.035,"Sph",800,0.015)
plot(v.ltpb, pl=T, model=m.ltpb)
# fit model parameters
m.ltpb.f <- fit.variogram(v.ltpb, m.ltpb)
plot(v.ltpb, pl=T, model=m.ltpb.f)

### Set  local kriging in the gstat object
g <- gstat(NULL, id = "ltpb", form = ltpb ~ dist, data=meuse, nmax=30)
g <- gstat(g, id = "ltzn", form = ltzn ~ dist, data=meuse, nmax=30)

v.cross <- variogram(g)
str(v.cross)
plot(v.cross, pl=F)

g <- gstat(g, id = "ltpb", model = m.ltpb.f, fill.all=T, nmax=30)

g <- fit.lmc(v.cross, g, correct.diagonal = 1.01)
plot(v.cross, model=g$model)

# Universal co-kriging
k.c <- predict(g, meuse.grid)
str(k.c)
plot.kresults(k.c, "ltpb", meuse, "lead", meuse, "UcK with Zn covariable")

#### try MC simulation
MC <- 100
sims <- predict(g, meuse.grid, nsim = MC)

	[[alternative HTML version deleted]]