[R-sig-Geo] gstat anisotropy : multipanel plotting of experiment and model
Edzer Pebesma
edzer.pebesma at uni-muenster.de
Thu Oct 24 16:07:27 CEST 2013
Dear Solazev, please study
http://www.r-project.org/posting-guide.html
carefully, and then try to formulate a question that will prompt a
useful answer.
On 10/24/2013 02:35 PM, solazef wrote:
> HI, thanks for the reply.
>
> looking at the example(plot.gstatVariogram) , the type of plot produced by
> plot(vgm2, model=model.2)
> (one panel for each direction, experiment and model overlapped) is exactly
> what am trying to get.
>
> Only in my case I don't have the direct output of the variogram function.
> Basically I first run the variogram function to my dataset:
>
> v <- variogram(dat0.z~1, data = dat, widht=width_d, alpha=hor_angle,
> tol.hor=15, cloud=T)
>
> with: hor_angle=seq(0,150,30)
>
> Then I do some processing of the cloud to extract the distance relative to
> one point only (not all pairs):
>
> vm <- v[((v$np%/% attr(v, ".BigInt"))+1)== 1 | ((v$np%% attr(v,
> ".BigInt"))+1)== 1, ]
>
> and split the vm according to the directions so to have six distinct
> variograms to feed into the autofit function:
>
> vm_split <- split(vm,vm$dir.hor)
> class (vm_split) <- c("gstatVariogram", "dataframe")
> vh =
> list(vm_split[[1]],vm_split[[2]],vm_split[[3]],vm_split[[4]],vm_split[[5]],vm_split[[6]])
>
> The output of the (modified) autofitVariogram fitting function is
> "auto_Sph_h". For plotting I run:
>
> for (i in 1:length(hor_angle)) {
> model.1 <- vgm(psill=sum(auto_Sph_h[[i]]$var_model[,2]),
> model="Sph",
> range=sum(auto_Sph_h[[i]]$var_model[,3]),
>
> nugget=sum(auto_Sph_h[[i]]$var_model[1,"psill"]),
> error=auto_Sph_h[[i]]$sserr)
>
> xyplot( gamma ~ dist | dir.hor
> ,data = vh[[i]]
> ,labels = "", shift =0.1,
> ,model = model.1
> ,direction = hor_angle
> ,xlim=..., ylim=..., xlab=..., ylab=..., main=....
> ,group.id = TRUE
> ,panel = automap:::autokrige.vgm.panel
> ,mode = "direct") }
>
>
> Now, this loop either does not show the model overlapped to the data or
> print the plots in six different pages. I have tried countless combinations
> of split.screen(), layout, par(mfrow), but nothing seems to work. Only by
> using "print" six times with the option "more=T" does the job, but the
> layout is not as nice.
>
> Hope you can suggest a way round.
> Thanks.
>
>
>
> --
> View this message in context: http://r-sig-geo.2731867.n2.nabble.com/gstat-anisotropy-multipanel-plotting-of-experiment-and-model-tp7584877p7584957.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251
83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 555 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20131024/415dc2c0/attachment.bin>
More information about the R-sig-Geo
mailing list