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)

})