Idea: the plot.vgm() function

[Nowosad]

gstat allows for plotting of variograms, and of variograms with models. However, it could be also useful to have an option to plot just a model alone, for example when using uncoditional simulation or for teaching purposes.

library(gstat)
library(sp)
data(meuse)
coordinates(meuse) = ~x+y

# variogram plot ----------------------------------------------------------
vario = variogram(log(zinc)~1, meuse)
plot(vario)

# variogram + model plot --------------------------------------------------
model = vgm(0.6, "Sph", 1000, nugget = 0.06)
plot(vario, model)        

# model plot? -------------------------------------------------------------
plot(model, sill = 0.7, cutoff = 1600)

[edzer]

As in

plot.variogramModel = function(x, ..., cutoff, type = 'l') {
    if (missing(cutoff))
        stop("parameter cutoff needs to be specified")
    plot(variogramLine(x, cutoff, ...), type = type)
}

?

[Nowosad]

Yes, this is exactly what I had in mind. Tiny suggestions - (a) change the line color to light blue, (b) change x and y labels to "distance" and "semiwariance", and (c) start the plot from 0 on the x and y axes to follow the convention of the previous plots.

library(gstat)
plot.variogramModel = function(x, ..., cutoff, type = 'l') {
  if (missing(cutoff))
    stop("parameter cutoff needs to be specified")
  plot(variogramLine(x, cutoff, ...), type = type)
}

# model -------------------------------------------------------------
model = vgm(0.6, "Sph", 1000, nugget = 0.06)
plot(model, cutoff = 1600)

^{Created on 2019-01-28 by the reprex package (v0.2.1)}

[edzer]

I'll look into; it should be similar to plot.gstatVariogram (using lattice, setting the x axis minimum to 0); you're right about the axis labels.

[Nowosad]

@edzer The function works great on a standard variogram.
Could it also work on directional variograms?

library(gstat)

# model -------------------------------------------------------------
model = vgm(0.6, "Sph", 1000, nugget = 0.06)
plot(model, cutoff = 1600)

# directional model -------------------------------------------------
model2 = vgm(0.6, "Sph", 1000, nugget = 0.06, anis = c(45, 0.4))
plot(model2, cutoff = 1600)

^{Created on 2019-01-28 by the reprex package (v0.2.1)}

[edzer]

Yes, if you specify the direction vector:

> plot(model2, cutoff = 1600, dir = c(.5 *sqrt(2),.5* sqrt(2),0))
> plot(model2, cutoff = 1600, dir = c(.5 *sqrt(2), -.5* sqrt(2),0))

see ?variogramLine

[Nowosad]

@edzer Would it be possible to add an alpha argument (like in, for example, variogram() function)?