R/variogramST.R

In gstat: Spatial and Spatio-Temporal Geostatistical Modelling, Prediction and Simulation

VgmFillNA <- function(x, boundaries) {
  # pads the sample variogram with NA rows where no data are available.
	n = length(boundaries) - 1
	ix = rep(NA, n)
	#ix[which(1:n %in% findInterval(x$dist, boundaries))] = 1:nrow(x)
	ix[ findInterval(x$dist, boundaries)] = 1:nrow(x)
	# x$b = boundaries[-1]
	x[ix,]
}

VgmAverage = function(ret, boundaries) {
	# take out NULL variograms:
	ret = ret[!sapply(ret, is.null)]
	# take care of missing rows...
	ret = lapply(ret, VgmFillNA, 
			boundaries = c(0, 1e-6 * boundaries[2], boundaries[-1]))
	# weighted average:
  # sum np, weighted sum of gamma and dist devided by summed np
	np = apply(do.call(cbind, lapply(ret, function(x) x$np)), 1, sum,
	           na.rm = TRUE)
	gamma = apply(do.call(cbind, lapply(ret, function(x) x$gamma*x$np)), 1, sum,
		na.rm = TRUE)/np
	dist = apply(do.call(cbind, lapply(ret, function(x) x$dist*x$np)), 1, sum,
		na.rm = TRUE)/np
	v = data.frame(np = np, dist = dist, gamma = gamma)
	class(v) = class(ret[[1]])
	attr(v, "boundaries") = attr(ret[[1]], "boundaries")
	v[is.na(v)] = NA
	v
}

StVgmLag = function(formula, data, dt, pseudo, ...) {
  dotLst <- list(...)
	.ValidObs = function(formula, data)
		!is.na(data[[as.character(as.list(formula)[[2]])]])
	d = dim(data)
	if (formula[[3]] != 1) { # there is a regression model:
		data$resid = residuals(lm(formula, data, na.action = na.exclude))
		formula = resid ~ 1
	}
	ret = vector("list", d[2] - dt)
	if (dt == 0) {
		for (i in 1:d[2]) {
			d0 = data[,i]
			valid = .ValidObs(formula, d0)
			if(sum(valid) <= 1)
				ret[[i]] <- NULL
			else {
				d0 = d0[valid,]
				ret[[i]] = variogram(formula, d0, ...)
			}
		}
	} else {
		for (i in 1:(d[2] - dt)) {
			d1 = data[, i]
			valid1 = .ValidObs(formula, d1)
			d2 = data[, i + dt]
			valid2 = .ValidObs(formula, d2)
			if(sum(valid1)==0 || sum(valid2)==0)
				ret[[i]] <- NULL
			else {
				d1 = d1[valid1,]
				d2 = d2[valid2,]
				obj = gstat(NULL, paste("D", i, sep=""), formula, d1, 
					set = list(zero_dist = 3), beta = 0)
				obj = gstat(obj, paste("D", i+dt, sep=""), formula, d2, 
					beta = 0)
				ret[[i]] = variogram(obj, cross = "ONLY", pseudo = pseudo, ...)
			}
		}
	}
  if(!is.null(dotLst$cloud)) {
    if(dotLst$cloud)
      ret <- do.call("rbind", ret)
      ret$id <- "var1"
      return(ret)
  } else {
	   return(VgmAverage(ret, dotLst$boundaries))
  }
}

variogramST = function(formula, locations, data, ..., tlags = 0:15, cutoff, 
                       width = cutoff/15, boundaries=seq(0,cutoff,width),
                       progress = interactive(), pseudo = TRUE, 
                       assumeRegular=FALSE, na.omit=FALSE) {
	if (missing(data))
		data = locations
	if(missing(cutoff)) {
		ll = !is.na(is.projected(data@sp)) && !is.projected(data@sp)
		cutoff <- spDists(t(data@sp@bbox), longlat = ll)[1,2]/3
	}
  
	if(is(data, "STIDF"))
		return(variogramST.STIDF(formula, data, tlags, cutoff, width, 
                             boundaries, progress, ...))
  
	stopifnot(is(data, "STFDF") || is(data, "STSDF"))
	it = index(data@time)
	if (assumeRegular || is.regular(zoo(matrix(1:length(it)), order.by = it),
                                  strict = TRUE)) {
		twidth = diff(it)[1]
		tlags = tlags[tlags <= min(max(tlags), length(unique(it)) - 1)]
	} else {
		warning("strictly irregular time steps were assumed to be regular")
		twidth = mean(diff(it))
	}
	ret = vector("list", length(tlags))
	obj = NULL
	t = twidth * tlags
	if (progress)
		pb = txtProgressBar(style = 3, max = length(tlags))
	for (dt in seq(along = tlags)) {
		ret[[dt]] = StVgmLag(formula, data, tlags[dt], pseudo = pseudo, 
                         boundaries = boundaries, ...)
		ret[[dt]]$id = paste("lag", dt - 1, sep="")
		if (progress)
			setTxtProgressBar(pb, dt)
	}
	if (progress)
		close(pb)
	# add time lag:
	v = do.call(rbind, ret)
	v$timelag = rep(t, sapply(ret, nrow))
	if (is(t, "yearmon"))
		class(v$timelag) = "yearmon"
    
	b = attr(ret[[min(length(tlags),2)]], "boundaries")
	b = c(0, b[2]/1e6, b[-1])
	# ix = findInterval(v$dist, b) will use all spacelags
	b = b[-2]
	# spacelags = c(0, b[-length(b)] + diff(b)/2) will use all spacelags
	v$spacelag = c(0, b[-length(b)] + diff(b)/2) # spacelags[ix] will use all spacelags

	v$avgDist <- v$dist * v$np
	for (lagId in unique(v$spacelag)) {
	  bool <- v$spacelag == lagId
	  v$avgDist[bool] <- sum(v$avgDist[bool], na.rm = TRUE) / sum(v$np[bool], na.rm = TRUE)
	}
	

  class(v) = c("StVariogram", "data.frame")
	if(na.omit)
    v <- na.omit(v)

  # setting attributes to allow krigeST to double check metrics
	attr(v$timelag,"units") <- attr(twidth,"units")
	if (isTRUE(!is.projected(data)))
		attr(v$spacelag, "units") = "km"
  
  return(v)
}

## very irregular data
variogramST.STIDF <- function (formula, data, tlags, cutoff, 
                               width, boundaries, progress, 
                               twindow, tunit) {
  ll = !is.na(is.projected(data@sp)) && !is.projected(data@sp)
  
  if (missing(cutoff))
    cutoff <- spDists(t(data@sp@bbox), longlat = ll)[1, 2]/3
  
  m = model.frame(terms(formula), as(data, "data.frame"))
  
  diffTime <- diff(index(data))
  timeScale <- units(diffTime)
  if(missing(tunit))
    warning(paste("The argument 'tunit' is missing: tlags are assumed to be given in ", timeScale, ".",sep=""))
  else {
    stopifnot(tunit %in% c("secs", "mins", "hours", "days", "weeks"))
    units(diffTime) <- tunit
    timeScale <- tunit
  }
  diffTime <- as.numeric(diffTime)
  if (missing(twindow)) {
    twindow <- round(2 * max(tlags, na.rm=TRUE)/mean(diffTime,na.rm=TRUE),0)
  }
    
  nData <- nrow(data)
  
  # re-using the order propertie of the time slot to only store the next "twindow" distances
  numTime <- as.numeric(index(data))
  diffTimeMat <- matrix(NA, nData, twindow)
  
  for (i in 1:nData) { # i <- 1
    diffTimeMat[i,1:min(nData,twindow)] <- cumsum(diffTime[i+0:(min(nData,twindow)-1)])
  }
  
  nSp <- length(boundaries)
  nTp <- length(tlags)
  
  distTp <- matrix(NA, nSp-1, nTp-1)
  distSp <- matrix(NA, nSp-1, nTp-1)
  gamma  <- matrix(NA, nSp-1, nTp-1)
  np     <- matrix(NA, nSp-1 ,nTp-1)
  
  # temporal selection
  if(progress)
    pb <- txtProgressBar(0, nTp-1, 0, style=3)
  for (i in 1:(nTp-1)) { # i <- 1
    ind <- which(diffTimeMat >= tlags[i] & diffTimeMat < tlags[i+1])
    if (length(ind) < 1) 
      next
    tmpInd <- matrix(NA,nrow=length(ind),4)
    tmpInd[,1] <- ind %% nData  # row number
    tmpInd[,2] <- (ind %/% nData)+1 # col number
    tmpInd[,3] <- apply(tmpInd[,1:2,drop=FALSE], 1, function(x) spDists(data@sp[x[1]], data@sp[x[2]+x[1],]))
    tmpInd[,4] <- diffTimeMat[tmpInd[,1:2, drop=FALSE]]
    
    # spatial selection
    for (j in 1:(nSp-1)) { # j <- 3
      indSp <- which(tmpInd[,3] >= boundaries[j] & tmpInd[,3] < boundaries[j+1])
      if (length(indSp) < 1)
        next
      distSp[j,i] <- mean(tmpInd[indSp,3])
      distTp[j,i] <- mean(tmpInd[indSp,4])
      
      indSp <- cbind(ind[indSp] %% nData, (ind[indSp] %/% nData)+1)
      np[j,i] <- length(indSp)
      gamma[j,i] <- 0.5*mean((data[indSp[,1],,colnames(m)[1]]@data[[1]] - data[indSp[,1]+indSp[,2],,colnames(m)[1]]@data[[1]])^2,
                             na.rm=TRUE)
    }
    if(progress)
      setTxtProgressBar(pb, value=i)
  }
  if(progress)
    close(pb)
  
  res <- data.frame(np=as.vector(np),
                    dist=as.vector(distSp),
                    gamma=as.vector(gamma),
                    id=paste("lag",rep(1:(nTp-1),each=nSp-1), sep=""),
                    timelag=rep(tlags[-nTp]+diff(tlags)/2,each=nSp-1),
                    spacelag=rep(boundaries[-nSp]+diff(boundaries)/2, nTp-1))
  
  res$avgDist <- res$dist * res$np
  for (lagId in unique(res$spacelag)) {
    bool <- res$spacelag == lagId
    res$avgDist[bool] <- sum(res$avgDist[bool], na.rm = TRUE) / sum(res$np[bool], na.rm = TRUE)
  }
  
  attr(res$timelag, "units") <- timeScale
  attr(res$spacelag, "units") <- ifelse(ll, "km", "m")
  class(res) <- c("StVariogram", "data.frame")
  
  return(res)
}



## plotting
plot.StVariogram = function(x, model=NULL, ..., col = bpy.colors(), xlab, ylab,
                            map = TRUE, convertMonths = FALSE, as.table = TRUE,
                            wireframe = FALSE, diff = FALSE, all=FALSE) {
	lst = list(...)
	if (!is.null(lst$col.regions))
		col = lst$col.regions
	if (is(x$timelag, "yearmon")) {
		if (convertMonths) {
			x$timelag = as.numeric(x$timelag) * 12
			attr(x$timelag, "units") = "months"
		} else
			attr(x$timelag, "units") = "years"
	}
	if (missing(xlab)) {
		xlab = "distance"
		u =  attr(x$spacelag, "units")
		if (!is.null(u))
			xlab = paste(xlab, " (", u, ")", sep="")
	}
	if (missing(ylab)) {
		ylab = "time lag"
		u = attr(x$timelag, "units")
		if (!is.null(u))
			ylab = paste(ylab, " (", u, ")", sep="")
	}
  
  # check for older spatio-temporal variograms and compute avgDist on demand
  if(is.null(x$avgDist)) {
    x$avgDist <- x$dist * x$np
    for (lagId in unique(x$spacelag)) {
      bool <- x$spacelag == lagId
      x$avgDist[bool] <- sum(x$avgDist[bool] / sum(x$np[bool], na.rm = TRUE), na.rm = TRUE)
    }
    
  }
  
  if(!is.null(model)) {
    if (is(model,"StVariogramModel"))
      model <- list(model)
    for (mod in model) {
      slag <- x$avgDist
      slag[slag == 0 & x$timelag == 0] <- sqrt(.Machine$double.eps)
      x[[mod$stModel]] <- variogramSurface(mod, data.frame(spacelag = slag,
                                                           timelag = x$timelag))$model
      if (diff)
        x[[mod$stModel]] <- x[[mod$stModel]] - x$gamma
    }
  }
	x0 = x # needed by wireframe()
	
	if (!is.null(model)) {
    modelNames  <- sapply(model, function(x) x$stModel)
    
    if(all &! diff)
      v0 <- x[,c("dist", "id", "avgDist", "timelag")]
    else
      v0 <- NULL
    
    for (i in modelNames)
      v0 <- rbind(v0, x[,c("dist", "id", "avgDist", "timelag")])

    if(all & !diff) {# we also need the sample
      v0$what = factor(c(rep("sample", nrow(x)), rep(modelNames, each=nrow(x))),
                       levels=c("sample", modelNames), ordered = TRUE)
      v0$gamma = c(x$gamma, unlist(x[,modelNames]))
    } else {
      v0$what = factor(rep(modelNames, each=nrow(x)),
                       levels=modelNames, ordered=TRUE)
      v0$gamma = c(unlist(x[,modelNames]))
    }
		
		x = v0
	}
	if (wireframe) { 
		if (!is.null(model)) {
      if (length(model) > 1)
        wireframe(gamma ~ avgDist*timelag | what, 
                  x, drape = TRUE, col.regions = col, 
                  xlab = xlab, ylab = ylab, as.table=as.table, ...)
      else 
        wireframe(as.formula(paste(model[[1]]$stModel,"~ avgDist*timelag")), 
                  x0, drape = TRUE, col.regions = col, 
                  xlab = xlab, ylab = ylab, as.table=as.table, ...)
		} else # without a model, plot only the sample variogram as a wireframe
			wireframe(gamma ~ avgDist * timelag, x0, drape = TRUE, col.regions = col,
				xlab = xlab, ylab = ylab, ...)
	} else if (map) {
		if (!is.null(model))
			f = gamma ~ avgDist + timelag | what
		else
			f = gamma ~ avgDist + timelag
		levelplot(f, x, xlab = xlab, ylab = ylab, col.regions = col, as.table=as.table, ...)
	} else { # not map, not wireplot
		if (!is.null(model))
			f = gamma ~ dist | what
		else
			f = gamma ~ dist
		x$id = factor(x$id, levels=unique(x$id))
		bp = bpy.colors(length(levels(x$id)))
		ps = list(superpose.line=list(col=bp), superpose.symbol=list(col=bp))
		xlim = c(0, max(x$dist) * 1.04)
		if (diff) {
		  xyplot(f, x, groups = x$id, type='b', xlim = xlim,
		         auto.key = list(space = "right"), xlab = xlab, 
		         par.settings = ps, as.table=as.table, ...)
		} else {
  		ylim = c(0, max(x$gamma) * 1.04)
  		xyplot(f, x, groups = x$id, type='b', ylim = ylim, xlim = xlim,
  		       auto.key = list(space = "right"), xlab = xlab, 
  		       par.settings = ps, as.table=as.table, ...)
    }
	}
}

print.StVariogramModel <- function(x, ...) {
  possComp <- c("space", "time", "joint")
  for(comp in possComp[possComp %in% names(x)]) {
    rownames(x[[comp]]) <- 1:nrow(x[[comp]])
    cat(paste(comp,"component: \n"))
    print(x[[comp]], ...)
  }
  
  possAddPar <- c("sill", "nugget", "stAni", "k")
  for(addPar in possAddPar[possAddPar %in% names(x)]) {
    cat(paste(addPar, ": ",x[[addPar]],"\n", sep=""))
  }
}