R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("geoR-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('geoR')

-------------------------------------------------------------
Functions for geostatistical data analysis
For an Introduction to geoR go to http://www.est.ufpr.br/geoR
geoR version 1.5-7 (built on 2005/06/07) is now loaded
-------------------------------------------------------------

> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "BoxCox"
> 
> ### * BoxCox
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: boxcox
> ### Title: The Box-Cox Transformed Normal Distribution
> ### Aliases: rboxcox dboxcox
> ### Keywords: distribution
> 
> ### ** Examples
> 
> ## Simulating data
> simul <- rboxcox(100, lambda=0.5, mean=10, sd=2)
> ##
> ## Comparing models with different lambdas,
> ## zero  means and unit variances
> curve(dboxcox(x, lambda=-1), 0, 8)
> for(lambda in seq(-.5, 1.5, by=0.5))
+   curve(dboxcox(x, lambda), 0, 8, add = TRUE)
> 
> 
> 
> cleanEx(); ..nameEx <- "InvChisquare"
> 
> ### * InvChisquare
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: InvChisquare
> ### Title: The (Scaled) Inverse Chi-Squared Distribution
> ### Aliases: dinvchisq rinvchisq
> ### Keywords: distribution
> 
> ### ** Examples
> 
> set.seed(1234); rinvchisq(5, df=2)
[1] 46.6089469  0.6693449  0.6358670  4.2780663  0.5423825
> set.seed(1234); 1/rchisq(5, df=2)
[1] 46.6089469  0.6693449  0.6358670  4.2780663  0.5423825
> 
> set.seed(1234); rinvchisq(5, df=2, scale=5)
[1] 466.089469   6.693449   6.358670  42.780663   5.423825
> set.seed(1234); 5*2/rchisq(5, df=2)
[1] 466.089469   6.693449   6.358670  42.780663   5.423825
> 
> ## inverse Chi-squared is a particular case
> x <- 1:10
> all.equal(dinvchisq(x, df=2), dinvchisq(x, df=2, scale=1/2))
[1] TRUE
> 
> 
> 
> cleanEx(); ..nameEx <- "Ksat"
> 
> ### * Ksat
> 
> flush(stderr()); flush(stdout())
> 
> ### Encoding: latin1
> 
> ### Name: Ksat
> ### Title: Saturated Hydraulic Conductivity
> ### Aliases: Ksat
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(Ksat)
> summary(Ksat)
Number of data points: 32 

Coordinates summary
       x   y
min  2.5 0.8
max 21.0 9.1

Distance summary
      min       max 
 0.509902 18.552897 

Borders summary
       x  y
min  0.0  0
max 22.5 13

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02000 0.09225 0.30500 0.87210 1.14300 4.30000 
> plot(Ksat, border=borders)
> 
> 
> 
> cleanEx(); ..nameEx <- "SIC"
> 
> ### * SIC
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: SIC
> ### Title: Spatial Interpolation Comparison data
> ### Aliases: sic SIC sic.100 sic.367 sic.all sic.some sic.borders
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(SIC)
> points(sic.100, borders=sic.borders)
> 
> 
> 
> cleanEx(); ..nameEx <- "as.geodata"
> 
> ### * as.geodata
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: as.geodata
> ### Title: Converts an Object to the Class "geodata"
> ### Aliases: as.geodata is.geodata
> ### Keywords: spatial classes manip
> 
> ### ** Examples
> 
> ## Not run: 
> ##D ## converting the data-set "topo" from the package MASS (VR's bundle)
> ##D ## to the geodata format:
> ##D require(MASS)
> ##D data(topo)
> ##D topo
> ##D topogeo <- as.geodata(topo)
> ##D names(topogeo)
> ##D topogeo
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "boxcox.fit"
> 
> ### * boxcox.fit
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: boxcox.fit
> ### Title: Parameter Estimation for the Box-Cox Transformation
> ### Aliases: boxcox.fit print.boxcox.fit plot.boxcox.fit lines.boxcox.fit
> ###   .negloglik.boxcox
> ### Keywords: regression models hplot
> 
> ### ** Examples
> 
> ## Simulating data
> simul <- rboxcox(100, lambda=0.5, mean=10, sd=2)
> ## Finding the ML estimates
> ml <- boxcox.fit(simul)
> ml
Fitted parameters:
    lambda       beta    sigmasq 
 0.5784082 12.3261397  5.6113236 

Convergence code returned by optim: 0
> ## Ploting histogram and fitted model
> plot(ml)
> ##
> ## Comparing models with different lambdas,
> ## zero  means and unit variances
> curve(dboxcox(x, lambda=-1), 0, 8)
> for(lambda in seq(-.5, 1.5, by=0.5))
+   curve(dboxcox(x, lambda), 0, 8, add = TRUE)
> ##
> ## Another example, now estimating lambda2
> ##
> simul <- rboxcox(100, lambda=0.5, mean=10, sd=2)
> ml <- boxcox.fit(simul, lambda2 = TRUE)
> ml
Fitted parameters:
      lambda      lambda2         beta      sigmasq 
-0.605558203 16.697691907  1.499075765  0.000375491 

Convergence code returned by optim: 0
> plot(ml)
> ##
> ## An example with a regression model
> ##
> if(require(MASS)){
+   data(trees)
+   boxcox.fit(object = trees[,3], xmat = trees[,1:2])
+   }
Loading required package: MASS
Fitted parameters:
    lambda      beta0      beta1      beta2    sigmasq 
 0.3065760 -2.7915276  0.4144820  0.0401038  0.0466829 

Convergence code returned by optim: 0
> 
> 
> 
> cleanEx(); ..nameEx <- "boxcox.geodata"
> 
> ### * boxcox.geodata
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: boxcox.geodata
> ### Title: Box-Cox transformation for geodata objects
> ### Aliases: boxcox.geodata
> ### Keywords: regression models hplot
> 
> ### ** Examples
> 
> if(require(MASS)){
+ data(wolfcamp)
+ boxcox(wolfcamp)
+ 
+ data(ca20)
+ boxcox(ca20, trend = ~altitude)
+ }
Loading required package: MASS
> 
> 
> 
> cleanEx(); ..nameEx <- "camg"
> 
> ### * camg
> 
> flush(stderr()); flush(stdout())
> 
> ### Encoding: latin1
> 
> ### Name: camg
> ### Title: Calcium and magnesium content in soil samples
> ### Aliases: camg
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(camg)
> plot(camg[-(1:2),])
> mg20 <- as.geodata(camg, data.col=6)
> plot(mg20)
> points(mg20)
> 
> 
> 
> cleanEx(); ..nameEx <- "cite.geoR"
> 
> ### * cite.geoR
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cite.geoR
> ### Title: Citing Package geoR in Publications
> ### Aliases: cite.geoR
> ### Keywords: misc
> 
> ### ** Examples
> 
> cite.geoR()

To cite geoR in publications, use

  RIBEIRO Jr., P.J. & DIGGLE, P.J. (2001) geoR: A package for
  geostatistical analysis. R-NEWS, Vol 1, No 2, 15-18. ISSN 1609-3631.

Please cite geoR when using it for data analysis!

A BibTeX entry for LaTeX users is

  @Article{,
     title	   = {{geoR}: a package for geostatistical analysis},
     author        = {Ribeiro Jr., P.J. and Diggle, P.J.},
     journal       = {R-NEWS},
     year	   = {2001},
     volume	   = {1},
     number	   = {2},
     pages	   = {15--18},
     issn          = {1609-3631},
     url           = {http://cran.R-project.org/doc/Rnews}
   }

> 
> 
> 
> cleanEx(); ..nameEx <- "coords.aniso"
> 
> ### * coords.aniso
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: coords.aniso
> ### Title: Geometric Anisotropy Correction
> ### Aliases: coords.aniso
> ### Keywords: spatial
> 
> ### ** Examples
> 
> op <- par(no.readonly = TRUE)
> par(mfrow=c(3,2))
> par(mar=c(2.5,0,0,0))
> par(mgp=c(2,.5,0))
> par(pty="s")
> ## Defining a set of coordinates
> coords <- expand.grid(seq(-1, 1, l=3), seq(-1, 1, l=5))
> plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
> text(coords[,1], coords[,2], 1:nrow(coords))
> ## Transforming coordinates according to some anisotropy parameters
> coordsA <- coords.aniso(coords, aniso.pars=c(0, 2))
> plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
> text(coordsA[,1], coordsA[,2], 1:nrow(coords))
> ##
> coordsB <- coords.aniso(coords, aniso.pars=c(pi/2, 2))
> plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
> text(coordsB[,1], coordsB[,2], 1:nrow(coords))
> ##
> coordsC <- coords.aniso(coords, aniso.pars=c(pi/4, 2))
> plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
> text(coordsC[,1], coordsC[,2], 1:nrow(coords))
> ##
> coordsD <- coords.aniso(coords, aniso.pars=c(3*pi/4, 2))
> plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
> text(coordsD[,1], coordsD[,2], 1:nrow(coords))
> ##
> coordsE <- coords.aniso(coords, aniso.pars=c(0, 5))
> plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
> text(coordsE[,1], coordsE[,2], 1:nrow(coords))
> ##
> par(op)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "coords2coords"
> 
> ### * coords2coords
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: coords2coords
> ### Title: Operations on Coordinates
> ### Aliases: coords2coords zoom.coords zoom.coords.default
> ###   zoom.coords.geodata rect.coords
> ### Keywords: spatial
> 
> ### ** Examples
> 
> foo <- matrix(c(4,6,6,4,2,2,4,4), nc=2)
> foo1 <- zoom.coords(foo, 2)
> foo1
     [,1] [,2]
[1,]    3    1
[2,]    7    1
[3,]    7    5
[4,]    3    5
> foo2 <- coords2coords(foo, c(6,10), c(6,10))
> foo2
     [,1] [,2]
[1,]    6    6
[2,]   10    6
[3,]   10   10
[4,]    6   10
> plot(1:10, 1:10, type="n")
> polygon(foo)
> polygon(foo1, lty=2)
> polygon(foo2, lwd=2)
> arrows(foo[,1], foo[,2],foo1[,1],foo1[,2], lty=2)
> arrows(foo[,1], foo[,2],foo2[,1],foo2[,2])
> legend(1,10, c("foo", "foo1 (zoom.coords)", "foo2 (coords2coords)"), lty=c(1,2,1), lwd=c(1,1,2), cex=1.7)
> 
> ## "zooming" part of The Gambia map
> data(gambia)
> gb <- gambia.borders/1000
> gd <- gambia[,1:2]/1000
> plot(gb, ty="l", asp=1, xlab="W-E (kilometres)", ylab="N-S (kilometres)")
> points(gd, pch=19, cex=0.5)
> r1b <- gb[gb[,1] < 420,]
> rc1 <- rect.coords(r1b, lty=2)
> 
> r1bn <- zoom.coords(r1b, 1.8, xoff=90, yoff=-90)
> rc2 <- rect.coords(r1bn, xz=1.05)
> segments(rc1[c(1,3),1],rc1[c(1,3),2],rc2[c(1,3),1],rc2[c(1,3),2], lty=3)
> 
> lines(r1bn)
> r1d <- gd[gd[,1] < 420,]
> r1dn <- zoom.coords(r1d, 1.7, xlim.o=range(r1b[,1],na.rm=TRUE), ylim.o=range(r1b[,2],na.rm=TRUE), xoff=90, yoff=-90)
> points(r1dn, pch=19, cex=0.5)
> text(450,1340, "Western Region", cex=1.5)
> 
> 
> 
> cleanEx(); ..nameEx <- "cov.spatial"
> 
> ### * cov.spatial
> 
> flush(stderr()); flush(stdout())
> 
> ### Encoding: latin1
> 
> ### Name: cov.spatial
> ### Title: Computes Value of the Covariance Function
> ### Aliases: cov.spatial .cor.number .check.cov.model
> ### Keywords: spatial models
> 
> ### ** Examples
> 
> #
> # Variogram models with the same "practical" range:
> #
> v.f <- function(x, ...){1-cov.spatial(x, ...)}
> #
> curve(v.f(x, cov.pars=c(1, .2)), from = 0, to = 1,
+       xlab = "distance", ylab = expression(gamma(h)),
+       main = "variograms with equivalent \"practical range\"")
> curve(v.f(x, cov.pars = c(1, .6), cov.model = "sph"), 0, 1,
+       add = TRUE, lty = 2)
> curve(v.f(x, cov.pars = c(1, .6/sqrt(3)), cov.model = "gau"),
+       0, 1, add = TRUE, lwd = 2)
> legend(0.5,.3, c("exponential", "spherical", "gaussian"),
+        lty=c(1,2,1), lwd=c(1,1,2))
> #
> # Matern models with equivalent "practical range"
> # and varying smoothness parameter
> #
> curve(v.f(x, cov.pars = c(1, 0.25), kappa = 0.5),from = 0, to = 1,
+       xlab = "distance", ylab = expression(gamma(h)), lty = 2,
+       main = "models with equivalent \"practical\" range")
> curve(v.f(x, cov.pars = c(1, 0.188), kappa = 1),from = 0, to = 1,
+       add = TRUE)      
> curve(v.f(x, cov.pars = c(1, 0.14), kappa = 2),from = 0, to = 1,
+       add = TRUE, lwd=2, lty=2)      
> curve(v.f(x, cov.pars = c(1, 0.117), kappa = 2),from = 0, to = 1,
+       add = TRUE, lwd=2)      
> legend(0.4,.4, c(expression(paste(kappa == 0.5, "  and  ",
+        phi == 0.250)),
+        expression(paste(kappa == 1, "  and  ", phi == 0.188)),
+        expression(paste(kappa == 2, "  and  ", phi == 0.140)),
+        expression(paste(kappa == 3, "  and  ", phi == 0.117))),
+        lty=c(2,1,2,1), lwd=c(1,1,2,2))
> # plotting a nested variogram model
> curve(v.f(x, cov.pars = rbind(c(.4, .2), c(.6,.3)),
+           cov.model = c("sph","exp")), 0, 1, ylab='nested model')
Warning in if (cov.model == "stable") cov.model <- "powered.exponential" : 
	 the condition has length > 1 and only the first element will be used
> 
> 
> 
> cleanEx(); ..nameEx <- "dup.coords"
> 
> ### * dup.coords
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: dup.coords
> ### Title: Locates duplicated coordinates
> ### Aliases: dup.coords dup.coords.default dup.coords.geodata
> ### Keywords: spatial manip
> 
> ### ** Examples
> 
> ## simulating data
> foo <- grf(30, cov.p=c(1, .3)) 
grf: simulation(s) on randomly chosen locations with  30  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0 
grf: covariance model 1 is: exponential(sigmasq=1, phi=0.3)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 1 
> ## forcing some duplicated locations
> foo$coords[4,] <- foo$coords[14,] <- foo$coords[24,] <- foo$coords[2,]
> foo$coords[17,] <- foo$coords[23,] <- foo$coords[8,]
> ## locations the duplicated coordinates
> dup.coords(foo)
[[1]]
[1]  2  4 14 24

[[2]]
[1]  8 17 23

> 
> 
> 
> cleanEx(); ..nameEx <- "elevation"
> 
> ### * elevation
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: elevation
> ### Title: Surface Elevations
> ### Aliases: elevation
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(elevation)
> summary(elevation)
Number of data points: 52 

Coordinates summary
      x   y
min 0.2 0.0
max 6.3 6.2

Distance summary
     min      max 
0.200000 8.275869 

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  690.0   787.5   830.0   827.1   873.0   960.0 
> plot(elevation)
> 
> 
> 
> cleanEx(); ..nameEx <- "gambia"
> 
> ### * gambia
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: gambia
> ### Title: Gambia Malaria Data
> ### Aliases: gambia gambia.borders gambia.map
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(gambia)
> plot(gambia.borders, type="l", asp=1)
> points(gambia[,1:2], pch=19)
> # a built-in function for a zoomed map
> gambia.map()
> # Building a "village-level" data frame
> ind <- paste("x",gambia[,1], "y", gambia[,2], sep="")
> village <- gambia[!duplicated(ind),c(1:2,7:8)]
> village$prev <- as.vector(tapply(gambia$pos, ind, mean))
> plot(village$green, village$prev)
> 
> 
> 
> cleanEx(); ..nameEx <- "grf"
> 
> ### * grf
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: grf
> ### Title: Simulation of Gaussian Random Fields
> ### Aliases: grf geoR2RF .grf.aux1 grfclass lines.grf
> ### Keywords: spatial datagen
> 
> ### ** Examples
> 
> sim1 <- grf(100, cov.pars = c(1, .25))
grf: simulation(s) on randomly chosen locations with  100  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0 
grf: covariance model 1 is: exponential(sigmasq=1, phi=0.25)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 1 
> # a display of simulated locations and values
> points(sim1)   
> # empirical and theoretical variograms
> plot(sim1)
variog: computing omnidirectional variogram
> #
> # a "smallish" simulation
> sim2 <- grf(441, grid = "reg", cov.pars = c(1, .25)) 
grf: generating grid  21  *  21  with  441  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0 
grf: covariance model 1 is: exponential(sigmasq=1, phi=0.25)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 1 
> image(sim2)
> ##
> ## 1-D simulations using the same seed and different noise/signal ratios
> ##
> set.seed(234)
> sim11 <- grf(100, ny=1, cov.pars=c(1, 0.25), nug=0)
simulations in 1D
grf: simulation(s) on randomly chosen locations with  100  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0 
grf: covariance model 1 is: exponential(sigmasq=1, phi=0.25)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 1 
> set.seed(234)
> sim12 <- grf(100, ny=1, cov.pars=c(0.75, 0.25), nug=0.25)
simulations in 1D
grf: simulation(s) on randomly chosen locations with  100  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0.25 
grf: covariance model 1 is: exponential(sigmasq=0.75, phi=0.25)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 1 
> set.seed(234)
> sim13 <- grf(100, ny=1, cov.pars=c(0.5, 0.25), nug=0.5)
simulations in 1D
grf: simulation(s) on randomly chosen locations with  100  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0.5 
grf: covariance model 1 is: exponential(sigmasq=0.5, phi=0.25)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 1 
> ##
> par.ori <- par(no.readonly = TRUE)
> par(mfrow=c(3,1), mar=c(3,3,.5,.5))
> yl <- range(c(sim11$data, sim12$data, sim13$data))
> image(sim11, type="l", ylim=yl)
> image(sim12, type="l", ylim=yl)
> image(sim13, type="l", ylim=yl)
> par(par.ori)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "head"
> 
> ### * head
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: head
> ### Title: Head observations in a regional confined aquifer
> ### Aliases: head
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(head)
> summary(head)
Number of data points: 29 

Coordinates summary
        x     y
min  4.32  3.38
max 16.27 11.41

Distance summary
      min       max 
 0.254951 12.197713 

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  584.0   676.0   782.0   830.1   998.0  1202.0 
> plot(head)
> 
> 
> 
> cleanEx(); ..nameEx <- "hist.krige.bayes"
> 
> ### * hist.krige.bayes
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: hist.krige.bayes
> ### Title: Plots Sample from Posterior Distributions
> ### Aliases: hist.krige.bayes
> ### Keywords: spatial dplot
> 
> ### ** Examples
> 
> ## See documentation for krige.bayes()
> 
> 
> 
> cleanEx(); ..nameEx <- "hoef"
> 
> ### * hoef
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: hoef
> ### Title: Data for spatial analysis of experiments
> ### Aliases: hoef
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(hoef)
> hoef.geo <- as.geodata(hoef, covar.col=4)
> summary(hoef)
       x1          x2         dat             trat         ut       
 Min.   :1   Min.   :1   Min.   :16.00   Min.   :1   Min.   :20.00  
 1st Qu.:2   1st Qu.:2   1st Qu.:21.00   1st Qu.:2   1st Qu.:23.00  
 Median :3   Median :3   Median :24.00   Median :3   Median :24.00  
 Mean   :3   Mean   :3   Mean   :24.88   Mean   :3   Mean   :24.08  
 3rd Qu.:4   3rd Qu.:4   3rd Qu.:29.00   3rd Qu.:4   3rd Qu.:25.00  
 Max.   :5   Max.   :5   Max.   :32.00   Max.   :5   Max.   :32.00  
> summary(hoef.geo)
Number of data points: 25 

Coordinates summary
    x1 x2
min  1  1
max  5  5

Distance summary
     min      max 
1.000000 5.656854 

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  16.00   21.00   24.00   24.88   29.00   32.00 

Covariates summary
      trat  
 Min.   :1  
 1st Qu.:2  
 Median :3  
 Mean   :3  
 3rd Qu.:4  
 Max.   :5  
> points(hoef.geo, cex.min=2, cex.max=2, pt.div="quintiles")
> 
> 
> 
> cleanEx(); ..nameEx <- "image.grf"
> 
> ### * image.grf
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: image.grf
> ### Title: Image, Contour or Perspective Plot of Simulated Gaussian Random
> ###   Field
> ### Aliases: image.grf contour.grf persp.grf
> ### Keywords: spatial dplot
> 
> ### ** Examples
> 
> # generating 4 simulations of a Gaussian random field
> sim <- grf(441, grid="reg", cov.pars=c(1, .25), nsim=4)
grf: generating grid  21  *  21  with  441  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0 
grf: covariance model 1 is: exponential(sigmasq=1, phi=0.25)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 4 
> op <- par(no.readonly = TRUE)
> par(mfrow=c(2,2), mar=c(3,3,1,1), mgp = c(2,1,0))
> for (i in 1:4)
+   image(sim, sim.n=i)
> par(op)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "image.krige.bayes"
> 
> ### * image.krige.bayes
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: image.krige.bayes
> ### Title: Plots Results of the Predictive Distribution
> ### Aliases: image.krige.bayes contour.krige.bayes persp.krige.bayes
> ###   .prepare.graph.krige.bayes
> ### Keywords: spatial
> 
> ### ** Examples
> 
> #See examples in the documentation for the function krige.bayes().
> 
> 
> 
> cleanEx(); ..nameEx <- "image.kriging"
> 
> ### * image.kriging
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: image.kriging
> ### Title: Image or Perspective Plot with Kriging Results
> ### Aliases: image.kriging contour.kriging persp.kriging
> ###   .prepare.graph.kriging plot.1d
> ### Keywords: spatial dplot
> 
> ### ** Examples
> 
> data(s100) 
> loci <- expand.grid(seq(0,1,l=31), seq(0,1,l=31))
> kc <- krige.conv(s100, loc=loci,
+                  krige=krige.control(cov.pars=c(1, .25)))
krige.conv: model with constant mean
krige.conv: Kriging performed using global neighbourhood 
> image(kc)
> contour(kc)
> image(kc)
> contour(kc, add=TRUE, nlev=21)
> persp(kc, theta=20, phi=20)
> contour(kc, filled=TRUE)
> contour(kc, filled=TRUE, color=terrain.colors)
> contour(kc, filled=TRUE, col=gray(seq(1,0,l=21)))
> # adding data locations
> image(kc, coords.data=s100$coords)
> contour(kc,filled=TRUE,coords.data=s100$coords,color=terrain.colors)
> #
> # now dealing with borders
> #
> bor <- matrix(c(.4,.1,.3,.9,.9,.7,.9,.7,.3,.2,.5,.8),
+               ncol=2)
> # plotting just inside borders
> image(kc, borders=bor)
Loading required package: splancs

Spatial Point Pattern Analysis Code in S-Plus

 Version 2 - Spatial and Space-Time analysis
> contour(kc, borders=bor)
> image(kc, borders=bor)
> contour(kc, borders=bor, add=TRUE)
> contour(kc, borders=bor, filled=TRUE, color=terrain.colors)
> # kriging just inside borders
> kc1 <- krige.conv(s100, loc=loci,
+                  krige=krige.control(cov.pars=c(1, .25)),
+                  borders=bor)
krige.conv: results will be returned only for prediction locations inside the borders
krige.conv: model with constant mean
krige.conv: Kriging performed using global neighbourhood 
> image(kc1)
> contour(kc1)
> # avoiding the borders 
> image(kc1, borders=NULL)
> contour(kc1, borders=NULL)
> 
> op <- par(no.readonly = TRUE)
> par(mfrow=c(1,2), mar=c(3,3,0,0), mgp=c(1.5, .8,0))
> image(kc)
> image(kc, val=sqrt(kc$krige.var))
> 
> # different ways to add the legends and pass arguments:
> image(kc, ylim=c(-0.2, 1), x.leg=c(0,1), y.leg=c(-0.2, -0.1))
> image(kc, val=kc$krige.var, ylim=c(-0.2, 1))
> legend.krige(y.leg=c(-0.2,-0.1), x.leg=c(0,1), val=sqrt(kc$krige.var))
> 
> image(kc, ylim=c(-0.2, 1), x.leg=c(0,1), y.leg=c(-0.2, -0.1), cex=1.5)
> image(kc, ylim=c(-0.2, 1), x.leg=c(0,1), y.leg=c(-0.2, -0.1), offset.leg=0.5)
> 
> image(kc, xlim=c(0, 1.2))
> legend.krige(x.leg=c(1.05,1.1), y.leg=c(0,1),kc$pred, vert=TRUE)
> image(kc, xlim=c(0, 1.2))
> legend.krige(x.leg=c(1.05,1.1), y.leg=c(0,1),kc$pred, vert=TRUE, off=1.5, cex=1.5)
> 
> par(op)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "isaaks"
> 
> ### * isaaks
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: isaaks
> ### Title: Data from Isaaks and Srisvastava's book
> ### Aliases: isaaks
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(isaaks)
> isaaks
$coords
     [,1] [,2]
[1,]   61  139
[2,]   63  140
[3,]   64  129
[4,]   68  128
[5,]   71  140
[6,]   73  141
[7,]   75  128

$data
[1] 477 696 227 646 606 791 783

$x0
     [,1] [,2]
[1,]   65  137

attr(,"class")
[1] "geodata"
> summary(isaaks)
Number of data points: 7 

Coordinates summary
    Coord.X Coord.Y
min      61     128
max      75     141

Distance summary
      min       max 
 2.236068 17.804494 

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  227.0   541.5   646.0   603.7   739.5   791.0 

Other elements in the geodata object
[1] "x0"
> plot(isaaks$coords, asp=1, type="n")
> text(isaaks$coords, as.character(isaaks$data))
> points(isaaks$x0, pch="?", cex=2, col=2)
> 
> 
> 
> cleanEx(); ..nameEx <- "krige.bayes"
> 
> ### * krige.bayes
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: krige.bayes
> ### Title: Bayesian Analysis for Gaussian Geostatistical Models
> ### Aliases: krige.bayes model.control prior.control post2prior
> ###   print.krige.bayes print.posterior.krige.bayes
> ### Keywords: spatial models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D # generating a simulated data-set
> ##D ex.data <- grf(50, cov.pars=c(10, .25))
> ##D #
> ##D # defining the grid of prediction locations:
> ##D ex.grid <- as.matrix(expand.grid(seq(0,1,l=21), seq(0,1,l=21)))
> ##D #
> ##D # computing posterior and predictive distributions
> ##D # (warning: the next command can be time demanding)
> ##D ex.bayes <- krige.bayes(ex.data, loc=ex.grid, prior =
> ##D                  prior.control(phi.discrete=seq(0, 2, l=21)))
> ##D #
> ##D # Prior and posterior for the parameter phi
> ##D plot(ex.bayes, type="h", tausq.rel = FALSE, col=c("red", "blue"))
> ##D #
> ##D # Plot histograms with samples from the posterior
> ##D par(mfrow=c(3,1))
> ##D hist(ex.bayes)
> ##D par(mfrow=c(1,1))
> ##D 
> ##D # Plotting empirical variograms and some Bayesian estimates:
> ##D # Empirical variogram
> ##D plot(variog(ex.data, max.dist = 1), ylim=c(0, 15))
> ##D # Since ex.data is a simulated data we can plot the line with the "true" model 
> ##D lines(ex.data)
> ##D # adding lines with summaries of the posterior of the binned variogram
> ##D lines(ex.bayes, summ = mean, lwd=1, lty=2)
> ##D lines(ex.bayes, summ = median, lwd=2, lty=2)
> ##D # adding line with summary of the posterior of the parameters
> ##D lines(ex.bayes, summary = "mode", post = "parameters")
> ##D 
> ##D # Plotting again the empirical variogram
> ##D plot(variog(ex.data, max.dist=1), ylim=c(0, 15))
> ##D # and adding lines with median and quantiles estimates
> ##D my.summary <- function(x){quantile(x, prob = c(0.05, 0.5, 0.95))}
> ##D lines(ex.bayes, summ = my.summary, ty="l", lty=c(2,1,2), col=1)
> ##D 
> ##D # Plotting some prediction results
> ##D op <- par(no.readonly = TRUE)
> ##D par(mfrow=c(2,2))
> ##D par(mar=c(3,3,1,1))
> ##D par(mgp = c(2,1,0))
> ##D image(ex.bayes, main="predicted values")
> ##D image(ex.bayes, val="variance", main="prediction variance")
> ##D image(ex.bayes, val= "simulation", number.col=1,
> ##D       main="a simulation from the \npredictive distribution")
> ##D image(ex.bayes, val= "simulation", number.col=2,
> ##D       main="another simulation from \nthe predictive distribution")
> ##D #
> ##D par(op)
> ## End(Not run)
> ##
> ## For a extended list of exemples of the usage of krige.bayes()
> ## see http://www.est.ufpr.br/geoR/geoRdoc/examples.krige.bayes.R
> ##
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "krige.conv"
> 
> ### * krige.conv
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: krige.conv
> ### Title: Spatial Prediction - Conventional Kriging
> ### Aliases: krige.conv krige.control
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Not run: 
> ##D data(s100)
> ##D # Defining a prediction grid
> ##D loci <- expand.grid(seq(0,1,l=21), seq(0,1,l=21))
> ##D # predicting by ordinary kriging
> ##D kc <- krige.conv(s100, loc=loci,
> ##D                  krige=krige.control(cov.pars=c(1, .25)))
> ##D # mapping point estimates and variances
> ##D par.ori <- par(no.readonly = TRUE)
> ##D par(mfrow=c(1,2), mar=c(3.5,3.5,1,0), mgp=c(1.5,.5,0))
> ##D image(kc, main="kriging estimates")
> ##D image(kc, val=sqrt(kc$krige.var), main="kriging std. errors")
> ##D # Now setting the output to simulate from the predictive
> ##D # (obtaining conditional simulations),
> ##D # and to compute quantile and probability estimators
> ##D s.out <- output.control(n.predictive = 1000, quant=0.9, thres=2)
> ##D set.seed(123)
> ##D kc <- krige.conv(s100, loc = loci,
> ##D          krige = krige.control(cov.pars = c(1,.25)),
> ##D          output = s.out)
> ##D par(mfrow=c(2,2))
> ##D image(kc, val=kc$sim[,1], main="a cond. simul.")
> ##D image(kc, val=kc$sim[,1], main="another cond. simul.")
> ##D image(kc, val=(1 - kc$prob), main="Map of P(Y > 2)")
> ##D image(kc, val=kc$quant, main="Map of y s.t. P(Y < y) = 0.9")
> ##D par(par.ori)
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "krweights"
> 
> ### * krweights
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: krweights
> ### Title: Computes kriging weights
> ### Aliases: krweights
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Figure 8.4 in Webster and Oliver (2001), see help(wo)
> data(wo)
> attach(wo)
> par(mfrow=c(2,2))
> plot(c(-10,130), c(-10,130), ty="n", asp=1)
> points(rbind(coords, x1))
> KC1 <- krige.control(cov.pars=c(0.382,90.53))
> w1 <- krweights(wo$coords, loc=x1, krige=KC1)
> text(coords[,1], 5+coords[,2], round(w1, dig=3))
> ##
> plot(c(-10,130), c(-10,130), ty="n", asp=1)
> points(rbind(coords, x1))
> KC2 <- krige.control(cov.pars=c(0.282,90.53), nug=0.1)
> w2 <- krweights(wo$coords, loc=x1, krige=KC2)
> text(coords[,1], 5+coords[,2], round(w2, dig=3))
> ##
> plot(c(-10,130), c(-10,130), ty="n", asp=1)
> points(rbind(coords, x1))
> KC3 <- krige.control(cov.pars=c(0.082,90.53), nug=0.3)
> w3 <- krweights(wo$coords, loc=x1, krige=KC3)
> text(coords[,1], 5+coords[,2], round(w3, dig=3))
> ##
> plot(c(-10,130), c(-10,130), ty="n", asp=1)
> points(rbind(coords, x1))
> KC4 <- krige.control(cov.pars=c(0,90.53), nug=0.382, micro=0.382)
> w4 <- krweights(wo$coords, loc=x1, krige=KC4)
> text(coords[,1], 5+coords[,2], round(w4, dig=3))
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "ksline"
> 
> ### * ksline
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ksline
> ### Title: Spatial Prediction - Conventional Kriging
> ### Aliases: ksline .ksline.aux.1
> ### Keywords: spatial
> 
> ### ** Examples
> 
> data(s100)
> loci <- expand.grid(seq(0,1,l=31), seq(0,1,l=31))
> kc <- ksline(s100, loc=loci, cov.pars=c(1, .25))
ksline: kriging location:  1 out of 961 
ksline: kriging location:  101 out of 961 
ksline: kriging location:  201 out of 961 
ksline: kriging location:  301 out of 961 
ksline: kriging location:  401 out of 961 
ksline: kriging location:  501 out of 961 
ksline: kriging location:  601 out of 961 
ksline: kriging location:  701 out of 961 
ksline: kriging location:  801 out of 961 
ksline: kriging location:  901 out of 961 
ksline: kriging location:  961 out of 961 
Kriging performed using global neighbourhood 
> par(mfrow=c(1,2))
> image(kc, main="kriging estimates")
> image(kc, val=sqrt(kc$krige.var), main="kriging std. errors")
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "legend.krige"
> 
> ### * legend.krige
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: legend.krige
> ### Title: Add a legend to a image with kriging results
> ### Aliases: legend.krige
> ### Keywords: spatial aplot
> 
> ### ** Examples
> 
> # See examples in the documentation for image.kriging
> 
> 
> 
> cleanEx(); ..nameEx <- "likfit"
> 
> ### * likfit
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: likfit
> ### Title: Likelihood Based Parameter Estimation for Gaussian Random Fields
> ### Aliases: likfit likfit.limits .negloglik.GRF
> ### Keywords: spatial models
> 
> ### ** Examples
> 
> ## Not run: 
> ##D data(s100)
> ##D ml <- likfit(s100, ini=c(0.5, 0.5), fix.nug = TRUE)
> ##D ml
> ##D summary(ml)
> ##D reml <- likfit(s100, ini=c(0.5, 0.5), fix.nug = TRUE, met = "REML")
> ##D summary(reml)
> ##D plot(variog(s100))
> ##D lines(ml)
> ##D lines(reml, lty = 2)
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "lines.krige.bayes"
> 
> ### * lines.krige.bayes
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lines.variomodel.krige.bayes
> ### Title: Adds a Bayesian Estimate of the Variogram to a Plot
> ### Aliases: lines.variomodel.krige.bayes
> ### Keywords: spatial aplot
> 
> ### ** Examples
> 
> #See examples in the documentation of the function krige.bayes().
> 
> 
> 
> cleanEx(); ..nameEx <- "lines.variogram.envelope"
> 
> ### * lines.variogram.envelope
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lines.variogram.envelope
> ### Title: Adds Envelopes Lines to a Variogram Plot
> ### Aliases: lines.variogram.envelope
> ### Keywords: spatial aplot
> 
> ### ** Examples
> 
> if(is.R()) data(s100)
> s100.vario <- variog(s100, max.dist = 1)
variog: computing omnidirectional variogram
> s100.ml <- likfit(s100, ini=c(.5, .5))
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optim.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.
> s100.mod.env <- variog.model.env(s100, obj.variog = s100.vario,
+    model = s100.ml) 
variog.env: generating 99 simulations (with  100 points each) using the function grf
variog.env: adding the mean or trend
variog.env: computing the empirical variogram for the 99 simulations
variog.env: computing the envelops
> s100.mc.env <- variog.mc.env(s100, obj.variog = s100.vario)
variog.env: generating 99 simulations by permutating data values
variog.env: computing the empirical variogram for the 99 simulations
variog.env: computing the envelops
> plot(s100.vario)
> lines(s100.mod.env)
> lines(s100.mc.env, lwd=2)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "lines.variomodel"
> 
> ### * lines.variomodel
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lines.variomodel
> ### Title: Adds a Line with a Variogram Model to a Variogram Plot
> ### Aliases: lines.variomodel lines.variomodel.default
> ### Keywords: spatial aplot
> 
> ### ** Examples
> 
> data(s100)
> # computing and ploting empirical variogram
> vario <- variog(s100, max.dist = 1)
variog: computing omnidirectional variogram
> plot(vario)
> # estimating parameters by weighted least squares
> vario.wls <- variofit(vario, ini = c(1, .3), fix.nugget = TRUE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
> # adding fitted model to the plot  
> lines(vario.wls)
> #
> # Ploting different variogram models
> plot(0:1, 0:1, type="n")
> lines.variomodel(cov.model = "exp", cov.pars = c(.7, .25), nug = 0.3, max.dist = 1) 
> # an alternative way to do this is:
> my.model <- list(cov.model = "exp", cov.pars = c(.7, .25), nugget = 0.3,
+ max.dist = 1) 
> lines.variomodel(my.model, lwd = 2)
> # now adding another model
> lines.variomodel(cov.m = "mat", cov.p = c(.7, .25), nug = 0.3,
+                  max.dist = 1, kappa = 1, lty = 2)
> # adding the so-called "nested" models
> # two exponential structures
> lines.variomodel(seq(0,1,l=101), cov.model="exp",
+                  cov.pars=rbind(c(0.6,0.15),c(0.4,0.25)), nug=0, col=2)
Warning in if (cov.model == "stable") cov.model <- "powered.exponential" : 
	 the condition has length > 1 and only the first element will be used
Warning in if (cov.model == "stable") cov.model <- "powered.exponential" : 
	 the condition has length > 1 and only the first element will be used
> ## exponential and spherical structures
> lines.variomodel(seq(0,1,l=101), cov.model=c("exp", "sph"),
+                  cov.pars=rbind(c(0.6,0.15), c(0.4,0.75)), nug=0, col=3)
Warning in if (cov.model == "stable") cov.model <- "powered.exponential" : 
	 the condition has length > 1 and only the first element will be used
Warning in if (cov.model == "stable") cov.model <- "powered.exponential" : 
	 the condition has length > 1 and only the first element will be used
> 
> 
> 
> cleanEx(); ..nameEx <- "lines.variomodel.grf"
> 
> ### * lines.variomodel.grf
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lines.variomodel.grf
> ### Title: Lines with True Variogram for Simulated Data
> ### Aliases: lines.variomodel.grf
> ### Keywords: spatial aplot
> 
> ### ** Examples
> 
> sim <- grf(100, cov.pars=c(1, .25)) # simulates data
grf: simulation(s) on randomly chosen locations with  100  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0 
grf: covariance model 1 is: exponential(sigmasq=1, phi=0.25)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 1 
> plot(variog(sim, max.dist=1))       # plot empirical variogram
variog: computing omnidirectional variogram
> 
> 
> 
> cleanEx(); ..nameEx <- "lines.variomodel.likGRF"
> 
> ### * lines.variomodel.likGRF
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lines.variomodel.likGRF
> ### Title: Adds a Variogram Line to a Variogram Plot
> ### Aliases: lines.variomodel.likGRF
> ### Keywords: spatial aplot
> 
> ### ** Examples
> 
> data(s100)
> # compute and plot empirical variogram
> vario <- variog(s100, max.dist = 1)
variog: computing omnidirectional variogram
> plot(vario)
> # estimate parameters
> vario.ml <- likfit(s100, ini = c(1, .3), fix.nugget = TRUE)
---------------------------------------------------------------
likfit: likelihood maximisation using the function optimize.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optimize.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.
> # adds fitted model to the plot  
> lines(vario.ml)
> 
> 
> 
> cleanEx(); ..nameEx <- "lines.variomodel.variofit"
> 
> ### * lines.variomodel.variofit
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lines.variomodel.variofit
> ### Title: Adds a Line with a Fitted Variogram Model to a Variogram Plot
> ### Aliases: lines.variomodel.variofit
> ### Keywords: spatial aplot
> 
> ### ** Examples
> 
> data(s100)
> # compute and plot empirical variogram
> vario <- variog(s100, max.dist = 1)
variog: computing omnidirectional variogram
> plot(vario)
> # estimate parameters
> vario.wls <- variofit(vario, ini = c(1, .3), fix.nugget = TRUE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
> # adds fitted model to the plot  
> lines(vario.wls)
> 
> 
> 
> cleanEx(); ..nameEx <- "locations.inside"
> 
> ### * locations.inside
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: locations.inside
> ### Title: Select prediction locations inside borders
> ### Aliases: locations.inside .check.coords
> ### Keywords: spatial
> 
> ### ** Examples
> 
> if(require(splancs)){
+   data(parana)
+   gr <- expand.grid(seq(150, 800, by=20), seq(70, 500, by=20))
+   plot(gr, asp=1, pch="+")
+   polygon(parana$borders)
+   gr.in <- locations.inside(gr, parana$borders)
+   points(gr.in, col=2, pch="+")
+ }
Loading required package: splancs

Spatial Point Pattern Analysis Code in S-Plus

 Version 2 - Spatial and Space-Time analysis
> 
> 
> 
> cleanEx(); ..nameEx <- "loglik.GRF"
> 
> ### * loglik.GRF
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: loglik.GRF
> ### Title: Log-Likelihood for a Gaussian Random Field
> ### Aliases: loglik.GRF
> ### Keywords: spatial
> 
> ### ** Examples
> 
> data(s100)
> loglik.GRF(s100, cov.pars=c(0.8, .25), nugget=0.2)
[1] -92.71374
> loglik.GRF(s100, cov.pars=c(0.8, .25), nugget=0.2, met="RML")
[1] -90.502
> 
> ## Computing the likelihood of a model fitted by ML
> s100.ml <- likfit(s100, ini=c(1, .5))
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optim.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.
> s100.ml
likfit: estimated model parameters:
    beta    tausq  sigmasq      phi 
"0.7766" "0.0000" "0.7517" "0.1827" 

likfit: maximised log-likelihood = -83.57
> s100.ml$loglik
[1] -83.56957
> loglik.GRF(s100, obj=s100.ml)
[1] -83.56957
> 
> ## Computing the likelihood of a variogram fitted model
> s100.v <- variog(s100, max.dist=1)
variog: computing omnidirectional variogram
> s100.vf <- variofit(s100.v, ini=c(1, .5))
variofit: weights used: npairs 
variofit: minimisation function used: optim 
> s100.vf
variofit: model parameters estimated by WLS (weighted least squares):
covariance model is: matern with fixed kappa = 0.5 (exponential)
parameter estimates:
  tausq sigmasq     phi 
 0.1848  1.4091  0.9620 

variofit: minimised weighted sum of squares = 23.1651
> loglik.GRF(s100, obj=s100.vf)
[1] -90.60537
> 
> 
> 
> cleanEx(); ..nameEx <- "matern"
> 
> ### * matern
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: matern
> ### Title: Computer Values of the Matern Correlation Function
> ### Aliases: matern
> ### Keywords: spatial
> 
> ### ** Examples
> 
> #
> # Models with fixed range and varying smoothness parameter
> #
> curve(matern(x, phi= 0.25, kappa = 0.5),from = 0, to = 1,
+       xlab = "distance", ylab = expression(rho(h)), lty = 2,
+       main=expression(paste("varying  ", kappa, "  and fixed  ", phi)))
> curve(matern(x, phi= 0.25, kappa = 1),from = 0, to = 1, add = TRUE)
> curve(matern(x, phi= 0.25, kappa = 2),from = 0, to = 1, add = TRUE,
+       lwd = 2, lty=2)
> curve(matern(x, phi= 0.25, kappa = 3),from = 0, to = 1, add = TRUE,
+       lwd = 2)
> legend(0.6,1, c(expression(kappa == 0.5), expression(kappa == 1),
+     expression(kappa == 2), expression(kappa == 3)),
+     lty=c(2,1,2,1), lwd=c(1,1,2,2))
> #
> # Correlations with equivalent "practical range"
> # and varying smoothness parameter
> #
> curve(matern(x, phi = 0.25, kappa = 0.5),from = 0, to = 1,
+       xlab = "distance", ylab = expression(gamma(h)), lty = 2,
+       main = "models with equivalent \"practical\" range")
> curve(matern(x, phi = 0.188, kappa = 1),from = 0, to = 1, add = TRUE)      
> curve(matern(x, phi = 0.14, kappa = 2),from = 0, to = 1,
+       add = TRUE, lwd=2, lty=2)      
> curve(matern(x, phi = 0.117, kappa = 2),from = 0, to = 1,
+       add = TRUE, lwd=2)      
> legend(0.4,1, c(expression(paste(kappa == 0.5, "  and  ",
+        phi == 0.250)),
+        expression(paste(kappa == 1, "  and  ", phi == 0.188)),
+        expression(paste(kappa == 2, "  and  ", phi == 0.140)),
+        expression(paste(kappa == 3, "  and  ", phi == 0.117))),
+        lty=c(2,1,2,1), lwd=c(1,1,2,2))
> 
> 
> 
> cleanEx(); ..nameEx <- "nearloc"
> 
> ### * nearloc
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: nearloc
> ### Title: Near location to a point
> ### Aliases: nearloc
> ### Keywords: spatial
> 
> ### ** Examples
> 
> set.seed(276)
> gr <- expand.grid(seq(0,1, l=11), seq(0,1, l=11))
> plot(gr, asp=1)
> pts <- matrix(runif(10), nc=2)
> points(pts, pch=19)
> near <- nearloc(points=pts, locations=gr)
> points(near, pch=19, col=2)
> rownames(near)
[1] "71" "46" "64" "38" "79"
> nearloc(points=pts, locations=gr, pos=TRUE)
[1] 71 46 64 38 79
> 
> 
> 
> cleanEx(); ..nameEx <- "parana"
> 
> ### * parana
> 
> flush(stderr()); flush(stdout())
> 
> ### Encoding: latin1
> 
> ### Name: parana
> ### Title: Rainfall Data from Parana State, Brasil
> ### Aliases: parana maijun
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(parana)
> summary(parana)
Number of data points: 143 

Coordinates summary
        east    north
min 150.1220  70.3600
max 768.5087 461.9681

Distance summary
     min      max 
  1.0000 619.4925 

Borders summary
        east    north
min 137.9873  46.7695
max 798.6256 507.9295

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  162.8   234.2   269.9   274.4   318.2   413.7 

Other elements in the geodata object
[1] "loci.paper"
> plot(parana, bor=borders)
> 
> 
> 
> cleanEx(); ..nameEx <- "pars.limits"
> 
> ### * pars.limits
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pars.limits
> ### Title: Set limits for the parameter values
> ### Aliases: pars.limits
> ### Keywords: spatial models
> 
> ### ** Examples
> 
> pars.limits(phi=c(0,10))
$phi
lower upper 
    0    10 

$sigmasq
lower upper 
    0   Inf 

$tausq.rel
lower upper 
    0   Inf 

$kappa
lower upper 
    0   Inf 

$lambda
lower upper 
   -3     3 

$psiR
lower upper 
    1   Inf 

$psiA
   lower    upper 
0.000000 6.283185 

> pars.limits(phi=c(0,10), sigmasq=c(0, 100))
$phi
lower upper 
    0    10 

$sigmasq
lower upper 
    0   100 

$tausq.rel
lower upper 
    0   Inf 

$kappa
lower upper 
    0   Inf 

$lambda
lower upper 
   -3     3 

$psiR
lower upper 
    1   Inf 

$psiA
   lower    upper 
0.000000 6.283185 

> 
> 
> 
> cleanEx(); ..nameEx <- "plot.geodata"
> 
> ### * plot.geodata
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.geodata
> ### Title: Exploratory Geostatistical Plots
> ### Aliases: plot.geodata
> ### Keywords: spatial dplot
> 
> ### ** Examples
> 
> require(geoR)
[1] TRUE
> data(s100)
> plot(s100)
> plot(s100, scatter3d=TRUE)
Loading required package: scatterplot3d
> plot(s100, qt.col=1)
> 
> data(ca20)                        
> plot(ca20, bor=borders)          # original data
> plot(ca20, trend=~altitude+area) # residuals from an external trend
> plot(ca20, trend='1st')          # residuals from a polynomial trend
> 
> data(SIC)
> plot(sic.100, bor=sic.borders)           # original data
> plot(sic.100, bor=sic.borders, lambda=0) # logarithm of the data
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.grf"
> 
> ### * plot.grf
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.grf
> ### Title: Plots Variograms for Simulated Data
> ### Aliases: plot.grf
> ### Keywords: spatial dplot
> 
> ### ** Examples
> 
> op <- par(no.readonly = TRUE)
> par(mfrow=c(2,1))
> sim1 <- grf(100, cov.pars=c(10, .25))
grf: simulation(s) on randomly chosen locations with  100  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0 
grf: covariance model 1 is: exponential(sigmasq=10, phi=0.25)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 1 
> # generates simulated data
> plot(sim1, plot.locations = TRUE)
variog: computing omnidirectional variogram
> #
> # plots the locations and the sample true variogram model
> #
> par(mfrow=c(1,1))
> sim2 <- grf(100, cov.pars=c(10, .25), nsim=10)
grf: simulation(s) on randomly chosen locations with  100  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0 
grf: covariance model 1 is: exponential(sigmasq=10, phi=0.25)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 10 
> # generates 10 simulated data
> plot(sim1)
variog: computing omnidirectional variogram
> # plots sample variograms for all simulations with the true model
> par(op)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "plot.krige.bayes"
> 
> ### * plot.krige.bayes
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.krige.bayes
> ### Title: Plots Prior and/or Posterior Distributions
> ### Aliases: plot.krige.bayes
> ### Keywords: spatial dplot aplot
> 
> ### ** Examples
> 
> ## See documentation for krige.bayes
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.proflik"
> 
> ### * plot.proflik
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.proflik
> ### Title: Plots Profile Likelihoods
> ### Aliases: plot.proflik .proflik.plot.aux1
> ### Keywords: spatial dplot
> 
> ### ** Examples
> 
> # see examples in the documentation for the function proflik()
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.variog4"
> 
> ### * plot.variog4
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.variog4
> ### Title: Plot Directional Variograms
> ### Aliases: plot.variog4
> ### Keywords: spatial dplot
> 
> ### ** Examples
> 
> if(is.R()) data(s100)
> s100.v4 <- variog4(s100, max.dist=1)
variog: computing variogram for direction = 0 degrees (0 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing variogram for direction = 45 degrees (0.785 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing variogram for direction = 90 degrees (1.571 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing variogram for direction = 135 degrees (2.356 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing omnidirectional variogram
> # Plotting variograms for the four directions
> plot(s100.v4)
> # changing plot options
> plot(s100.v4, lwd=2)
> plot(s100.v4, lty=1, col=c("darkorange", "darkblue", "darkgreen","darkviolet"))
> plot(s100.v4, lty=1, lwd=2)
> # including the omnidirectional variogram
> plot(s100.v4, omni=TRUE)
> # variograms on different plots
> plot(s100.v4, omni=TRUE, same=FALSE)
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.variogram"
> 
> ### * plot.variogram
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.variogram
> ### Title: Plot Empirical Variogram
> ### Aliases: plot.variogram
> ### Keywords: spatial dplot
> 
> ### ** Examples
> 
> op <- par(no.readonly = TRUE)
> sim <- grf(100, cov.pars=c(1, .2)) # simulates data
grf: simulation(s) on randomly chosen locations with  100  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0 
grf: covariance model 1 is: exponential(sigmasq=1, phi=0.2)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 1 
> vario <- variog(sim, max.dist=1)   # computes sample variogram
variog: computing omnidirectional variogram
> par(mfrow=c(2,2))
> plot(vario)                     # the sample variogram
> plot(vario, scaled = TRUE)      # the scaled sample variogram
> plot(vario, max.dist = 1)       # limiting the maximum distance
> plot(vario, pts.range = c(1,3)) # points sizes proportional to number of pairs
> par(op)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "plot.xvalid"
> 
> ### * plot.xvalid
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.xvalid
> ### Title: Plot Cross-Validation Results
> ### Aliases: plot.xvalid
> ### Keywords: spatial dplot
> 
> ### ** Examples
> 
> if(is.R()) data(s100)
> wls <- variofit(variog(s100, max.dist = 1), ini = c(.5, .5), fix.n = TRUE)
variog: computing omnidirectional variogram
variofit: weights used: npairs 
variofit: minimisation function used: optim 
> xvl <- xvalid(s100, model = wls)
xvalid: number of data locations       = 100
xvalid: number of validation locations = 100
xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 
xvalid: end of cross-validation
> #
> op <- par(no.readonly = TRUE)
> par(mfcol = c(3,2))
> par(mar = c(3,3,0,1))
> par(mgp = c(2,1,0))
> plot(xvl, error = FALSE, ask = FALSE)
> plot(xvl, std.err = FALSE, ask = FALSE)
> par(op)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "points.geodata"
> 
> ### * points.geodata
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: points.geodata
> ### Title: Plots Spatial Locations and Data Values
> ### Aliases: points.geodata
> ### Keywords: spatial dplot aplot
> 
> ### ** Examples
> 
> data(s100)
> op <- par(no.readonly = TRUE)
> par(mfrow=c(2,2), mar=c(3,3,1,1), mgp = c(2,1,0))
> points(s100, xlab="Coord X", ylab="Coord Y")
> points(s100, xlab="Coord X", ylab="Coord Y", pt.divide="rank.prop")
> points(s100, xlab="Coord X", ylab="Coord Y", cex.max=1.7,
+                col=gray(seq(1, 0.1, l=100)), pt.divide="equal")
> points(s100, pt.divide="quintile", xlab="Coord X", ylab="Coord Y")
> par(op)
> 
> data(ca20)
> points(ca20, pt.div='quartile', x.leg=4900, y.leg=5850, bor=borders)
> 
> par(mfrow=c(1,2), mar=c(3,3,1,1), mgp = c(2,1,0))
> points(s100, main="Original data")
> points(s100, permute=TRUE, main="Permuting locations")
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "polygrid"
> 
> ### * polygrid
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: polygrid
> ### Title: Coordinates of Points Inside a Polygon
> ### Aliases: polygrid
> ### Keywords: spatial
> 
> ### ** Examples
> 
> if(require(splancs)){
+  poly <- matrix(c(.2, .8, .7, .1, .2, .1, .2, .7, .7, .1), ncol=2)
+  plot(0:1, 0:1, type="n")
+  lines(poly)
+  poly.in <- polygrid(seq(0,1,l=11), seq(0,1,l=11), poly, vec=TRUE)
+  points(poly.in$xy)
+ }
Loading required package: splancs

Spatial Point Pattern Analysis Code in S-Plus

 Version 2 - Spatial and Space-Time analysis
> 
> 
> 
> cleanEx(); ..nameEx <- "proflik"
> 
> ### * proflik
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: proflik
> ### Title: Computes Profile Likelihoods
> ### Aliases: proflik .proflik.aux0 .proflik.aux1 .proflik.aux10
> ###   .proflik.aux11 .proflik.aux1.1 .proflik.aux12 .proflik.aux13
> ###   .proflik.aux14 .proflik.aux15 .proflik.aux16 .proflik.aux17
> ###   .proflik.aux18 .proflik.aux19 .proflik.aux2 .proflik.aux20
> ###   .proflik.aux21 .proflik.aux21.1 .proflik.aux22 .proflik.aux23
> ###   .proflik.aux24 .proflik.aux27 .proflik.aux28 .proflik.aux30
> ###   .proflik.aux3 .proflik.aux31 .proflik.aux32 .proflik.aux33
> ###   .proflik.aux4 .proflik.aux5 .proflik.aux6 .proflik.aux7 .proflik.aux8
> ###   .proflik.aux9 .proflik.cov .proflik.lambda .proflik.main
> ### Keywords: spatial
> 
> ### ** Examples
> 
> op <- par(no.readonly=TRUE)
> data(s100)
> ml <- likfit(s100, ini=c(.5, .5), fix.nug=TRUE)
---------------------------------------------------------------
likfit: likelihood maximisation using the function optimize.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optimize.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.
> ## a first atempt to find reasonable values for the x-axis:
> prof <- proflik(ml, s100, sill.values=seq(0.5, 1.5, l=4),
+                 range.val=seq(0.1, .5, l=4))
proflik: computing profile likelihood for the sill
proflik: computing profile likelihood for the range
> par(mfrow=c(1,2))
> plot(prof)
> ## a nicer setting 
> ## Not run: 
> ##D prof <- proflik(ml, s100, sill.values=seq(0.45, 2, l=11),
> ##D                 range.val=seq(0.1, .55, l=11))
> ##D plot(prof)
> ##D ## to include 2-D profiles use:
> ##D ## (commented because this is time demanding)
> ##D #prof <- proflik(ml, s100, sill.values=seq(0.45, 2, l=11),
> ##D #                range.val=seq(0.1, .55, l=11), uni.only=FALSE)
> ##D #par(mfrow=c(2,2))
> ##D #plot(prof, nlevels=16)
> ## End(Not run)
> par(op)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "s100"
> 
> ### * s100
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: s100 and s121
> ### Title: Simulated Data-Sets which Illustrate the Usage of the Package
> ###   geoR
> ### Aliases: s100 s121
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(s100)
> plot(s100)
> data(s121)
> plot(s121, type="l")
variog: computing omnidirectional variogram
> 
> 
> 
> cleanEx(); ..nameEx <- "s256i"
> 
> ### * s256i
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: s256i
> ### Title: Simulated Data-Set which Illustrate the Usage of krige.bayes
> ### Aliases: s256i
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(s256i)
> points(s256i, pt.div="quintiles", cex.min=1, cex.max=1)
> 
> 
> 
> cleanEx(); ..nameEx <- "sample.prior"
> 
> ### * sample.prior
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sample.prior
> ### Title: Sample the prior distribution
> ### Aliases: sample.prior
> ### Keywords: spatial distribution
> 
> ### ** Examples
> 
> sample.prior(50, prior=prior.control(beta.prior = "normal", beta = .5, beta.var.std=0.1, sigmasq.prior="sc", sigmasq=1.2, df.sigmasq= 2, phi.prior="rec", phi.discrete = seq(0,2, l=21)))
          beta     sigmasq phi tausq.rel
1   0.89886163   2.6024803 0.1         0
2  -0.24408489   9.1759846 0.2         0
3   0.04407973   5.4781082 0.4         0
4   0.52738361   3.4560252 1.5         0
5  -0.09717768   2.7909610 0.1         0
6   1.02248018   8.2075105 1.4         0
7   0.01641619   1.4256574 1.7         0
8   1.64364010   4.9503164 0.6         0
9   0.36202871   0.7593258 0.5         0
10  1.03401494   1.0124395 0.1         0
11  0.27130155   3.0735244 0.1         0
12 -0.20064421   5.1928559 0.1         0
13  0.50892281   1.2357263 0.7         0
14  0.52623925   9.1204622 0.2         0
15  0.20001456   0.3187772 0.9         0
16  0.90283186   1.4614089 0.3         0
17 -0.46695764   7.4591581 0.7         0
18  0.87571825  12.5324545 2.0         0
19  0.62039902   0.5920984 0.2         0
20  0.52817573   0.4165437 0.9         0
21  0.42139241   4.3787932 1.6         0
22  0.60789889   2.9791318 0.1         0
23 -0.03849183   2.5389406 0.6         0
24  0.07140521   2.8472925 0.1         0
25 -0.65008602  10.6628527 0.1         0
26  0.84714741   0.4826844 0.2         0
27  0.78931415   0.3730962 0.1         0
28  0.59629247   1.3441412 0.2         0
29 -0.71286224   9.6775717 1.3         0
30  0.49800799   2.8620133 0.2         0
31  7.02208383 186.1476051 0.3         0
32  0.28846174   1.9774917 0.5         0
33  0.92640169   2.8557686 0.3         0
34  0.04666853   2.1662647 0.1         0
35  1.07913871   7.0575998 1.1         0
36  0.25181742   0.6741756 0.6         0
37  0.10969221   1.0041523 1.0         0
38 -0.05491361   3.3629856 0.1         0
39  0.19992265   1.1902679 0.8         0
40  0.13887872   1.5724394 0.2         0
41  0.70730978   0.7821307 1.1         0
42  0.54276118   3.8955796 0.6         0
43  0.39116288   1.1301406 0.9         0
44  0.22259288   0.6518885 0.4         0
45 -0.80114490  16.4030519 0.4         0
46  0.29339031   0.7241295 1.0         0
47  0.13310876   1.0736328 0.1         0
48  0.39653470   0.5329613 0.3         0
49  0.66567575   1.2334432 0.8         0
50  0.05363083   1.4297604 0.7         0
> 
> 
> 
> cleanEx(); ..nameEx <- "soil"
> 
> ### * soil
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: soil
> ### Title: Soil chemistry properties data set
> ### Aliases: soil
> ### Keywords: datasets spatial
> 
> ### ** Examples
> 
> data(soil)
> ctc <- as.geodata(soil, data.col=16)
> plot(ctc)
> 
> 
> 
> cleanEx(); ..nameEx <- "subarea"
> 
> ### * subarea
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: subarea
> ### Title: Selects a Subarea from a Geodata Object
> ### Aliases: subarea
> ### Keywords: spatial
> 
> ### ** Examples
> 
> foo <- matrix(c(4,6,6,4,2,2,4,4), nc=2)
> foo1 <- zoom.coords(foo, 2)
> foo1
     [,1] [,2]
[1,]    3    1
[2,]    7    1
[3,]    7    5
[4,]    3    5
> foo2 <- coords2coords(foo, c(6,10), c(6,10))
> foo2
     [,1] [,2]
[1,]    6    6
[2,]   10    6
[3,]   10   10
[4,]    6   10
> plot(1:10, 1:10, type="n")
> polygon(foo)
> polygon(foo1, lty=2)
> polygon(foo2, lwd=2)
> arrows(foo[,1], foo[,2],foo1[,1],foo1[,2], lty=2)
> arrows(foo[,1], foo[,2],foo2[,1],foo2[,2])
> legend(1,10, c("foo", "foo1 (zoom.coords)", "foo2 (coords2coords)"), lty=c(1,2,1), lwd=c(1,1,2), cex=1.7)
> 
> ## "zooming" part of The Gambia map
> data(gambia)
> gb <- gambia.borders/1000
> gd <- gambia[,1:2]/1000
> plot(gb, ty="l", asp=1, xlab="W-E (kilometres)", ylab="N-S (kilometres)")
> points(gd, pch=19, cex=0.5)
> r1b <- gb[gb[,1] < 420,]
> rc1 <- rect.coords(r1b, lty=2)
> 
> r1bn <- zoom.coords(r1b, 1.8, xoff=90, yoff=-90)
> rc2 <- rect.coords(r1bn, xz=1.05)
> segments(rc1[c(1,3),1],rc1[c(1,3),2],rc2[c(1,3),1],rc2[c(1,3),2], lty=3)
> 
> lines(r1bn)
> r1d <- gd[gd[,1] < 420,]
> r1dn <- zoom.coords(r1d, 1.7, xlim.o=range(r1b[,1],na.rm=TRUE), ylim.o=range(r1b[,2],na.rm=TRUE), xoff=90, yoff=-90)
> points(r1dn, pch=19, cex=0.5)
> text(450,1340, "Western Region", cex=1.5)
> 
> if(require(geoRglm)){
+ data(rongelap)
+ points(rongelap, bor=borders)
+ ## zooming the western area
+ rongwest <- subarea(rongelap, xlim=c(-6300, -4800))
+ points(rongwest, bor=borders)
+ ## now zooming in the same plot
+ points(rongelap, bor=borders)
+ rongwest.z <- zoom.coords(rongwest, xzoom=3, xoff=2000, yoff=3000)
+ points(rongwest.z, bor=borders, add=TRUE)
+ rect.coords(rongwest$sub, quiet=TRUE)
+ rect.coords(rongwest.z$sub, quiet=TRUE)
+ }
Loading required package: geoRglm
-----------------------------------------------------------
geoRglm - a package for generalised linear spatial models
geoRglm version 0.8-11 (2005-06-07) is now loaded
-----------------------------------------------------------

> 
> 
> 
> cleanEx(); ..nameEx <- "summary.geodata"
> 
> ### * summary.geodata
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.geodata
> ### Title: Summaries for geodata object
> ### Aliases: summary.geodata print.summary.geodata
> ### Keywords: univar spatial
> 
> ### ** Examples
> 
> data(s100)
> summary(s100)
Number of data points: 100 

Coordinates summary
        Coord.X    Coord.Y
min 0.005638006 0.01091027
max 0.983920544 0.99124979

Distance summary
        min         max 
0.007640962 1.278175109 

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-1.1680  0.2730  1.1050  0.9307  1.6100  2.8680 

Other elements in the geodata object
[1] "cov.model" "nugget"    "cov.pars"  "kappa"     "lambda"   
> 
> data(ca20)
> summary(ca20)
Number of data points: 178 

Coordinates summary
    east north
min 4957  4829
max 5961  5720

Distance summary
       min        max 
  43.01163 1138.11774 

Borders summary
    east north
min 4920  4800
max 5990  5800

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  21.00   43.00   50.50   50.68   58.00   78.00 

Covariates summary
    altitude     area   
 Min.   :3.300   1: 14  
 1st Qu.:5.200   2: 48  
 Median :5.650   3:116  
 Mean   :5.524          
 3rd Qu.:6.000          
 Max.   :6.600          

Other elements in the geodata object
[1] "reg1" "reg2" "reg3"
> 
> 
> 
> cleanEx(); ..nameEx <- "summary.likGRF"
> 
> ### * summary.likGRF
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.likGRF
> ### Title: Summarizes Parameter Estimation Results for Gaussian Random
> ###   Fields
> ### Aliases: summary.likGRF print.summary.likGRF print.likGRF
> ### Keywords: spatial print
> 
> ### ** Examples
> 
> # See examples for the function likfit()
> 
> 
> 
> cleanEx(); ..nameEx <- "summary.variofit"
> 
> ### * summary.variofit
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.variofit
> ### Title: Summarize Results of Variogram Estimation
> ### Aliases: summary.variofit print.summary.variofit print.variofit
> ### Keywords: spatial
> 
> ### ** Examples
> 
> data(s100)
> s100.vario <- variog(s100, max.dist=1)
variog: computing omnidirectional variogram
> wls <- variofit(s100.vario, ini=c(.5, .5), fix.nugget = TRUE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
> wls
variofit: model parameters estimated by WLS (weighted least squares):
covariance model is: matern with fixed kappa = 0.5 (exponential)
fixed value for tausq =  0 
parameter estimates:
sigmasq     phi 
 1.1894  0.4721 

variofit: minimised weighted sum of squares = 28.7342
> summary(wls)
$pmethod
[1] "WLS (weighted least squares)"

$cov.model
[1] "matern"

$spatial.component
  sigmasq       phi 
1.1894204 0.4720877 

$spatial.component.extra
kappa 
  0.5 

$nugget.component
tausq 
    0 

$fix.nugget
[1] TRUE

$fix.kappa
[1] TRUE

$sum.of.squares
   value 
28.73421 

$estimated.pars
  sigmasq       phi 
1.1894204 0.4720877 

$weights
[1] "npairs"

$call
variofit(vario = s100.vario, ini.cov.pars = c(0.5, 0.5), fix.nugget = TRUE)

attr(,"class")
[1] "summary.variomodel"
> 
> 
> 
> cleanEx(); ..nameEx <- "tce"
> 
> ### * tce
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: tce
> ### Title: TCE concentrations in groundwater in a vertical cross section
> ### Aliases: tce
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(tce)
> summary(tce)
Number of data points: 56 

Coordinates summary
      x   y
min   0 -90
max 370 -45

Distance summary
     min      max 
  5.0000 371.6517 

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    2.5    41.5   434.0  4719.0  2843.0 67700.0 
> summary(tce, lambda=0)
Number of data points: 56 

Coordinates summary
      x   y
min   0 -90
max 370 -45

Distance summary
     min      max 
  5.0000 371.6517 

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    2.5    41.5   434.0  4719.0  2843.0 67700.0 

Transformed Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.9163  3.7250  6.0520  5.8550  7.9370 11.1200 
> plot(tce)
> points(tce)
> points(tce, lambda=0)
> 
> 
> 
> cleanEx(); ..nameEx <- "trend.spatial"
> 
> ### * trend.spatial
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: trend.spatial
> ### Title: Builds the Trend Matrix
> ### Aliases: trend.spatial
> ### Keywords: spatial
> 
> ### ** Examples
> 
> data(SIC)
> # a first order polynomial trend
> trend.spatial("1st", sic.100)[1:5,]
     [,1]     [,2]      [,3]
[1,]    1 29.52739  80.71854
[2,]    1 33.77939  99.52954
[3,]    1 46.80639 102.58454
[4,]    1 48.71439 121.45354
[5,]    1 49.31639 113.65554
> # a second order polynomial trend
> trend.spatial("2nd", sic.100)[1:5,]
     [,1]     [,2]      [,3]      [,4]      [,5]     [,6]
[1,]    1 29.52739  80.71854  871.8668  6515.483 2383.408
[2,]    1 33.77939  99.52954 1141.0473  9906.130 3362.047
[3,]    1 46.80639 102.58454 2190.8382 10523.588 4801.612
[4,]    1 48.71439 121.45354 2373.0919 14750.963 5916.535
[5,]    1 49.31639 113.65554 2432.1064 12917.582 5605.081
> # a trend with a covariate
> trend.spatial(~altitude, sic.100)[1:5,]
     [,1] [,2]
[1,]    1  682
[2,]    1  813
[3,]    1  436
[4,]    1  833
[5,]    1  579
> # a first degree trend plus a covariate
> trend.spatial(~coords+altitude, sic.100)[1:5,]
     [,1]     [,2]      [,3] [,4]
[1,]    1 29.52739  80.71854  682
[2,]    1 33.77939  99.52954  813
[3,]    1 46.80639 102.58454  436
[4,]    1 48.71439 121.45354  833
[5,]    1 49.31639 113.65554  579
> # with produces the same as
> trend.spatial("1st", sic.100, add=~altitude)[1:5,]
     [,1]     [,2]      [,3] [,4]
[1,]    1 29.52739  80.71854  682
[2,]    1 33.77939  99.52954  813
[3,]    1 46.80639 102.58454  436
[4,]    1 48.71439 121.45354  833
[5,]    1 49.31639 113.65554  579
> # and yet another exemple
> trend.spatial("2nd", sic.100, add=~altitude)[1:5,]
     [,1]     [,2]      [,3]      [,4]      [,5]     [,6] [,7]
[1,]    1 29.52739  80.71854  871.8668  6515.483 2383.408  682
[2,]    1 33.77939  99.52954 1141.0473  9906.130 3362.047  813
[3,]    1 46.80639 102.58454 2190.8382 10523.588 4801.612  436
[4,]    1 48.71439 121.45354 2373.0919 14750.963 5916.535  833
[5,]    1 49.31639 113.65554 2432.1064 12917.582 5605.081  579
> 
> 
> 
> cleanEx(); ..nameEx <- "variofit"
> 
> ### * variofit
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: variofit
> ### Title: Variogram Based Parameter Estimation
> ### Aliases: variofit .loss.vario
> ### Keywords: spatial
> 
> ### ** Examples
> 
> data(s100)
> vario100 <- variog(s100, max.dist=1)
variog: computing omnidirectional variogram
> ini.vals <- expand.grid(seq(0,1,l=5), seq(0,1,l=5))
> ols <- variofit(vario100, ini=ini.vals, fix.nug=TRUE, wei="equal")
variofit: weights used: equal 
variofit: minimisation function used: nls 
variofit: searching for best initial value ... selected values:
              sigmasq phi    tausq kappa
initial.value "1"     "0.25" "0"   "0.5"
status        "est"   "est"  "fix" "fix"
loss value: 0.167887026209413 
> summary(ols)
$pmethod
[1] "OLS (ordinary least squares)"

$cov.model
[1] "matern"

$spatial.component
  sigmasq       phi 
1.1066543 0.4003639 

$spatial.component.extra
kappa.NA 
     0.5 

$nugget.component
tausq.NA 
       0 

$fix.nugget
[1] TRUE

$fix.kappa
[1] TRUE

$sum.of.squares
    value 
0.1027018 

$estimated.pars
  sigmasq       phi 
1.1066543 0.4003639 

$weights
[1] "equal"

$call
variofit(vario = vario100, ini.cov.pars = ini.vals, fix.nugget = TRUE, 
    weights = "equal")

attr(,"class")
[1] "summary.variomodel"
> wls <- variofit(vario100, ini=ini.vals, fix.nug=TRUE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
variofit: searching for best initial value ... selected values:
              sigmasq phi    tausq kappa
initial.value "1"     "0.25" "0"   "0.5"
status        "est"   "est"  "fix" "fix"
loss value: 68.37535548371 
> summary(wls)
$pmethod
[1] "WLS (weighted least squares)"

$cov.model
[1] "matern"

$spatial.component
  sigmasq       phi 
1.1894198 0.4720876 

$spatial.component.extra
kappa.NA 
     0.5 

$nugget.component
tausq.NA 
       0 

$fix.nugget
[1] TRUE

$fix.kappa
[1] TRUE

$sum.of.squares
   value 
28.73421 

$estimated.pars
  sigmasq       phi 
1.1894198 0.4720876 

$weights
[1] "npairs"

$call
variofit(vario = vario100, ini.cov.pars = ini.vals, fix.nugget = TRUE)

attr(,"class")
[1] "summary.variomodel"
> plot(vario100)
> lines(wls)
> lines(ols, lty=2)
> 
> ## Don't show:
> vr <- variog(s100, max.dist=1)
variog: computing omnidirectional variogram
> ## OLS#
> o1 <- variofit(vr, ini = c(.5, .5), fix.nug=TRUE, wei = "equal")
variofit: weights used: equal 
variofit: minimisation function used: nls 
> o2 <- variofit(vr, ini = c(.5, .5), wei = "equal")
variofit: weights used: equal 
variofit: minimisation function used: nls 
> o3 <- variofit(vr, ini = c(.5, .5), fix.nug=TRUE,
+       fix.kappa = FALSE, wei = "equal")
variofit: weights used: equal 
variofit: minimisation function used: nls 
> #o4 <- variofit(vr, ini = c(.5, .5), fix.kappa = FALSE, wei = "equal")
> ## WLS
> w1 <- variofit(vr, ini = c(.5, .5), fix.nug=TRUE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
> w2 <- variofit(vr, ini = c(.5, .5))
variofit: weights used: npairs 
variofit: minimisation function used: optim 
> w3 <- variofit(vr, ini = c(.5, .5), fix.nug=TRUE, fix.kappa = FALSE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
> w4 <- variofit(vr, ini = c(.5, .5), fix.kappa = FALSE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
> ## End Don't show
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "variog"
> 
> ### * variog
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: variog
> ### Title: Compute Empirical Variograms
> ### Aliases: variog .rfm.bin .define.bins
> ### Keywords: spatial smooth robust
> 
> ### ** Examples
> 
> # Loading data:
> data(s100)
> #
> # computing variograms:
> #
> # binned variogram
> vario.b <- variog(s100, max.dist=1)
variog: computing omnidirectional variogram
> # variogram cloud
> vario.c <- variog(s100, max.dist=1, op="cloud")
variog: computing omnidirectional variogram
> #binned variogram and stores the cloud
> vario.bc <- variog(s100, max.dist=1, bin.cloud=TRUE)
variog: computing omnidirectional variogram
> # smoothed variogram
> vario.s <- variog(s100, max.dist=1, op="sm", band=0.2)
variog: computing omnidirectional variogram
> #
> #
> # plotting the variograms:
> par(mfrow=c(2,2))
> plot(vario.b, main="binned variogram") 
> plot(vario.c, main="variogram cloud")
> plot(vario.bc, bin.cloud=TRUE, main="clouds for binned variogram")  
> plot(vario.s, main="smoothed variogram") 
> 
> # computing a directional variogram
> vario.0 <- variog(s100, max.dist=1, dir=0, tol=pi/8)
variog: computing variogram for direction = 0 degrees (0 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
> plot(vario.b, type="l", lty=2)
> lines(vario.0)
> legend(0, 1.2, legend=c("omnidirectional", expression(0 * degree)), lty=c(2,1))
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "variog.mc.env"
> 
> ### * variog.mc.env
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: variog.mc.env
> ### Title: Envelops for Empirical Variograms Based on Permutation
> ### Aliases: variog.mc.env
> ### Keywords: spatial nonparametric
> 
> ### ** Examples
> 
> if(is.R()) data(s100)
> s100.vario <- variog(s100, max.dist=1)
variog: computing omnidirectional variogram
> s100.env <- variog.mc.env(s100, obj.var = s100.vario)
variog.env: generating 99 simulations by permutating data values
variog.env: computing the empirical variogram for the 99 simulations
variog.env: computing the envelops
> plot(s100.vario, envelope = s100.env)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "variog.model.env"
> 
> ### * variog.model.env
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: variog.model.env
> ### Title: Envelops for Empirical Variograms Based on Model Parameters
> ### Aliases: variog.model.env
> ### Keywords: spatial
> 
> ### ** Examples
> 
> if(is.R()) data(s100)
> s100.ml <- likfit(s100, ini = c(0.5, 0.5), fix.nugget = TRUE)
---------------------------------------------------------------
likfit: likelihood maximisation using the function optimize.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optimize.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.
> s100.vario <- variog(s100, max.dist = 1)
variog: computing omnidirectional variogram
> s100.env <- variog.model.env(s100, obj.v = s100.vario,
+                              model.pars = s100.ml)
variog.env: generating 99 simulations (with  100 points each) using the function grf
variog.env: adding the mean or trend
variog.env: computing the empirical variogram for the 99 simulations
variog.env: computing the envelops
> plot(s100.vario, env = s100.env)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "variog4"
> 
> ### * variog4
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: variog4
> ### Title: Computes Directional Variograms
> ### Aliases: variog4
> ### Keywords: spatial
> 
> ### ** Examples
> 
> data(s100)
> var4 <- variog4(s100, max.dist=1)
variog: computing variogram for direction = 0 degrees (0 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing variogram for direction = 45 degrees (0.785 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing variogram for direction = 90 degrees (1.571 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing variogram for direction = 135 degrees (2.356 radians)
        tolerance angle = 22.5 degrees (0.393 radians)
variog: computing omnidirectional variogram
> plot(var4)
> 
> 
> 
> cleanEx(); ..nameEx <- "wo"
> 
> ### * wo
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: wo
> ### Title: Kriging example data from Webster and Oliver
> ### Aliases: wo
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(wo)
> attach(wo)
> par(mfrow=c(1,2))
> plot(c(-10,130), c(-10,130), ty="n", asp=1)
> points(rbind(coords, x1))
> text(coords[,1], 5+coords[,2], format(data))
> text(x1[1]+5, x1[2]+5, "?", col=2)
> plot(c(-10,130), c(-10,130), ty="n", asp=1)
> points(rbind(coords, x2))
> text(coords[,1], 5+coords[,2], format(data))
> text(x2[1]+5, x2[2]+5, "?", col=2)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "wolfcamp"
> 
> ### * wolfcamp
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: wolfcamp
> ### Title: Wolfcamp Aquifer Data
> ### Aliases: wolfcamp wolf
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(wolfcamp)
> summary(wolfcamp)
Number of data points: 85 

Coordinates summary
      Coord.X   Coord.Y
min -233.7217 -145.7884
max  181.5314  136.4061

Distance summary
        min         max 
  0.3669819 436.2067085 

Data summary
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  312.1   471.8   547.7   610.3   774.2  1088.0 
> plot(wolfcamp)
> 
> 
> 
> cleanEx(); ..nameEx <- "xvalid"
> 
> ### * xvalid
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: xvalid
> ### Title: Cross-validation using kriging
> ### Aliases: xvalid summary.xvalid print.summary.xvalid
> ### Keywords: spatial
> 
> ### ** Examples
> 
> data(s100)
> #
> # Maximum likelihood estimation
> #
> s100.ml <- likfit(s100, ini = c(.5, .5), fix.nug = TRUE)
---------------------------------------------------------------
likfit: likelihood maximisation using the function optimize.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optimize.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.
> #
> # Weighted least squares estimation
> #
> s100.var <- variog(s100, max.dist = 1)
variog: computing omnidirectional variogram
> s100.wls <- variofit(s100.var, ini = c(.5, .5), fix.nug = TRUE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
> #
> # Now, performing cross-validation without reestimating the model
> #
> s100.xv.ml <- xvalid(s100, model = s100.ml)
xvalid: number of data locations       = 100
xvalid: number of validation locations = 100
xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 
xvalid: end of cross-validation
> s100.xv.wls <- xvalid(s100, model = s100.wls)
xvalid: number of data locations       = 100
xvalid: number of validation locations = 100
xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 
xvalid: end of cross-validation
> ##
> ## Plotting results
> ##
> par.ori <- par(no.readonly = TRUE)
> ##
> par(mfcol=c(5,2), mar=c(2.3,2.3,.5,.5), mgp=c(1.3, .6, 0))
> plot(s100.xv.ml)
> par(mfcol=c(5,2))
> plot(s100.xv.wls)
> ##
> par(par.ori)
> #
> ## Don't show:
> set.seed(234)
> data(s100)
>   vr <- variog(s100, max.dist=1)
variog: computing omnidirectional variogram
>   ## OLS#
>   o1 <- variofit(vr, ini=c(.5, .5), fix.nug=TRUE, wei="equal")
variofit: weights used: equal 
variofit: minimisation function used: nls 
>   xvo1 <- xvalid(s100, model=o1, variog.obj=vr, loc=sample(1:100,5))
xvalid: number of data locations       = 100
xvalid: number of validation locations = 5
xvalid: performing cross-validation at location ... 75, 78, 2, 76, 7, 
xvalid: end of cross-validation
>   o2 <- variofit(vr, ini=c(.5, .5), wei = "equal")
variofit: weights used: equal 
variofit: minimisation function used: nls 
>   xvo2 <- xvalid(s100, model=o2, variog.obj=vr, loc=sample(1:100,5))
xvalid: number of data locations       = 100
xvalid: number of validation locations = 5
xvalid: performing cross-validation at location ... 65, 93, 71, 90, 28, 
xvalid: end of cross-validation
>   o3 <- variofit(vr, ini=c(.5, .5), fix.nug=TRUE, fix.kappa = FALSE, wei = "equal")
variofit: weights used: equal 
variofit: minimisation function used: nls 
>   xvo3 <- xvalid(s100, model=o3, variog.obj=vr, loc=sample(1:100,5))
xvalid: number of data locations       = 100
xvalid: number of validation locations = 5
xvalid: performing cross-validation at location ... 56, 55, 58, 57, 1, 
xvalid: end of cross-validation
>   #o4 <- variofit(vr, ini=c(.5, .5), fix.kappa = FALSE, wei = "equal")
>   #xvo4 <- xvalid(s100, model=o4, variog.obj=vr, loc=sample(1:100,5))
>   ## WLS
>   w1 <- variofit(vr, ini=c(.5, .5), fix.nug=TRUE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
>   xvw1 <- xvalid(s100, model=w1, variog.obj=vr, loc=sample(1:100,5))
xvalid: number of data locations       = 100
xvalid: number of validation locations = 5
xvalid: performing cross-validation at location ... 45, 32, 73, 14, 84, 
xvalid: end of cross-validation
>   w2 <- variofit(vr, ini=c(.5, .5))
variofit: weights used: npairs 
variofit: minimisation function used: optim 
>   xvw2 <- xvalid(s100, model=w2, variog.obj=vr, loc=sample(1:100,5))
xvalid: number of data locations       = 100
xvalid: number of validation locations = 5
xvalid: performing cross-validation at location ... 53, 58, 85, 60, 48, 
xvalid: end of cross-validation
>   w3 <- variofit(vr, ini=c(.5, .5), fix.nug=TRUE, fix.kappa = FALSE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
>   xvw3 <- xvalid(s100, model=w3, variog.obj=vr, loc=sample(1:100,5))
xvalid: number of data locations       = 100
xvalid: number of validation locations = 5
xvalid: performing cross-validation at location ... 38, 69, 19, 82, 85, 
xvalid: end of cross-validation
>   w4 <- variofit(vr, ini=c(.5, .5), fix.kappa = FALSE)
variofit: weights used: npairs 
variofit: minimisation function used: optim 
>   xvw4 <- xvalid(s100, model=w4, variog.obj=vr, loc=sample(1:100,5))
xvalid: number of data locations       = 100
xvalid: number of validation locations = 5
xvalid: performing cross-validation at location ... 62, 25, 18, 69, 52, 
xvalid: end of cross-validation
>   # ML
>   m1 <- likfit(s100, ini=c(.5,.5), fix.nug=FALSE, cov.model="mat", kappa=1)
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optim.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.
>   xvm1 <- xvalid(s100, model=m1, loc=sample(1:100,5))
xvalid: number of data locations       = 100
xvalid: number of validation locations = 5
xvalid: performing cross-validation at location ... 69, 70, 16, 50, 35, 
xvalid: end of cross-validation
>   ap <- grf(10, cov.pars=c(1, .25))
grf: simulation(s) on randomly chosen locations with  10  points
grf: process with  1  covariance structure(s)
grf: nugget effect is: tausq= 0 
grf: covariance model 1 is: exponential(sigmasq=1, phi=0.25)
grf: decomposition algorithm used is:  cholesky 
grf: End of simulation procedure. Number of realizations: 1 
>   xvm2 <- xvalid(s100, model=m1, locations.xvalid=ap$coords, data.xvalid=ap$data)
xvalid: cross-validation to be performed on locations provided by the user
xvalid: number of data locations       = 100
xvalid: number of validation locations = 10
xvalid: performing validation at the locations provided
xvalid: end of cross-validation
>   bor <- s100$coords[chull(s100$coords),]
>   par(mfcol=c(5,2),mar=c(3,3,1,0),mgp=c(2,.5,0))  
>   plot(xvm2, borders=bor)
>   # RML
>   rm1 <- likfit(s100, ini=c(.5,.5), fix.nug=FALSE, met = "REML", cov.model="mat", kappa=1)
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
         arguments for the maximisation function.
        For further details see documentation for optim.
likfit: It is highly advisable to run this function several
        times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.
>   xvrm1 <- xvalid(s100, model=m1, loc=sample(1:100,5))
xvalid: number of data locations       = 100
xvalid: number of validation locations = 5
xvalid: performing cross-validation at location ... 59, 54, 33, 82, 90, 
xvalid: end of cross-validation
> ## End Don't show
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  71.08 19.29 92.86 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')