# R/est.variograms.R In geospt: Geostatistical Analysis and Design of Optimal Spatial Sampling Networks

```assign("est.variograms",
function(point.obj,pair.obj,a1,a2,trim) {

#  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
#  \$trimmed.mean - the median trimmed estimator
#  \$n - the number of pairs in the lag
#  library(sgeostat, pos=which(search()=="package:gstat")+1)

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))
trimmed.mean <- 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 & trimmed mean, 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]))
trimmed.mean[i] <- (mean(abs(bo[[i]])^.5, trim=trim, 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,trimmed.mean,n=n)
class(o.variogram) <- c("variogram","data.frame")
return(o.variogram)

})
```

## Try the geospt package in your browser

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

geospt documentation built on May 29, 2017, 9:35 p.m.