[R-sig-Geo] Help: R-function to plot all CSGS realization variograms in single a plot

[Edzer Pebesma]

Zia, time to find out about for-loops, and lists:

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

# Variogram log Zn

lzn.vgm = variogram(log(zinc) ~ 1, meuse)
lzn.fit = fit.variogram(lzn.vgm, model = vgm(1, "Sph", 900, 1))

#Conditional simulation

nsim = 100
lzn.sim = krige(log(zinc) ~ 1, meuse, meuse.grid, model = lzn.fit,
    nmax = 30, nsim = nsim)

# Variogram of all realizations

m = out = list()
for (i in 1:nsim) {
    s = paste("sim", i, sep="")
    f = as.formula(paste(s, "~1"))
    v = variogram(f, lzn.sim)
    v$id = rep(s, nrow(v))
    out[[s]] = v
    m[[s]] = fit.variogram(v, lzn.fit)
}
plot(do.call(rbind, out), m, layout=c(10,10), skip = FALSE,
    scales = list(y = list(relation = "same")))


Zia Ahmed wrote:
> Hi,
> I am trying to plot 100 variograms from 100 CSGS realizations  in a
> single plot with the observed  variogram model. I used following codes
> to plot them  together. But problem is this code is: I have to create
> variogram for all realizations,  ie  100  times I  have to repeat this 
> code!
> This their  any  way  to run following codes with in a loop to create a
> plot of variograms of all realizations with observed variogram model.
> Help will  be appreciated. Thanks
> Zia
> 
> 
> library(gstat)
> data(meuse)
> data(meuse.grid)
> coordinates(meuse) = ~x + y
> coordinates(meuse.grid) = ~x + y
> 
> # Variogram log Zn
> 
> lzn.vgm = variogram(log(zinc) ~ 1, meuse)
> lzn.fit = fit.variogram(lzn.vgm, model = vgm(1, "Sph", 900, 1))
> 
> #Conditional simulation
> 
> lzn.sim = krige(log(zinc) ~ 1, meuse, meuse.grid, model = lzn.fit,
> nmax = 30, nsim = 100)
> 
> # Variogram of all relizations
> 
> var.sim1<-variogram(sim1~1, lzn.sim)    # Realization 1
> var.sim2<-variogram(sim2~1, lzn.sim)    # Realization 2
> var.sim3<-variogram(sim3~1, lzn.sim)    # Realization 3
> var.sim4<-variogram(sim4~1, lzn.sim)    # Realization 4
> .
> .
> .
> var.sim100<-variogram(sim100~1, lzn.sim)    # Realization 100
> 
> # Vriogram of observed data and all relaizations
> ##################################################
> 
> plot(lzn.vgm$gamma ~ lzn.vgm$dist, xlim = c(0, max(lzn.vgm$dist) * 1.05),
> ylim = c(0, max(lzn.vgm$gamma) * 1.2), pch = 21, col = "black",cex=.8,
> cex.axis = 0.85, cex.lab=.8, xlab =  "Distance (m)", ylab = "Semivariance")
> 
> lines(variogramLine(fit.variogram(lzn.vgm,  lzn.fit),  # Sim1
> maxdist = max(lzn.vgm$dist)), col = "red", lwd=2)
> 
> lines(variogramLine(fit.variogram(var.sim1, lzn.fit),  # sim2
> maxdist = max(lzn.vgm$dist)), col = "black", lty=1)
> 
> lines(variogramLine(fit.variogram(var.sim2, lzn.fit), # Sim3
> maxdist = max(lzn.vgm$dist)), col = "black", lty=1)
> 
> lines(variogramLine(fit.variogram(var.sim3, lzn.fit), # Sim4
> maxdist = max(lzn.vgm$dist)), col = "black", lty=1)
> .
> .
> .
> .
> lines(variogramLine(fit.variogram(var.sim100, lzn.fit), # Sim100
> maxdist = max(lzn.vgm$dist)), col = "black", lty=1)
> 
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de