[R-sig-Geo] gstat anisotropy : multipanel plotting of experiment and model

[solazef]

HI, thanks for the reply.

looking at the example(plot.gstatVariogram) , the type of plot produced by 
plot(vgm2, model=model.2)
(one panel for each direction, experiment and model overlapped)  is exactly
what  am trying to get.

Only in my case I don't have the direct output of the variogram function.
Basically I first run the variogram function to my dataset:

v <- variogram(dat0.z~1, data = dat, widht=width_d, alpha=hor_angle,
tol.hor=15, cloud=T)

with: hor_angle=seq(0,150,30)

Then I do some processing of the cloud to extract the distance relative to
one point only (not all pairs):

vm <- v[((v$np%/% attr(v, ".BigInt"))+1)== 1 | ((v$np%%  attr(v,
".BigInt"))+1)== 1, ]
    
 and split the vm according to the directions so to have six distinct
variograms to feed into the autofit function:
       
    vm_split <- split(vm,vm$dir.hor) 
    class (vm_split) <- c("gstatVariogram", "dataframe")
    vh =
list(vm_split[[1]],vm_split[[2]],vm_split[[3]],vm_split[[4]],vm_split[[5]],vm_split[[6]])
 
The output of the (modified) autofitVariogram fitting function is
"auto_Sph_h". For plotting I run:

for (i in 1:length(hor_angle)) { 
              model.1 <- vgm(psill=sum(auto_Sph_h[[i]]$var_model[,2]),
model="Sph", 
                              range=sum(auto_Sph_h[[i]]$var_model[,3]),
                             
nugget=sum(auto_Sph_h[[i]]$var_model[1,"psill"]), 
                              error=auto_Sph_h[[i]]$sserr)
        
          xyplot( gamma ~ dist | dir.hor
          ,data = vh[[i]] 
          ,labels = "", shift =0.1,  
          ,model = model.1
          ,direction = hor_angle 
          ,xlim=..., ylim=..., xlab=..., ylab=..., main=....
          ,group.id = TRUE           
           ,panel = automap:::autokrige.vgm.panel 
         ,mode = "direct")  }


 Now, this loop either does not show the model overlapped to the data or
print the plots in six different pages. I have tried countless combinations
of split.screen(), layout, par(mfrow), but nothing seems to work. Only by
using "print" six times with the option "more=T" does the job, but the
layout is not as nice.  

Hope you can suggest a way round.
Thanks.



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/gstat-anisotropy-multipanel-plotting-of-experiment-and-model-tp7584877p7584957.html
Sent from the R-sig-geo mailing list archive at Nabble.com.