[R-sig-Geo] plotting 100 variograms on the same plot

[Paul Hiemstra]

Laura Poggio wrote:
> Dear all,
> I would like to plot 100 variograms on the same plot, possibly using a loop.
>
> At the moment I am using autofitVariogram in the package Automap to fit the
> variograms automatically.
> I found some code calling the lattice package (
> http://tolstoy.newcastle.edu.au/R/e7/help/09/06/0940.html), but it is giving
> some problems for high numbers of variograms.
>
> Thank you very much in advance!
>
> Laura
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>   
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

-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +3130 274 3113 Mon-Tue
Phone:  +3130 253 5773 Wed-Fri
http://intamap.geo.uu.nl/~paul