0

I have generated a 3D unconditional normal scores simulation using gstat package and would now like to check the variogram reproduction. How do a specify the calculation of the variogram on the grid?

#-- Create the grid
grd<-expand.grid(1:100,1:100,1:50)
names(grd)<-c("x","y","z")

#-- Define the variogram model
MyNug<-0.02
MySill<-0.98
MyRng<-70
MyAnis1<-1.0
MyAnis2<-0.1
vmdl<-vgm(nugget=MyNug,
          psill=MySill,
          range= MyRng, 
          model="Sph",
          anis=c(90,0,0,MyAnis1,MyAnis2))

#--- Prepare the simulation pars for an unconditional simulation
SimPar<-gstat(formula=p~1,
              locations=~x+y+z,
              dummy=T,
              beta=0,
              nmin=8,
              nmax=16,
              model=vmdl)

#--- Run  simulations
MyNsim=1
Sim<- predict(SimPar, newdata=grd, nsim=MyNsim,debug=-1)
summary(Sim)

#--- Prepare the variograms of the the simulation
#--- Convert Sim to a spatial object
SimS<-Sim
coordinates(SimS)<-~x+y+z
SimGrd<-grd
coordinates(SimGrd)<-~x+y+z
VarSim<-variogram(sim1~1,SimS,alpha=0, beta=0,grid=SimGrd)

At this point I get an error as follows

Error in variogram.default(y, locations, X, trend.beta = beta, grid = grid,  : 
  formal argument "grid" matched by multiple actual arguments

Any guidance or and example of variogram calculation using gridded data appreciated (note the simulation takes a while to run as there are 1 million nodes)

Markm0705
  • 1,340
  • 1
  • 13
  • 31

1 Answers1

2

When you do a

> gridded(SimS)=TRUE
> class(SimS)
[1] "SpatialPixelsDataFrame"
attr(,"package")
[1] "sp"
> VarSim<-variogram(sim1~1,SimS,alpha=0, beta=0)
Error: length of grid topology 9 unrecognized

this indicates that efficient variogram computation, using knowledge of the grid structure, is only implemented for 2D gridded data, not 3D.

> gridded(SimS)=FALSE
> VarSim<-variogram(sim1~1,SimS,alpha=0, beta=0)

works, but will be slow!

Edzer Pebesma
  • 3,814
  • 16
  • 26
  • Thank you for the prompt answer. It is a shame the code only deals efficiently with 2D data. Clearly it will be more efficient to do this planned investigation in GsLib using the Gam program – Markm0705 Sep 01 '15 at 13:25