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

[Zia Ahmed]

Thank you so much! It works.
However, I want to plot them in one figure with observed variogram model 
.  Is it possible? I know it will be looked very clumsy!
Thanks again
Zia

Edzer Pebesma wrote:
> 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
>>     
>
>   
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20100318/5bde31d0/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: zua3.vcf
Type: text/x-vcard
Size: 281 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20100318/5bde31d0/attachment.vcf>