[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