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("sgeostat-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('sgeostat')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "est.variogram"
> 
> ### * est.variogram
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: est.variogram
> ### Title: Variogram Estimator
> ### Aliases: est.variogram
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> # prepare variables from other example pages if they are not already there:
> if(length(ls(pat="maas.pair"))==0){example(pair)}

pair> if (length(ls(pat = "maas.point")) == 0) {
    example(point)
}

point> data(maas)

point> maas.point <- point(maas)

pair> maas.pair <- pair(maas.point, num.lags = 10, maxdist = 2000)
..........................................................................................................................................................

pair> maas.pair25 <- pair(maas.point, num.lags = 10, type = "anisotropic", 
    theta = 25, maxdist = 500)
..............................
NOTE: Number of pairs in lag  1 :  0 
NOTE: Number of pairs in lag  2 :  1 
NOTE: Number of pairs in lag  3 :  3 
NOTE: Number of pairs in lag  4 :  12 
NOTE: Number of pairs in lag  5 :  10 
NOTE: Number of pairs in lag  6 :  12 
NOTE: Number of pairs in lag  7 :  10 
NOTE: Number of pairs in lag  8 :  13 
NOTE: Number of pairs in lag  9 :  10 
NOTE: Number of pairs in lag  10 :  19 

> ## End Don't show
> maas.v<-est.variogram(maas.point,maas.pair,'zinc') 
> 
> 
> 
> cleanEx(); ..nameEx <- "fit.trend"
> 
> ### * fit.trend
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fit.trend
> ### Title: Fit polynomial trend functions
> ### Aliases: fit.trend
> ### Keywords: spatial spatial
> 
> ### ** Examples
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "fit.variogram"
> 
> ### * fit.variogram
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fit.variogram
> ### Title: Variogram Model Fit
> ### Aliases: fit.variogram fit.exponential fit.wave fit.gaussian
> ###   fit.spherical fit.linear
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> # prepare variables from other example pages if they are not already there:
> if(length(ls(pat="maas.v"))==0){example(est.variogram)}

est.vr> if (length(ls(pat = "maas.pair")) == 0) {
    example(pair)
}

pair> if (length(ls(pat = "maas.point")) == 0) {
    example(point)
}

point> data(maas)

point> maas.point <- point(maas)

pair> maas.pair <- pair(maas.point, num.lags = 10, maxdist = 2000)
..........................................................................................................................................................

pair> maas.pair25 <- pair(maas.point, num.lags = 10, type = "anisotropic", 
    theta = 25, maxdist = 500)
..............................
NOTE: Number of pairs in lag  1 :  0 
NOTE: Number of pairs in lag  2 :  1 
NOTE: Number of pairs in lag  3 :  3 
NOTE: Number of pairs in lag  4 :  12 
NOTE: Number of pairs in lag  5 :  10 
NOTE: Number of pairs in lag  6 :  12 
NOTE: Number of pairs in lag  7 :  10 
NOTE: Number of pairs in lag  8 :  13 
NOTE: Number of pairs in lag  9 :  10 
NOTE: Number of pairs in lag  10 :  19 


est.vr> maas.v <- est.variogram(maas.point, maas.pair, "zinc")
> ## End Don't show
> #
> # automatic fit:
> #
> maas.vmod<-fit.gaussian(maas.v,c0=60000,cg=110000,ag=800,plot.it=TRUE,
+ iterations=30)
Initial parameter estimates:  60000 110000 800 
Iteration: 1 
Gradient vector:  5499.075 -37894.01 -497.3869 
New parameter estimates:  65499.08 72106 302.6131 

rse.dif =  249665571 (rse = 249665571 )  ;  parm.dist =  38294.16 

Iteration: 2 
Gradient vector:  1877.678 14818.09 180.7598 
New parameter estimates:  67376.75 86924.08 483.373 

rse.dif =  1099461343 (rse = 1349126914 )  ;  parm.dist =  14937.67 

Iteration: 3 
Gradient vector:  -6669.642 8207.31 -37.26603 
New parameter estimates:  60707.11 95131.4 446.1069 

rse.dif =  -32082819 (rse = 1317044095 )  ;  parm.dist =  10575.7 

Iteration: 4 
Gradient vector:  909.8936 -484.5986 13.38276 
New parameter estimates:  61617 94646.8 459.4897 

rse.dif =  17615795 (rse = 1334659891 )  ;  parm.dist =  1030.981 

Iteration: 5 
Gradient vector:  -290.5093 194.8503 -2.573651 
New parameter estimates:  61326.5 94841.64 456.916 

rse.dif =  -1765886 (rse = 1332894005 )  ;  parm.dist =  349.8127 

Iteration: 6 
Gradient vector:  56.4257 -33.54592 0.5844782 
New parameter estimates:  61382.92 94808.1 457.5005 

rse.dif =  690040.2 (rse = 1333584045 )  ;  parm.dist =  65.64701 

Iteration: 7 
Gradient vector:  -12.63262 7.648973 -0.1273216 
New parameter estimates:  61370.29 94815.75 457.3732 

rse.dif =  -140598.6 (rse = 1333443446 )  ;  parm.dist =  14.76841 

Iteration: 8 
Gradient vector:  2.7564 -1.660847 0.02796193 
New parameter estimates:  61373.05 94814.08 457.4012 

rse.dif =  31404.83 (rse = 1333474851 )  ;  parm.dist =  3.218219 

Iteration: 9 
Gradient vector:  -0.6050773 0.3649581 -0.006129619 
New parameter estimates:  61372.44 94814.45 457.395 

rse.dif =  -6859.57 (rse = 1333467992 )  ;  parm.dist =  0.7066474 

Iteration: 10 
Gradient vector:  0.1326532 -0.07999283 0.001344229 
New parameter estimates:  61372.57 94814.37 457.3964 

rse.dif =  1505.506 (rse = 1333469497 )  ;  parm.dist =  0.1549114 

Iteration: 11 
Gradient vector:  -0.02909032 0.01754299 -0.0002947644 
New parameter estimates:  61372.54 94814.39 457.3961 

rse.dif =  -330.0719 (rse = 1333469167 )  ;  parm.dist =  0.0339719 

Iteration: 12 
Gradient vector:  0.006378995 -0.003846825 6.463759e-05 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  72.38277 (rse = 1333469239 )  ;  parm.dist =  0.007449417 

Iteration: 13 
Gradient vector:  -0.001398820 0.0008435545 -1.417403e-05 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -15.8723 (rse = 1333469224 )  ;  parm.dist =  0.001633549 

Iteration: 14 
Gradient vector:  0.0003067399 -0.0001849785 3.108150e-06 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  3.480562 (rse = 1333469227 )  ;  parm.dist =  0.0003582123 

Iteration: 15 
Gradient vector:  -6.726341e-05 4.0563e-05 -6.815702e-07 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -0.7632339 (rse = 1333469226 )  ;  parm.dist =  7.855054e-05 

Iteration: 16 
Gradient vector:  1.474986e-05 -8.894858e-06 1.494581e-07 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  0.1673656 (rse = 1333469227 )  ;  parm.dist =  1.722496e-05 

Iteration: 17 
Gradient vector:  -3.234433e-06 1.950522e-06 -3.277396e-08 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -0.03670073 (rse = 1333469226 )  ;  parm.dist =  3.777191e-06 

Iteration: 18 
Gradient vector:  7.092709e-07 -4.277348e-07 7.186811e-09 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  0.008048058 (rse = 1333469226 )  ;  parm.dist =  8.282952e-07 

Iteration: 19 
Gradient vector:  -1.555217e-07 9.37892e-08 -1.575938e-09 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -0.001764536 (rse = 1333469226 )  ;  parm.dist =  1.816207e-07 

Iteration: 20 
Gradient vector:  3.410008e-08 -2.055502e-08 3.456296e-10 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  0.0003864765 (rse = 1333469226 )  ;  parm.dist =  3.982316e-08 

Iteration: 21 
Gradient vector:  -7.49142e-09 4.52347e-09 -7.593299e-11 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -8.416176e-05 (rse = 1333469226 )  ;  parm.dist =  8.755045e-09 

Iteration: 22 
Gradient vector:  1.650187e-09 -9.971753e-10 1.676294e-11 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  1.716614e-05 (rse = 1333469226 )  ;  parm.dist =  1.932973e-09 

Iteration: 23 
Gradient vector:  -3.678616e-10 2.341807e-10 -3.730106e-12 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -2.622604e-06 (rse = 1333469226 )  ;  parm.dist =  4.380867e-10 

Iteration: 24 
Gradient vector:  1.004182e-10 -7.590225e-11 8.078534e-13 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -1.430511e-06 (rse = 1333469226 )  ;  parm.dist =  1.251828e-10 

Iteration: 25 
Gradient vector:  -2.121985e-11 2.328697e-11 -2.071375e-14 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  2.861023e-06 (rse = 1333469226 )  ;  parm.dist =  3.637979e-11 

Iteration: 26 
Gradient vector:  -1.061002e-11 5.692926e-12 -6.3925e-14 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -9.536743e-07 (rse = 1333469226 )  ;  parm.dist =  7.27618e-12 

Convergence achieved by sums of squares.
Final parameter estimates:  61372.55 94814.38 457.3961 

> #
> # iterations=0, means no fit, intended for "subjective" fit
> #
> maas.vmod.fixed<-fit.variogram("gaussian",maas.v,nugget=60000,sill=110000,range=800,plot.it=TRUE,iterations=0)
Convergence not achieved!
> 
> 
> 
> cleanEx(); ..nameEx <- "identify.point"
> 
> ### * identify.point
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: identify.point
> ### Title: Identify points on a Point Object
> ### Aliases: identify.point
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> # prepare variables from other example pages if they are not already there:
> if(length(ls(pat="maas.point"))==0){example(point)}

point> data(maas)

point> maas.point <- point(maas)
> ## End Don't show
> plot(maas.point)
> # use indices as labels:
> identify(maas.point)
numeric(0)
> # use values as labels:
> identify(maas.point,v="zinc")
numeric(0)
> 
> 
> 
> cleanEx(); ..nameEx <- "in.chull"
> 
> ### * in.chull
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: in.chull
> ### Title: Convex hull test
> ### Aliases: in.chull
> ### Keywords: spatial
> 
> ### ** Examples
> 
> in.chull(c(0,1),c(0,1),c(0,1,0,-1),c(-1,0,1,0))
[1]  TRUE FALSE
> # should give: TRUE FALSE
> 
> 
> 
> cleanEx(); ..nameEx <- "in.polygon"
> 
> ### * in.polygon
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: in.polygon
> ### Title: In-Polygon test
> ### Aliases: in.polygon
> ### Keywords: spatial
> 
> ### ** Examples
> 
> in.polygon(c(0,1),c(0,1),c(0,1,0,-1),c(-1,0,1,0))
[1]  TRUE FALSE
> # should give: TRUE FALSE
> 
> 
> 
> cleanEx(); ..nameEx <- "krige"
> 
> ### * krige
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: krige
> ### Title: Kriging
> ### Aliases: krige
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> # prepare variables from other example pages if they are not already there:
> if(length(ls(pat="maas.vmod"))==0){example(fit.variogram)}

ft.vrg> if (length(ls(pat = "maas.v")) == 0) {
    example(est.variogram)
}

est.vr> if (length(ls(pat = "maas.pair")) == 0) {
    example(pair)
}

pair> if (length(ls(pat = "maas.point")) == 0) {
    example(point)
}

point> data(maas)

point> maas.point <- point(maas)

pair> maas.pair <- pair(maas.point, num.lags = 10, maxdist = 2000)
..........................................................................................................................................................

pair> maas.pair25 <- pair(maas.point, num.lags = 10, type = "anisotropic", 
    theta = 25, maxdist = 500)
..............................
NOTE: Number of pairs in lag  1 :  0 
NOTE: Number of pairs in lag  2 :  1 
NOTE: Number of pairs in lag  3 :  3 
NOTE: Number of pairs in lag  4 :  12 
NOTE: Number of pairs in lag  5 :  10 
NOTE: Number of pairs in lag  6 :  12 
NOTE: Number of pairs in lag  7 :  10 
NOTE: Number of pairs in lag  8 :  13 
NOTE: Number of pairs in lag  9 :  10 
NOTE: Number of pairs in lag  10 :  19 


est.vr> maas.v <- est.variogram(maas.point, maas.pair, "zinc")

ft.vrg> maas.vmod <- fit.gaussian(maas.v, c0 = 60000, cg = 110000, 
    ag = 800, plot.it = TRUE, iterations = 30)
Initial parameter estimates:  60000 110000 800 
Iteration: 1 
Gradient vector:  5499.075 -37894.01 -497.3869 
New parameter estimates:  65499.08 72106 302.6131 

rse.dif =  249665571 (rse = 249665571 )  ;  parm.dist =  38294.16 

Iteration: 2 
Gradient vector:  1877.678 14818.09 180.7598 
New parameter estimates:  67376.75 86924.08 483.373 

rse.dif =  1099461343 (rse = 1349126914 )  ;  parm.dist =  14937.67 

Iteration: 3 
Gradient vector:  -6669.642 8207.31 -37.26603 
New parameter estimates:  60707.11 95131.4 446.1069 

rse.dif =  -32082819 (rse = 1317044095 )  ;  parm.dist =  10575.7 

Iteration: 4 
Gradient vector:  909.8936 -484.5986 13.38276 
New parameter estimates:  61617 94646.8 459.4897 

rse.dif =  17615795 (rse = 1334659891 )  ;  parm.dist =  1030.981 

Iteration: 5 
Gradient vector:  -290.5093 194.8503 -2.573651 
New parameter estimates:  61326.5 94841.64 456.916 

rse.dif =  -1765886 (rse = 1332894005 )  ;  parm.dist =  349.8127 

Iteration: 6 
Gradient vector:  56.4257 -33.54592 0.5844782 
New parameter estimates:  61382.92 94808.1 457.5005 

rse.dif =  690040.2 (rse = 1333584045 )  ;  parm.dist =  65.64701 

Iteration: 7 
Gradient vector:  -12.63262 7.648973 -0.1273216 
New parameter estimates:  61370.29 94815.75 457.3732 

rse.dif =  -140598.6 (rse = 1333443446 )  ;  parm.dist =  14.76841 

Iteration: 8 
Gradient vector:  2.7564 -1.660847 0.02796193 
New parameter estimates:  61373.05 94814.08 457.4012 

rse.dif =  31404.83 (rse = 1333474851 )  ;  parm.dist =  3.218219 

Iteration: 9 
Gradient vector:  -0.6050773 0.3649581 -0.006129619 
New parameter estimates:  61372.44 94814.45 457.395 

rse.dif =  -6859.57 (rse = 1333467992 )  ;  parm.dist =  0.7066474 

Iteration: 10 
Gradient vector:  0.1326532 -0.07999283 0.001344229 
New parameter estimates:  61372.57 94814.37 457.3964 

rse.dif =  1505.506 (rse = 1333469497 )  ;  parm.dist =  0.1549114 

Iteration: 11 
Gradient vector:  -0.02909032 0.01754299 -0.0002947644 
New parameter estimates:  61372.54 94814.39 457.3961 

rse.dif =  -330.0719 (rse = 1333469167 )  ;  parm.dist =  0.0339719 

Iteration: 12 
Gradient vector:  0.006378995 -0.003846825 6.463759e-05 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  72.38277 (rse = 1333469239 )  ;  parm.dist =  0.007449417 

Iteration: 13 
Gradient vector:  -0.001398820 0.0008435545 -1.417403e-05 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -15.8723 (rse = 1333469224 )  ;  parm.dist =  0.001633549 

Iteration: 14 
Gradient vector:  0.0003067399 -0.0001849785 3.108150e-06 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  3.480562 (rse = 1333469227 )  ;  parm.dist =  0.0003582123 

Iteration: 15 
Gradient vector:  -6.726341e-05 4.0563e-05 -6.815702e-07 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -0.7632339 (rse = 1333469226 )  ;  parm.dist =  7.855054e-05 

Iteration: 16 
Gradient vector:  1.474986e-05 -8.894858e-06 1.494581e-07 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  0.1673656 (rse = 1333469227 )  ;  parm.dist =  1.722496e-05 

Iteration: 17 
Gradient vector:  -3.234433e-06 1.950522e-06 -3.277396e-08 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -0.03670073 (rse = 1333469226 )  ;  parm.dist =  3.777191e-06 

Iteration: 18 
Gradient vector:  7.092709e-07 -4.277348e-07 7.186811e-09 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  0.008048058 (rse = 1333469226 )  ;  parm.dist =  8.282952e-07 

Iteration: 19 
Gradient vector:  -1.555217e-07 9.37892e-08 -1.575938e-09 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -0.001764536 (rse = 1333469226 )  ;  parm.dist =  1.816207e-07 

Iteration: 20 
Gradient vector:  3.410008e-08 -2.055502e-08 3.456296e-10 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  0.0003864765 (rse = 1333469226 )  ;  parm.dist =  3.982316e-08 

Iteration: 21 
Gradient vector:  -7.49142e-09 4.52347e-09 -7.593299e-11 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -8.416176e-05 (rse = 1333469226 )  ;  parm.dist =  8.755045e-09 

Iteration: 22 
Gradient vector:  1.650187e-09 -9.971753e-10 1.676294e-11 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  1.716614e-05 (rse = 1333469226 )  ;  parm.dist =  1.932973e-09 

Iteration: 23 
Gradient vector:  -3.678616e-10 2.341807e-10 -3.730106e-12 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -2.622604e-06 (rse = 1333469226 )  ;  parm.dist =  4.380867e-10 

Iteration: 24 
Gradient vector:  1.004182e-10 -7.590225e-11 8.078534e-13 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -1.430511e-06 (rse = 1333469226 )  ;  parm.dist =  1.251828e-10 

Iteration: 25 
Gradient vector:  -2.121985e-11 2.328697e-11 -2.071375e-14 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  2.861023e-06 (rse = 1333469226 )  ;  parm.dist =  3.637979e-11 

Iteration: 26 
Gradient vector:  -1.061002e-11 5.692926e-12 -6.3925e-14 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -9.536743e-07 (rse = 1333469226 )  ;  parm.dist =  7.27618e-12 

Convergence achieved by sums of squares.
Final parameter estimates:  61372.55 94814.38 457.3961 


ft.vrg> maas.vmod.fixed <- fit.variogram("gaussian", maas.v, 
    nugget = 60000, sill = 110000, range = 800, plot.it = TRUE, 
    iterations = 0)
Convergence not achieved!
> ## End Don't show
> # a single point:
> prdpnt <- point(data.frame(list(x=180000,y=331000)))
> prdpnt <- krige(prdpnt, maas.point, 'zinc', maas.vmod)

Using all points.
  Preparing the kriging system matrix...
  Inverting the matrix (and it's a big one)...
  Predicting.
> prdpnt

Point object: structure(list(x = 180000, y = 331000, do = TRUE, zhat = 179.126792748158,  
 
Point object:     sigma2hat = 72787.0311334009), .Names = c("x", "y", "do",  
 
Point object: "zhat", "sigma2hat"), row.names = "1", class = c("point", "data.frame" 
 
Point object: )) 

   Locations: 1

   Attributes:
      x
      y
      do
      zhat
      sigma2hat

> 
> # kriging on a grid (slow!)
> grid <- list(x=seq(min(maas$x),max(maas$x),by=100),
+              y=seq(min(maas$y),max(maas$y),by=100))
> grid$xr <- range(grid$x)
> grid$xs <- grid$xr[2] - grid$xr[1]
> grid$yr <- range(grid$y)
> grid$ys <- grid$yr[2] - grid$yr[1]
> grid$max <- max(grid$xs, grid$ys)
> grid$xy <- data.frame(cbind(c(matrix(grid$x, length(grid$x), length(grid$y))),
+              c(matrix(grid$y, length(grid$x), length(grid$y), byrow=TRUE))))
> colnames(grid$xy) <- c("x", "y")
> grid$point <- point(grid$xy)
> data(maas.bank)
> grid$krige <- krige(grid$point,maas.point,'zinc',maas.vmod,
+                     maxdist=1000,extrap=FALSE,border=maas.bank)

Using points within 1000 units of prediction points.
  Predicting....................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................
> op <- par(no.readonly = TRUE)
> par(pty="s")
> plot(grid$xy, type="n", xlim=c(grid$xr[1], grid$xr[1]+grid$max),
+                     ylim=c(grid$yr[1], grid$yr[1]+grid$max))
> image(grid$x,grid$y,
+       matrix(grid$krige$zhat,length(grid$x),length(grid$y)),
+       add=TRUE)
> contour(grid$x,grid$y,
+         matrix(grid$krige$zhat,length(grid$x),length(grid$y)),
+         add=TRUE)
> data(maas.bank)
> lines(maas.bank$x,maas.bank$y,col="blue")
> par(op)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "lagplot"
> 
> ### * lagplot
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lagplot
> ### Title: Lag Scatter Plot
> ### Aliases: lagplot
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> # prepare variables from other example pages if they are not already there:
> if(length(ls(pat="maas.pair"))==0){example(pair)}

pair> if (length(ls(pat = "maas.point")) == 0) {
    example(point)
}

point> data(maas)

point> maas.point <- point(maas)

pair> maas.pair <- pair(maas.point, num.lags = 10, maxdist = 2000)
..........................................................................................................................................................

pair> maas.pair25 <- pair(maas.point, num.lags = 10, type = "anisotropic", 
    theta = 25, maxdist = 500)
..............................
NOTE: Number of pairs in lag  1 :  0 
NOTE: Number of pairs in lag  2 :  1 
NOTE: Number of pairs in lag  3 :  3 
NOTE: Number of pairs in lag  4 :  12 
NOTE: Number of pairs in lag  5 :  10 
NOTE: Number of pairs in lag  6 :  12 
NOTE: Number of pairs in lag  7 :  10 
NOTE: Number of pairs in lag  8 :  13 
NOTE: Number of pairs in lag  9 :  10 
NOTE: Number of pairs in lag  10 :  19 

> ## End Don't show
> opar <- par(ask = interactive() && .Device == "X11")
> lagplot(maas.point,maas.pair,'zinc')
> # with identifying pairs:
> lagplot(maas.point,maas.pair,'zinc',lag=2,query.a='zinc')
Identify "from" points...
Identify "to" points...numeric(0)
> par(opar)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "pair"
> 
> ### * pair
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pair
> ### Title: Pair Object
> ### Aliases: pair
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> # prepare variables from other example pages if they are not already there:
> if(length(ls(pat="maas.point"))==0){example(point)}

point> data(maas)

point> maas.point <- point(maas)
> ## End Don't show
> maas.pair   <- pair(maas.point,num.lags=10,maxdist=2000)
..........................................................................................................................................................
> maas.pair25 <- pair(maas.point,num.lags=10,type='anisotropic',
+                       theta=25,maxdist=500)
..............................
NOTE: Number of pairs in lag  1 :  0 
NOTE: Number of pairs in lag  2 :  1 
NOTE: Number of pairs in lag  3 :  3 
NOTE: Number of pairs in lag  4 :  12 
NOTE: Number of pairs in lag  5 :  10 
NOTE: Number of pairs in lag  6 :  12 
NOTE: Number of pairs in lag  7 :  10 
NOTE: Number of pairs in lag  8 :  13 
NOTE: Number of pairs in lag  9 :  10 
NOTE: Number of pairs in lag  10 :  19 

> 
> 
> 
> cleanEx(); ..nameEx <- "plot.point"
> 
> ### * plot.point
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.point
> ### Title: Plot Point Objects
> ### Aliases: plot.point
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> if(length(ls(pat="maas.point"))==0){example(point)}

point> data(maas)

point> maas.point <- point(maas)
> ## End Don't show
> plot(maas.point) 
> plot(maas.point,v='zinc') 
> plot(maas.point,v='zinc',xlab='easting',ylab='northing',axes=TRUE,legend.pos=4)
> # plot additionally the maas bank:
> data(maas.bank)
> lines(maas.bank)
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.variogram"
> 
> ### * plot.variogram
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.variogram
> ### Title: Plot Variogram
> ### Aliases: plot.variogram
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> # prepare variables from other example pages if they are not already there:
> if(length(ls(pat="maas.vmod"))==0){example(fit.variogram)}

ft.vrg> if (length(ls(pat = "maas.v")) == 0) {
    example(est.variogram)
}

est.vr> if (length(ls(pat = "maas.pair")) == 0) {
    example(pair)
}

pair> if (length(ls(pat = "maas.point")) == 0) {
    example(point)
}

point> data(maas)

point> maas.point <- point(maas)

pair> maas.pair <- pair(maas.point, num.lags = 10, maxdist = 2000)
..........................................................................................................................................................

pair> maas.pair25 <- pair(maas.point, num.lags = 10, type = "anisotropic", 
    theta = 25, maxdist = 500)
..............................
NOTE: Number of pairs in lag  1 :  0 
NOTE: Number of pairs in lag  2 :  1 
NOTE: Number of pairs in lag  3 :  3 
NOTE: Number of pairs in lag  4 :  12 
NOTE: Number of pairs in lag  5 :  10 
NOTE: Number of pairs in lag  6 :  12 
NOTE: Number of pairs in lag  7 :  10 
NOTE: Number of pairs in lag  8 :  13 
NOTE: Number of pairs in lag  9 :  10 
NOTE: Number of pairs in lag  10 :  19 


est.vr> maas.v <- est.variogram(maas.point, maas.pair, "zinc")

ft.vrg> maas.vmod <- fit.gaussian(maas.v, c0 = 60000, cg = 110000, 
    ag = 800, plot.it = TRUE, iterations = 30)
Initial parameter estimates:  60000 110000 800 
Iteration: 1 
Gradient vector:  5499.075 -37894.01 -497.3869 
New parameter estimates:  65499.08 72106 302.6131 

rse.dif =  249665571 (rse = 249665571 )  ;  parm.dist =  38294.16 

Iteration: 2 
Gradient vector:  1877.678 14818.09 180.7598 
New parameter estimates:  67376.75 86924.08 483.373 

rse.dif =  1099461343 (rse = 1349126914 )  ;  parm.dist =  14937.67 

Iteration: 3 
Gradient vector:  -6669.642 8207.31 -37.26603 
New parameter estimates:  60707.11 95131.4 446.1069 

rse.dif =  -32082819 (rse = 1317044095 )  ;  parm.dist =  10575.7 

Iteration: 4 
Gradient vector:  909.8936 -484.5986 13.38276 
New parameter estimates:  61617 94646.8 459.4897 

rse.dif =  17615795 (rse = 1334659891 )  ;  parm.dist =  1030.981 

Iteration: 5 
Gradient vector:  -290.5093 194.8503 -2.573651 
New parameter estimates:  61326.5 94841.64 456.916 

rse.dif =  -1765886 (rse = 1332894005 )  ;  parm.dist =  349.8127 

Iteration: 6 
Gradient vector:  56.4257 -33.54592 0.5844782 
New parameter estimates:  61382.92 94808.1 457.5005 

rse.dif =  690040.2 (rse = 1333584045 )  ;  parm.dist =  65.64701 

Iteration: 7 
Gradient vector:  -12.63262 7.648973 -0.1273216 
New parameter estimates:  61370.29 94815.75 457.3732 

rse.dif =  -140598.6 (rse = 1333443446 )  ;  parm.dist =  14.76841 

Iteration: 8 
Gradient vector:  2.7564 -1.660847 0.02796193 
New parameter estimates:  61373.05 94814.08 457.4012 

rse.dif =  31404.83 (rse = 1333474851 )  ;  parm.dist =  3.218219 

Iteration: 9 
Gradient vector:  -0.6050773 0.3649581 -0.006129619 
New parameter estimates:  61372.44 94814.45 457.395 

rse.dif =  -6859.57 (rse = 1333467992 )  ;  parm.dist =  0.7066474 

Iteration: 10 
Gradient vector:  0.1326532 -0.07999283 0.001344229 
New parameter estimates:  61372.57 94814.37 457.3964 

rse.dif =  1505.506 (rse = 1333469497 )  ;  parm.dist =  0.1549114 

Iteration: 11 
Gradient vector:  -0.02909032 0.01754299 -0.0002947644 
New parameter estimates:  61372.54 94814.39 457.3961 

rse.dif =  -330.0719 (rse = 1333469167 )  ;  parm.dist =  0.0339719 

Iteration: 12 
Gradient vector:  0.006378995 -0.003846825 6.463759e-05 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  72.38277 (rse = 1333469239 )  ;  parm.dist =  0.007449417 

Iteration: 13 
Gradient vector:  -0.001398820 0.0008435545 -1.417403e-05 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -15.8723 (rse = 1333469224 )  ;  parm.dist =  0.001633549 

Iteration: 14 
Gradient vector:  0.0003067399 -0.0001849785 3.108150e-06 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  3.480562 (rse = 1333469227 )  ;  parm.dist =  0.0003582123 

Iteration: 15 
Gradient vector:  -6.726341e-05 4.0563e-05 -6.815702e-07 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -0.7632339 (rse = 1333469226 )  ;  parm.dist =  7.855054e-05 

Iteration: 16 
Gradient vector:  1.474986e-05 -8.894858e-06 1.494581e-07 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  0.1673656 (rse = 1333469227 )  ;  parm.dist =  1.722496e-05 

Iteration: 17 
Gradient vector:  -3.234433e-06 1.950522e-06 -3.277396e-08 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -0.03670073 (rse = 1333469226 )  ;  parm.dist =  3.777191e-06 

Iteration: 18 
Gradient vector:  7.092709e-07 -4.277348e-07 7.186811e-09 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  0.008048058 (rse = 1333469226 )  ;  parm.dist =  8.282952e-07 

Iteration: 19 
Gradient vector:  -1.555217e-07 9.37892e-08 -1.575938e-09 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -0.001764536 (rse = 1333469226 )  ;  parm.dist =  1.816207e-07 

Iteration: 20 
Gradient vector:  3.410008e-08 -2.055502e-08 3.456296e-10 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  0.0003864765 (rse = 1333469226 )  ;  parm.dist =  3.982316e-08 

Iteration: 21 
Gradient vector:  -7.49142e-09 4.52347e-09 -7.593299e-11 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -8.416176e-05 (rse = 1333469226 )  ;  parm.dist =  8.755045e-09 

Iteration: 22 
Gradient vector:  1.650187e-09 -9.971753e-10 1.676294e-11 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  1.716614e-05 (rse = 1333469226 )  ;  parm.dist =  1.932973e-09 

Iteration: 23 
Gradient vector:  -3.678616e-10 2.341807e-10 -3.730106e-12 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -2.622604e-06 (rse = 1333469226 )  ;  parm.dist =  4.380867e-10 

Iteration: 24 
Gradient vector:  1.004182e-10 -7.590225e-11 8.078534e-13 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -1.430511e-06 (rse = 1333469226 )  ;  parm.dist =  1.251828e-10 

Iteration: 25 
Gradient vector:  -2.121985e-11 2.328697e-11 -2.071375e-14 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  2.861023e-06 (rse = 1333469226 )  ;  parm.dist =  3.637979e-11 

Iteration: 26 
Gradient vector:  -1.061002e-11 5.692926e-12 -6.3925e-14 
New parameter estimates:  61372.55 94814.38 457.3961 

rse.dif =  -9.536743e-07 (rse = 1333469226 )  ;  parm.dist =  7.27618e-12 

Convergence achieved by sums of squares.
Final parameter estimates:  61372.55 94814.38 457.3961 


ft.vrg> maas.vmod.fixed <- fit.variogram("gaussian", maas.v, 
    nugget = 60000, sill = 110000, range = 800, plot.it = TRUE, 
    iterations = 0)
Convergence not achieved!
> ## End Don't show
> # two plots 
> oldpar <- par(mfrow=c(2,1))
> plot(maas.v) 
> plot(maas.v,var.mod.obj=maas.vmod)
> par(oldpar)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "point"
> 
> ### * point
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: point
> ### Title: Point Object
> ### Aliases: point
> ### Keywords: spatial
> 
> ### ** Examples
> 
> data(maas)
> maas.point <- point(maas) 
> 
> 
> 
> cleanEx(); ..nameEx <- "print.pair"
> 
> ### * print.pair
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: print.pair
> ### Title: Pairs Object Description
> ### Aliases: print.pair
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> # prepare variables from other example pages if they are not already there:
> if(length(ls(pat="maas.pair"))==0){example(pair)}

pair> if (length(ls(pat = "maas.point")) == 0) {
    example(point)
}

point> data(maas)

point> maas.point <- point(maas)

pair> maas.pair <- pair(maas.point, num.lags = 10, maxdist = 2000)
..........................................................................................................................................................

pair> maas.pair25 <- pair(maas.point, num.lags = 10, type = "anisotropic", 
    theta = 25, maxdist = 500)
..............................
NOTE: Number of pairs in lag  1 :  0 
NOTE: Number of pairs in lag  2 :  1 
NOTE: Number of pairs in lag  3 :  3 
NOTE: Number of pairs in lag  4 :  12 
NOTE: Number of pairs in lag  5 :  10 
NOTE: Number of pairs in lag  6 :  12 
NOTE: Number of pairs in lag  7 :  10 
NOTE: Number of pairs in lag  8 :  13 
NOTE: Number of pairs in lag  9 :  10 
NOTE: Number of pairs in lag  10 :  19 

> ## End Don't show
> print(maas.pair)

Pair object: maas.pair 

      Type:             isotropic
      Number of pairs:  8370
      Number of lags:   10
      Max h:            1999.867

> # gives:
> # Pairs object: maas.pair 
> #
> #      Type:             isotropic
> #      Number of pairs:  8370
> #      Number of lags:   10
> #      Max h:            1999.867
> 
> 
> 
> cleanEx(); ..nameEx <- "print.point"
> 
> ### * print.point
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: print.point
> ### Title: Point Object Description
> ### Aliases: print.point
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> if(length(ls(pat="maas.point"))==0){example(point)}

point> data(maas)

point> maas.point <- point(maas)
> ## End Don't show
> print.point(maas.point)

Point object: maas.point 

   Locations: 155

   Attributes:
      x
      y
      zinc

> # gives
> # Point object: maas.point 
> #
> #   Locations: 155
> #
> #   Attributes:
> #      x
> #      y
> #      zinc
> 
> 
> 
> cleanEx(); ..nameEx <- "spacebox"
> 
> ### * spacebox
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: spacebox
> ### Title: Boxplot of Variogram Cloud
> ### Aliases: spacebox
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> # prepare variables from other example pages if they are not already there:
> if(length(ls(pat="maas.pair"))==0){example(pair)}

pair> if (length(ls(pat = "maas.point")) == 0) {
    example(point)
}

point> data(maas)

point> maas.point <- point(maas)

pair> maas.pair <- pair(maas.point, num.lags = 10, maxdist = 2000)
..........................................................................................................................................................

pair> maas.pair25 <- pair(maas.point, num.lags = 10, type = "anisotropic", 
    theta = 25, maxdist = 500)
..............................
NOTE: Number of pairs in lag  1 :  0 
NOTE: Number of pairs in lag  2 :  1 
NOTE: Number of pairs in lag  3 :  3 
NOTE: Number of pairs in lag  4 :  12 
NOTE: Number of pairs in lag  5 :  10 
NOTE: Number of pairs in lag  6 :  12 
NOTE: Number of pairs in lag  7 :  10 
NOTE: Number of pairs in lag  8 :  13 
NOTE: Number of pairs in lag  9 :  10 
NOTE: Number of pairs in lag  10 :  19 

> ## End Don't show
> spacebox(maas.point,maas.pair,'zinc')
> 
> 
> 
> cleanEx(); ..nameEx <- "spacecloud"
> 
> ### * spacecloud
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: spacecloud
> ### Title: Variogram Cloud
> ### Aliases: spacecloud
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ## Don't show:
> # prepare variables from other example pages if they are not already there:
> if(length(ls(pat="maas.pair"))==0){example(pair)}

pair> if (length(ls(pat = "maas.point")) == 0) {
    example(point)
}

point> data(maas)

point> maas.point <- point(maas)

pair> maas.pair <- pair(maas.point, num.lags = 10, maxdist = 2000)
..........................................................................................................................................................

pair> maas.pair25 <- pair(maas.point, num.lags = 10, type = "anisotropic", 
    theta = 25, maxdist = 500)
..............................
NOTE: Number of pairs in lag  1 :  0 
NOTE: Number of pairs in lag  2 :  1 
NOTE: Number of pairs in lag  3 :  3 
NOTE: Number of pairs in lag  4 :  12 
NOTE: Number of pairs in lag  5 :  10 
NOTE: Number of pairs in lag  6 :  12 
NOTE: Number of pairs in lag  7 :  10 
NOTE: Number of pairs in lag  8 :  13 
NOTE: Number of pairs in lag  9 :  10 
NOTE: Number of pairs in lag  10 :  19 

> ## End Don't show
> opar <- par(ask = interactive() && .Device == "X11")
> spacecloud(maas.point,maas.pair,'zinc')
> # identify some points:
> spacecloud(maas.point,maas.pair,'zinc',query.a='zinc')
Identify "from" points...
Identify "to" points...
> par(opar)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  8.34 2.42 11.13 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')