Skip to content

plotting 100 variograms on the same plot

5 messages · Paul Hiemstra, Laura Poggio

Original
Paul Hiemstra
#
Laura Poggio wrote:
Hi Laura,

The trick is to organize your data, an example using the meuse dataset:

library(lattice)

# Make a dataset with 100 fitted variogram models
# Don't know how you've organized your 100 models,
# But I put them in a list.
data(meuse)
coordinates(meuse) =~ x+y
av_list = lapply(1:100, function(x) autofitVariogram(zinc~1, meuse))

# Now extract the sample variograms, add an extra
# identification column to them and rbind them together
sv_list = do.call("rbind", lapply(1:length(av_list), function(x) {
sv = av_list[[x]]$exp_var
sv$identification = as.character(x)
return(sv)
}))

# Make a list of fitted models
model_list = lapply(av_list, function(x) x$var_model)

# Make the plot, notice the formula part!
# I use panel.number(), the current panel number,
# to extract the appropriate variogram model
xyplot(gamma ~ dist | identification, sv_list,
panel = function(...) {
panel.xyplot(...) # points of the sample var
ret = variogramLine(model_list[[panel.number()]], maxdist = 
max(sv_list$dist))
llines(ret$dist, ret$gamma, col = "red")
})

cheers,
Paul
Paul Hiemstra
#
Paul Hiemstra wrote:
...and ofcourse you need to do:

library(automap)

  
    

      
      
    

  


  


      

    

    

      
      

        
      

    

    
2 days later

    

      
      

        


  

        

      Paul Hiemstra
    

    

      
      #
    

  


  

      
Hi,

I think the only thing you need to change is the call to xyplot:

# So instead of "| identification" use groups = identification
xyplot(gamma ~ dist, groups = identification, sv_list,
panel = function(...) {
panel.xyplot(...) # points of the sample var
ret = variogramLine(model_list[[panel.number()]], maxdist = 
max(sv_list$dist))
llines(ret$dist, ret$gamma, col = "red")
})

This example above doesn't illustrate this point clearly, but it should 
work. You can choose to not show the sample variogram but only the 
models, comment the line "panel.xyplot(...)" in the call to xyplot to do 
this.

# So instead of "| identification" use groups = identification
xyplot(gamma ~ dist, groups = identification, sv_list,
panel = function(...) {
#panel.xyplot(...) # points of the sample var
ret = variogramLine(model_list[[panel.number()]], maxdist = 
max(sv_list$dist))
llines(ret$dist, ret$gamma, col = "red")
})

cheers,
Paul


Laura Poggio schreef: