# R/est.variogram.R In sgeostat: An Object-Oriented Framework for Geostatistical Modeling in S+

```#
assign("est.variogram",
function(point.obj,pair.obj,a1,a2) {

#  est.variogram takes a "point" object, point.obj, and a "pair" object, pair.obj,
#  calculates empirical variogram estimates.

#  The result is an object of type "variogram" with 4 components: \$lags,
#  \$classic, \$robust, and \$n.
#
#  \$lags - lag category
#  \$bins - distance bins for plotting
#  \$classic - the classic variogram estimator
#  \$robust - the robust variogram estimator
#  \$med - the median estimator
#  \$n - the number of pairs in the lag

if (!inherits(point.obj,"point")) stop('Point.obj must be of class, "point".\n')

if (!inherits(pair.obj,"pair")) stop('Pair.obj must be of class, "pair".\n')

if(missing(a1)) stop('Must enter at least one attribute.\n')
if(missing(a2)) a2 <- a1

a1 <- point.obj[[match(a1,names(point.obj))]]
a2 <- point.obj[[match(a2,names(point.obj))]]

# Allocate some space...
lags    <- sort(unique(pair.obj\$lags))
classic <- rep(0,length(lags))
robust  <- rep(0,length(lags))
med     <- rep(0,length(lags))
n       <- rep(0,length(lags))

diff <- a1[pair.obj\$from]-a2[pair.obj\$to]
bo   <- split(diff,pair.obj\$lags)

# this fails sometimes:
#  tmp<-unique(pair.obj\$lags[-which.na(pair.obj\$lags)])
# so do this:
if (any(is.na(pair.obj\$lags)))
tmp <- unique(pair.obj\$lags[-which.na(pair.obj\$lags)])
else
tmp <- unique(pair.obj\$lags)

for (i in c(1:length(tmp))) {
#  for (i in unique(pair.obj\$lags)) {
n[i] <- length(bo[[i]][!is.na(bo[[i]])])

#   classic, see Matheron
classic[i] <- sum((bo[[i]])^2,na.rm=TRUE) / n[i]

#   robust & med, see Cressie, 1990
robust[i] <- (sum(abs(bo[[i]])^.5,na.rm=TRUE) / n[i] )^4 / (0.457 + (0.494/n[i]))
med[i] <- (median(abs(bo[[i]])^.5,na.rm=TRUE))^4 / (0.457 + (0.494/n[i]))
#    med[i] <- median(abs(bo[[i]]),na.rm=TRUE)
}
o.variogram <- data.frame(lags,bins=c(pair.obj\$bins,recursive=TRUE),classic,
robust,med,n=n)
class(o.variogram) <- c("variogram","data.frame")
return(o.variogram)

})
```

## Try the sgeostat package in your browser

Any scripts or data that you put into this service are public.

sgeostat documentation built on May 29, 2017, 9:04 a.m.