R/est.variogram.R

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