Variogram Estimator

Description

Calculate empirical variogram estimates.

An object of class variogram contains empirical variogram estimates generated from a point object and a pair object. A variogram object is stored as a data frame containing six columns: lags, bins, classic, robust, med, and n. The length of each vector is equal to the number of lags in the pair object used to create the variogram object, say l. The lags vector contains the lag numbers for each lag, beginning with one (1) and going to the number of lags (l). The bins vector contains the spatial midpoint of each lag. The classic, robust, and med vectors contain the classical,

\gamma_{c}(h)=\frac{1}{n}\sum_{(i,j) \in N(h)}(z(x_{i})-z(x_{j}))^{2}

robust,

\gamma_{m}(h)=\frac{(\frac{1}{n}\sum_{(i,j) \in N(h)} (\sqrt{|z(x_{i})-z(x_{j})|}))^{4}}{0.457 + \frac{0.494}{n}}

and median

\gamma_{m}(h)=\frac{(\mbox{median}_{(i,j) \in N(h)} (\sqrt{|z(x_{i})-z(x_{j})|}))^{4}}{0.457 + \frac{0.494}{|N(h)|}}

variogram estimates for each lag, respectively (see Cressie, 1993, p. 75). The n vector contains the number |N(h)| of pairs of points in each lag N(h).

Usage

est.variogram(point.obj, pair.obj, a1, a2) 

Arguments

Value

A variogram object:

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

point, pair

Examples

maas.v<-est.variogram(maas.point,maas.pair,'zinc') 

Fit polynomial trend functions

Description

Fits a polynomial trend function to a point object. Similar to functions in B. Ripleys spatial library.

Usage

fit.trend(point.obj, at, np=2, plot.it=TRUE)

Arguments

Value

...

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

Variogram Model Fit

Description

Fit variogram models (exponential, spherical, gaussian, linear) to empirical variogram estimates.

An object of class variogram.model represents a fitted variogram model generated by fitting a function to a variogram object. A variogram.model object is composed of a list consisting of a vector of parameters, parameters, and a semi-variogram model function, model.

Usage

fit.variogram(model="exponential", v.object, nugget, sill,
              range, slope, ...)
fit.exponential(v.object, c0, ce, ae, type='c', 
                iterations=10, tolerance=1e-06, echo=FALSE, plot.it=FALSE, weighted=TRUE)
fit.gaussian(v.object, c0, cg, ag, type='c', 
             iterations=10, tolerance=1e-06, echo=FALSE,  plot.it=FALSE, weighted=TRUE)
fit.spherical(v.object, c0, cs, as, type='c', iterations=10,
             tolerance=1e-06, echo=FALSE,  plot.it=FALSE, weighted=TRUE,
             delta=0.1, verbose=TRUE)
fit.wave(v.object, c0, cw, aw, type='c', 
         iterations=10, tolerance=1e-06, echo=FALSE,  plot.it=FALSE, weighted=TRUE)
fit.linear(v.object, type='c', plot.it=FALSE,iterations=1, c0=0, cl=1)

Arguments

Value

A variogram.model object:

Note

fit.exponential, fit.gaussian and fit.wave use an iterative, Gauss-Newton fitting algorithm to fit to an exponential or gaussian variogram model to empirical variogram estimates. fit.spherical uses the same algorithm but with differential quotients in place of first derivatives. When weighted is TRUE, the regression is weighted by n(h)/gamma(h)^2 where the numerator is the number of pairs of points in a given lag.

Setting iterations to 0 means no fit procedure is applied. Thus parameter values from external sources can be plugged into a variogram model object.

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

est.variogram

Examples

#
# automatic fit:
#
maas.vmod<-fit.gaussian(maas.v,c0=60000,cg=110000,ag=800,plot.it=TRUE,
                        iterations=30)
#
# 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)

Identify points on a Point Object

Description

Plot variable values next to locations after the plot.point() function.

Usage

## S3 method for class 'point'
identify(x, v, ...) 

Arguments

Value

An integer vector containing the indexes of the identified points.

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

plot.point

Examples

plot(maas.point)
# use indices as labels:
identify(maas.point)
# use values as labels:
identify(maas.point,v="zinc")

Convex hull test

Description

Checks if points are in the interior of a convex hull.

Usage

in.chull(x0, y0, x, y)

Arguments

Details

Uses a simple points-in-polygon check combined with the chull function.

Value

Author(s)

Albrecht Gebhardt <agebhard@uni-klu.ac.at>

References

Follows an idea from algorithm 112 from CACM (available at http://www.netlib.org/tomspdf/112.pdf)

See Also

in.convex.hull, chull

Examples

in.chull(c(0,1),c(0,1),c(0,1,0,-1),c(-1,0,1,0))
# should give: TRUE FALSE

In-Polygon test

Description

Checks if points are in the interior of a polygon.

Usage

in.polygon(x0, y0, x, y)

Arguments

Details

Uses a simple points-in-polygon check combined with the polygon function.

Polygon is closed automatically.

Value

Author(s)

Albrecht Gebhardt <agebhard@uni-klu.ac.at>

References

Follows an idea from algorithm 112 from CACM (available at http://www.netlib.org/tomspdf/112.pdf)

See Also

in.convex.hull, polygon, in.chull

Examples

in.polygon(c(0,1),c(0,1),c(0,1,0,-1),c(-1,0,1,0))
# should give: TRUE FALSE

Kriging

Description

Carry out spatial prediction (or kriging).

Usage

krige(s, point.obj, at, var.mod.obj, maxdist=NULL, extrap=FALSE, border) 

Arguments

Value

A point object which is a copy of the s object with two new variables, zhat and sigma2hat, which are, repspectively, the predicted value and the kriging variance.

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

est.variogram,fit.variogram

Examples

# a single point:
prdpnt <- point(data.frame(list(x=180000,y=331000)))
prdpnt <- krige(prdpnt, maas.point, 'zinc', maas.vmod)
prdpnt

# 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)
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)

Lag Scatter Plot

Description

Create a spatially lagged scatter plot, e.g. plot z(s) versus z(s+h), where h is a lag in a pair object.

Usage

lagplot(point.obj, pair.obj, a1, a2, lag=1, std=FALSE, query.a, xlim, ylim)

Arguments

Value

NULL

Note

When query.a is entered, the user will be prompted to identify points on the display device. Because each point in the plot represents a pair of locations, the user must identify each point twice, once for the "from" point and once for the "to" point. Querying is ended by pressing the middle mouse button on the mouse while the cursor is in the display window.

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

point, pair

Examples

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')
par(opar)

maas- zinc measurements

Description

Zinc measurements as groundwater quality variable.

Usage

maas

Value

list with components x,y and zinc.

References

gstat E.J Pebesma (E.J.Pebesma@frw.uva.nl) http://www.frw.uva.nl/~pebesma/gstat/

See Also

maas.bank

maas.bank - coordinates

Description

Coordinates of maas bank. To be used together with maas.

Usage

maas.bank

Value

list with components x andy..

References

gstat E.J Pebesma (E.J.Pebesma@frw.uva.nl) http://www.frw.uva.nl/~pebesma/gstat/

See Also

maas

Pair Object

Description

Create a pair object from a point object.

A pair object contains information defining pairs of points contained in a point object. A pair object is a list containing five vectors: from, to, lags, dist, and bins. The length of each of these vectors (except bins) is equal to the number of pairs of points being represented, say k. The vectors from and to contain pointers into the vectors of a point object, pointing to each member of the pair of points (e.g., from[k] points to si and to[k] points to sj). The vector dist contains the distance between the pairs of points. The vector lags contains the lag number to which each pair of points has been assigned. The vector bins contains the spatial midpoint between each lag and is used for plotting.

Usage

pair(point.obj,num.lags=10,type='isotropic', theta=0, dtheta=5, maxdist) 

Arguments

Value

A pair object:

NOTE

Name of this function changed from pairs to pair to avoid conflicts with R's pairs function!!

Note

When creating an anisotropic pair object, the assumption is that the direction, as well as the distance, between pairs of points is important in describing the variation. Using the theta and dtheta arguments, pairs of points that meet direction requirements can be selected. A pair of points will be included when the angle between the positive x axis and the vector formed by the pair of points falls within the tolerance angle given by (theta-dtheta,theta+dtheta)

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

point

Examples

maas.pair   <- pair(maas.point,num.lags=10,maxdist=2000)
maas.pair25 <- pair(maas.point,num.lags=10,type='anisotropic',
                      theta=25,maxdist=500)

Plot Point Objects

Description

Plot the spatial locations in a point object, optionally coloring by quantile.

Usage

## S3 method for class 'point'
plot(x, v, legend.pos=0,axes=TRUE,xlab='',ylab='', add=FALSE, ...) 

Arguments

Value

NULL

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

point

Examples

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)

Plot Variogram

Description

Plot empirical variogram estimates, optionally plotting a fitted variogram model.

Usage

## S3 method for class 'variogram'
plot(x, var.mod.obj, title.str,ylim, type='c',N=FALSE, ...) 

Arguments

Value

NULL

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

est.variogram

Examples

# two plots 
oldpar <- par(mfrow=c(2,1))
plot(maas.v) 
plot(maas.v,var.mod.obj=maas.vmod)
par(oldpar)

Point Object

Description

Create an object of class point from a data frame.

An object of class point represents the observed data of a spatial process. This includes the spatial location of sampling sites and the values observed at those sites. A point object is stored as a data frame. The data frame must contain one column for the X coordinate and one column for the Y coordinate of each point, as well as any number of columns representing data observed at the points.

Usage

point(dframe, x='x', y='y') 

Arguments

Value

A point object:

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

point

Examples

data(maas)
maas.point <- point(maas) 

Pairs Object Description

Description

Print descriptive information about a pair object.

Usage

## S3 method for class 'pair'
print(x,...) 

Arguments

Value

NULL

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

pair

Examples

print(maas.pair)
# gives:
# Pairs object: maas.pair 
#
#      Type:             isotropic
#      Number of pairs:  8370
#      Number of lags:   10
#      Max h:            1999.867

Point Object Description

Description

Print descriptive information about a point object.

Usage

## S3 method for class 'point'
print(x,...) 

Arguments

Value

NULL

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

point

Examples

print.point(maas.point)
# gives
# Point object: maas.point 
#
#   Locations: 155
#
#   Attributes:
#      x
#      y
#      zinc

Internal sgeostat functions

Description

Internal sgeostat functions

Details

These functions are not intended to be called by the user.

The krige function interfaces to krige.*, pair to pair.* and fit.trend to trend.*.

Boxplot of Variogram Cloud

Description

Create boxplots of square-root or squared differences between variable values at pairs of points versus the distance between the points.

Usage

spacebox(point.obj, pair.obj, a1, a2, type='r')

Arguments

Value

NULL

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

point, pair

Examples

spacebox(maas.point,maas.pair,'zinc')

Variogram Cloud

Description

Create a scatter plot of square-root or squared differences between variable values at pairs of points versus the distance between the points.

Usage

spacecloud(point.obj, pair.obj, a1, a2, type='r', query.a, ...)

Arguments

Value

NULL

Note

When query.a is entered, the user will be prompted to identify points on the display device. Because each point in the plot represents a pair of locations, the user must identify each point twice, once for the "from" point and once for the "to" point. Querying is ended by pressing the middle mouse button on the mouse while the cursor is in the display window.

References

http://www.gis.iastate.edu/SGeoStat/homepage.html

See Also

point, pair

Examples

opar <- par(ask = interactive() && .Device == "X11")
spacecloud(maas.point,maas.pair,'zinc')
# identify some points:
spacecloud(maas.point,maas.pair,'zinc',query.a='zinc')
par(opar)