R/stVariogramModels.R

Defines functions insertParMetric insertParSimpleSumMetric insertParSumMetric insertParProdSum insertParProdSumOld insertParSeparable extractParNames extractPar insertPar fit.StVariogram covSurfMetric covMetric vgmMetric covSurfSimpleSumMetric covSimpleSumMetric vgmSimpleSumMetric covSurfSumMetric covSumMetric vgmSumMetric covSurfProdSum covProdSum vgmProdSum covSurfProdSumOld covProdSumOld vgmProdSumOld covSurfSeparable covSeparable vgmSeparable variogramSurface vgmST

Documented in extractPar extractParNames fit.StVariogram variogramSurface vgmST

# constructiong spatio-temporal variogram models
vgmST <- function(stModel, ..., space, time, joint, sill, k, nugget, stAni, 
                  temporalUnit) {
  stopifnot(is.character(stModel) && length(stModel)==1)
  
  old.stModel <- stModel
  stModel <- strsplit(stModel, "_")[[1]][1]
  
  if (stModel == "productSum" & !missing(sill))
    stop("The sill argument for the product-sum model has been removed 
due a change in notation of the spatio-temporal models. This 
affects as well how the spatial and temporal variograms are parameterised. 
Re-fit your model or use \"productSumOld\" instead.")
  
  if(!missing(sill))
    if(sill <= 0) stop("\"sill\" must be positive.")
  if(!missing(k))
    if(k <= 0) stop("\"k\" must be positive.")
  if(!missing(nugget))
    if(nugget < 0) stop("\"nugget\" must be non-negative.")
  if(!missing(stAni))
    if(stAni <= 0) stop("\"stAni\" must be positive.")
  
  vgmModel <- switch(stModel,
                     separable = list(space = space, time = time, sill = sill),
                     productSum = list(space = space, time = time, k = k),
                     productSumOld = list(space = space, time = time,
                                          sill = sill, nugget = nugget),
                     sumMetric = list(space = space, time = time, 
                                      joint = joint, stAni = stAni),
                     simpleSumMetric = list(space = space, time = time, 
                                            joint = joint, nugget = nugget, 
                                            stAni = stAni),
                     metric = list(joint = joint, stAni = stAni),
                     stop(paste("model", stModel, "unknown")))
  
  vgmModel$stModel <- old.stModel
  
  if (!missing(temporalUnit))
    attr(vgmModel, "temporal unit") = temporalUnit
  
  class(vgmModel) <- c("StVariogramModel", "list")
  
  vgmModel
}

# calculating spatio-temporal variogram surfaces
variogramSurface <- function(model, dist_grid, covariance=FALSE) {
  stopifnot(inherits(model, "StVariogramModel"))
  stopifnot(all(c("spacelag", "timelag") %in% colnames(dist_grid)))
  
  if (covariance) {
    switch(strsplit(model$stModel, "_")[[1]][1],
           separable=covSurfSeparable(model, dist_grid),
           productSum=covSurfProdSum(model, dist_grid),
           productSumOld=covSurfProdSumOld(model, dist_grid),
           sumMetric=covSurfSumMetric(model, dist_grid),
           simpleSumMetric=covSurfSimpleSumMetric(model, dist_grid),
           metric=covSurfMetric(model, dist_grid),
           stop("Only \"separable\", \"productSum\", \"sumMetric\", \"simpleSumMetric\" and \"metric\" are implemented."))
  } else {
    switch(strsplit(model$stModel, "_")[[1]][1],
           separable=vgmSeparable(model, dist_grid),
           productSum=vgmProdSum(model, dist_grid),
           productSumOld=vgmProdSumOld(model, dist_grid),
           sumMetric=vgmSumMetric(model, dist_grid),
           simpleSumMetric=vgmSimpleSumMetric(model, dist_grid),
           metric=vgmMetric(model, dist_grid),
           stop("Only \"separable\", \"productSum\", \"sumMetric\", \"simpleSumMetric\" and \"metric\" are implemented."))
  }
}

################################
## separable model: C_s * C_t ##
################################

vgmSeparable <- function(model, dist_grid) {
  vs = variogramLine(model$space, dist_vector=dist_grid$spacelag)$gamma
  vt = variogramLine(model$time,  dist_vector=dist_grid$timelag)$gamma
  
  cbind(dist_grid, "gamma" = model$sill*(vs+vt-vs*vt))
}

covSeparable <- function(x, y, model, separate) {  
  if(missing(separate))
    separate <- inherits(x, "STF") & inherits(y, "STF") & length(x) > 1 & length(y) > 1
  
  # the STF case
  if (inherits(x, "STF") && inherits(y, "STF")) {
    # calculate all spatial and temporal distances
    ds = spDists(x@sp, y@sp)
    dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    debug_time_unit(units(dt))
    dt <- as(dt, "matrix")
    
    # compose the cov-matrix
    Sm = variogramLine(model$space, covariance = TRUE, dist_vector = ds)*model$sill
    Tm = variogramLine(model$time, covariance = TRUE, dist_vector = dt)
    
    if (separate)
      return(list(Sm = Sm, Tm = Tm))
    else
      return(Tm %x% Sm) # kronecker product
  } 
  
  # separate makes only sense if both of x and y inherit STF
  if (separate)
    stop("An efficient inversion by separating the covarinace model is only possible if both of \"x\" and \"y\" inherit \"STF\"")
  
  # the STI case
  if (inherits(x, "STI") || inherits(y, "STI")) {
    # make sure that now both are of type STI
    x <- as(x, "STI")
    y <- as(y, "STI")
    
    # calculate all spatial and temporal distances
    ds = spDists(x@sp, y@sp)
    dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    debug_time_unit(units(dt))
    dt <- as(dt, "matrix")
    
    # compose the cov-matrix
    Sm = variogramLine(model$space, covariance = TRUE, dist_vector = ds)*model$sill
    Tm = variogramLine(model$time, covariance = TRUE, dist_vector = dt)
    
    return(Sm * Tm)
  }
  
  # the remaining cases, none of x and y is STI nor are both STF
  # make sure both are of type STS
  x <- as(x, "STS")
  y <- as(y, "STS")
  
  # calculate all spatial and temporal distances
  ds = spDists(x@sp, y@sp)
  dt = abs(outer(index(x@time), index(y@time), "-"))
  if(!is.null(attr(model,"temporal unit")))
    units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
  debug_time_unit(units(dt))
  dt <- as(dt, "matrix")
  
  # re-arrange the spatial and temporal distances
  sMat <- matrix(NA, nrow(x@index), nrow(y@index))
  tMat <- matrix(NA, nrow(x@index), nrow(y@index))
  for(r in 1:nrow(x@index)) {
    sMat[r,] <- ds[x@index[r,1], y@index[,1]]
    tMat[r,] <- dt[x@index[r,2], y@index[,2]]
  }
  
  # compose the cov-matrix
  Sm = variogramLine(model$space, covariance = TRUE, dist_vector = sMat)*model$sill
  Tm = variogramLine(model$time, covariance = TRUE, dist_vector = tMat)
  
  return(Sm * Tm)  
}

# covariance for the circulant embedding in ST
covSurfSeparable <- function(model, dist_grid) {
  Sm = variogramLine(model$space, covariance = TRUE, dist_vector = dist_grid$spacelag)$gamma*model$sill
  Tm = variogramLine(model$time, covariance = TRUE, dist_vector = dist_grid$timelag)$gamma
  
  cbind(dist_grid, "gamma" = Tm * Sm)
}

###########################################
## productSum model: C_s*C_t + C_s + C_t ##
###########################################

vgmProdSumOld <- function(model, dist_grid) {
  .Deprecated("vgmProdSum", package = "gstat", 
              msg="The former product-sum model is dprecited, consider to refit the new model specification",
              old = "vgmProdSumOld")
  
  vs = variogramLine(model$space, dist_vector=dist_grid$spacelag)$gamma
  vt = variogramLine(model$time, dist_vector=dist_grid$timelag)$gamma
  vn <- rep(model$nugget, length(vs))
  vn[vs == 0 & vt == 0] <- 0
  
  k <- (sum(model$space$psill)+sum(model$time$psill)-(model$sill+model$nugget))/(sum(model$space$psill)*sum(model$time$psill))
  
  if (k <= 0 | k > 1/max(rev(model$space$psill)[1], rev(model$time$psill)[1])) 
    k <- 10^6*abs(k) # distorting the model to let optim "hopefully" find suitable parameters
  
  cbind(dist_grid, "gamma" = as.vector(vs+vt-k*vs*vt+vn))
}

covProdSumOld <- function(x, y, model) {
  .Deprecated("covProdSum", package = "gstat", 
              msg="The former product-sum model is dprecited, consider to refit the new model specification",
              old = "covProdSumOld")
  
  stopifnot(inherits(x, c("STF", "STS", "STI")) & inherits(y, c("STF", "STS", "STI")))
  
  # double check model for validity, i.e. k:
  k <- (sum(model$space$psill)+sum(model$time$psill)-model$sill)/(sum(model$space$psill)*sum(model$time$psill))
  if (k <= 0 | k > 1/max(model$space$psill[model$space$model!="Nug"], 
                         model$time$psill[model$time$model!="Nug"]))
    stop(paste("k (",k,") is non-positive or too large: no valid model!",sep=""))
  
  # the STF case
  if (inherits(x, "STF") & inherits(y, "STF")) {
    # calculate all spatial and temporal distances
    ds = spDists(x@sp, y@sp)
    dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    debug_time_unit(units(dt))
    dt <- as(dt, "matrix")
    
    # compose the cov-matrix
    vs = variogramLine(model$space, dist_vector = ds, covariance = TRUE)
    vt = variogramLine(model$time, dist_vector = dt, covariance = TRUE)
    
    return(model$sill-(vt %x% matrix(1,nrow(vs),ncol(vs)) + matrix(1,nrow(vt),ncol(vt)) %x% vs - k * vt %x% vs))
  } 
  
  # the STI case
  if(inherits(x, "STI") | inherits(y, "STI")) {
    # make sure that now both are of type STI
    x <- as(x, "STI")
    y <- as(y, "STI")
    
    # calculate all spatial and temporal distances
    ds = spDists(x@sp, y@sp)
    dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    debug_time_unit(units(dt))
    dt <- as(dt, "matrix")
    
    # compose the cov-matrix
    vs = variogramLine(model$space, dist_vector = ds, covariance = TRUE)
    vt = variogramLine(model$time, dist_vector = dt, covariance = TRUE)
    
    return(model$sill-(vt + vs - k * vt * vs))
  }
  
  # the remaining cases, none of x and y is STI nor are both STF
  # make sure both are of type STS
  x <- as(x, "STS")
  y <- as(y, "STS")
  
  # calculate all spatial and temporal distances
  ds = spDists(x@sp, y@sp)
  dt = abs(outer(index(x@time), index(y@time), "-"))
  if(!is.null(attr(model,"temporal unit")))
    units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
  debug_time_unit(units(dt))
  dt <- as(dt, "matrix")
  
  # re-arrange the spatial and temporal distances
  sMat <- matrix(NA, nrow(x@index), nrow(y@index))
  tMat <- matrix(NA, nrow(x@index), nrow(y@index))
  for(r in 1:nrow(x@index)) {
    sMat[r,] <- ds[x@index[r,1], y@index[,1]]
    tMat[r,] <- dt[x@index[r,2], y@index[,2]]
  }
  
  # compose the cov-matrix
  vs = variogramLine(model$space, dist_vector = sMat, covariance = TRUE)
  vt = variogramLine(model$time, dist_vector = tMat, covariance = TRUE)
  
  return(model$sill-(vt + vs - k * vt * vs))
}

# covariance for the circulant embedding in ST
covSurfProdSumOld <- function(model, dist_grid) {
  .Deprecated("covSurfProdSum", package = "gstat", 
              msg="The former product-sum model is dprecited, consider to refit the new model specification",
              old = "covSurfProdSumOld")
  
  vs = variogramLine(model$space, dist_vector = dist_grid$spacelag, covariance = TRUE)$gamma
  vt = variogramLine(model$time,  dist_vector = dist_grid$timelag, covariance = TRUE)$gamma
  
  k <- (sum(model$space$psill)+sum(model$time$psill)-(model$sill+model$nugget))/(sum(model$space$psill)*sum(model$time$psill))
  
  cbind(dist_grid, "gamma" = model$sill-(vt + vs - k * vt * vs))
}

vgmProdSum <- function(model, dist_grid) {
  if(!is.null(model$sill)) # backwards compatibility
    vgmProdSumOld(model, dist_grid)
  vs = variogramLine(model$space, dist_vector=dist_grid$spacelag)$gamma
  vt = variogramLine(model$time, dist_vector=dist_grid$timelag)$gamma
  
  sill_s <- sum(model$space$psill)
  sill_t <- sum(model$time$psill)
  k <- model$k
  
  cbind(dist_grid, "gamma" = as.vector((k*sill_t+1)*vs + (k*sill_s+1)*vt-k*vs*vt))
}

covProdSum <- function(x, y, model) {
  stopifnot(inherits(x, c("STF", "STS", "STI")) & inherits(y, c("STF", "STS", "STI")))
  if(!is.null(model$sill)) # backwards compatibility
    covProdSumOld(x, y, model)
  
  # the STF case
  if (inherits(x, "STF") & inherits(y, "STF")) {
    # calculate all spatial and temporal distances
    ds = spDists(x@sp, y@sp)
    dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    debug_time_unit(units(dt))
    dt <- as(dt, "matrix")
    
    # compose the cov-matrix
    vs = variogramLine(model$space, dist_vector = ds, covariance =TRUE)
    vt = variogramLine(model$time, dist_vector = dt, covariance =TRUE)
    
    return(vt %x% matrix(1,nrow(vs),ncol(vs)) + matrix(1,nrow(vt),ncol(vt)) %x% vs + model$k * vt %x% vs)
  } 
  
  # the STI case
  if(inherits(x, "STI") | inherits(y, "STI")) {
    # make sure that now both are of type STI
    x <- as(x, "STI")
    y <- as(y, "STI")
    
    # calculate all spatial and temporal distances
    ds = spDists(x@sp, y@sp)
    dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    debug_time_unit(units(dt))
    dt <- as(dt, "matrix")
    
    # compose the cov-matrix
    vs = variogramLine(model$space, dist_vector = ds, covariance=TRUE)
    vt = variogramLine(model$time, dist_vector = dt, covariance=TRUE)
    
    return(vt + vs + model$k * vt * vs)
  }
  
  # the remaining cases, none of x and y is STI nor are both STF
  # make sure both are of type STS
  x <- as(x, "STS")
  y <- as(y, "STS")
  
  # calculate all spatial and temporal distances
  ds = spDists(x@sp, y@sp)
  dt = abs(outer(index(x@time), index(y@time), "-"))
  if(!is.null(attr(model,"temporal unit")))
    units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
  debug_time_unit(units(dt))
  dt <- as(dt, "matrix")
  
  # re-arrange the spatial and temporal distances
  sMat <- matrix(NA, nrow(x@index), nrow(y@index))
  tMat <- matrix(NA, nrow(x@index), nrow(y@index))
  for(r in 1:nrow(x@index)) {
    sMat[r,] <- ds[x@index[r,1], y@index[,1]]
    tMat[r,] <- dt[x@index[r,2], y@index[,2]]
  }
  
  # compose the cov-matrix
  vs = variogramLine(model$space, dist_vector = sMat, covariance = TRUE)
  vt = variogramLine(model$time, dist_vector = tMat, covariance = TRUE)
  
  return(vt + vs + model$k * vt * vs)
}

# covariance for the circulant embedding in ST
covSurfProdSum <- function(model, dist_grid) {
  vs = variogramLine(model$space, dist_vector = dist_grid$spacelag, covariance = TRUE)$gamma
  vt = variogramLine(model$time, dist_vector = dist_grid$timelag, covariance = TRUE)$gamma
  
  cbind(dist_grid, "gamma" = vt + vs + model$k * vt * vs)
}

#########################################################
# sumMetric model: C_s + C_t + C_st (Gerard Heuvelink) ##
#########################################################

vgmSumMetric <- function(model, dist_grid) {
  vs = variogramLine(model$space, dist_vector=dist_grid$spacelag)$gamma
  vt = variogramLine(model$time,  dist_vector=dist_grid$timelag)$gamma
  h = sqrt(dist_grid$spacelag^2 + (model$stAni * as.numeric(dist_grid$timelag))^2)
  vst = variogramLine(model$joint, dist_vector=h)$gamma
  
  cbind(dist_grid, "gamma" = vs + vt + vst)
}

covSumMetric <- function(x, y, model) {
  stopifnot(inherits(x, c("STF", "STS", "STI")) & inherits(y, c("STF", "STS", "STI")))
  
  # the STF case
  if (inherits(x, "STF") & inherits(y, "STF")) {
    # calculate all spatial and temporal distances
    ds = spDists(x@sp, y@sp)
    dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    debug_time_unit(units(dt))
    dt <- as(dt, "matrix")
    
    # compose the cov-matrix
    Sm = variogramLine(model$space, covariance = TRUE, dist_vector = ds)
    Tm = variogramLine(model$time, covariance = TRUE, dist_vector = dt)
    
    h  = sqrt((matrix(1,nrow(dt),ncol(dt)) %x% ds)^2 
              + (model$stAni * dt %x% matrix(1,nrow(ds),ncol(ds)))^2)
    Mm = variogramLine(model$joint, covariance = TRUE, dist_vector = h)
    
    return(matrix(1,nrow(Tm),ncol(Tm)) %x% Sm + Tm %x% matrix(1,nrow(Sm),ncol(Sm)) + Mm)
  } 
  
  # the STI case
  if(inherits(x, "STI") | inherits(y, "STI")) {
    # make sure that now both are of type STI
    x <- as(x, "STI")
    y <- as(y, "STI")
    
    # calculate all spatial and temporal distances
    ds = spDists(x@sp, y@sp)
    dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    debug_time_unit(units(dt))
    dt <- as(dt, "matrix")
    
    # compose the cov-matrix
    Sm = variogramLine(model$space, covariance = TRUE, dist_vector = ds)
    Tm = variogramLine(model$time, covariance = TRUE, dist_vector = dt)
    
    h  = sqrt(ds^2 + (model$stAni * dt)^2)
    Mm = variogramLine(model$joint, covariance = TRUE, dist_vector = h)
    
    return(Sm + Tm + Mm)
  }
  
  # the remaining cases, none of x and y is STI nor are both STF
  # make sure both are of type STS
  x <- as(x, "STS")
  y <- as(y, "STS")
  
  # calculate all spatial and temporal distances
  ds = spDists(x@sp, y@sp)
  dt = abs(outer(index(x@time), index(y@time), "-"))
  if(!is.null(attr(model,"temporal unit")))
    units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
  debug_time_unit(units(dt))
  dt <- as(dt, "matrix")
  
  # re-arrange the spatial and temporal distances
  sMat <- matrix(NA, nrow(x@index), nrow(y@index))
  tMat <- matrix(NA, nrow(x@index), nrow(y@index))
  for(r in 1:nrow(x@index)) {
    sMat[r,] <- ds[x@index[r,1], y@index[,1]]
    tMat[r,] <- dt[x@index[r,2], y@index[,2]]
  }
  
  # compose the cov-matrix
  Sm = variogramLine(model$space, covariance = TRUE, dist_vector = sMat)
  Tm = variogramLine(model$time, covariance = TRUE, dist_vector = tMat)
  
  h  = sqrt(sMat^2 + (model$stAni * tMat)^2)
  Mm = variogramLine(model$joint, covariance = TRUE, dist_vector = h)
  
  return(Sm + Tm + Mm)
}

# covariance for the circulant embedding in ST
covSurfSumMetric <- function(model, dist_grid) {
  Sm = variogramLine(model$space, covariance = TRUE, dist_vector = dist_grid$spacelag)$gamma
  Tm = variogramLine(model$time, covariance = TRUE, dist_vector = dist_grid$timelag)$gamma
  
  h  = sqrt(dist_grid$spacelag^2 + (model$stAni * dist_grid$timelag)^2)
  Mm = variogramLine(model$joint, covariance = TRUE, dist_vector = h)$gamma
  
  cbind(dist_grid, "gamma" = Sm + Tm + Mm)
}

################################
## simplified sumMetric model ##
################################

vgmSimpleSumMetric <- function(model, dist_grid) {
  vs = variogramLine(model$space, dist_vector=dist_grid$spacelag)$gamma
  vt = variogramLine(model$time,  dist_vector=dist_grid$timelag)$gamma
  
  h = sqrt(dist_grid$spacelag^2 + (model$stAni * as.numeric(dist_grid$timelag))^2)
  
  vm = variogramLine(model$joint, dist_vector=h)$gamma
  vn <- variogramLine(vgm(model$nugget, "Nug", 0), dist_vector=h)$gamma
  
  cbind(dist_grid, "gamma" = vs + vt + vm + vn)
}

covSimpleSumMetric <- function(x, y, model) {
  modelNew <- vgmST("sumMetric", 
                    space=model$space, 
                    time=model$time,
                    joint=vgm(model$joint$psill[model$joint$model != "Nug"],
                              model$joint$model[model$joint$model != "Nug"],
                              model$joint$range[model$joint$model != "Nug"], 
                              model$nugget),
                    stAni=model$stAni)
  if (!is.null(attr(model,"temporal unit")))
    attr(modelNew,"temporal unit") <- attr(model,"temporal unit")
  covSumMetric(x, y, modelNew) 
}

# covariance for the circulant embedding in ST
covSurfSimpleSumMetric <- function(model, dist_grid) {
  modelNew <- vgmST("sumMetric", 
                    space=model$space, 
                    time=model$time,
                    joint=vgm(model$joint$psill[model$joint$model != "Nug"],
                              model$joint$model[model$joint$model != "Nug"],
                              model$joint$range[model$joint$model != "Nug"], 
                              model$nugget),
                    stAni=model$stAni)
  
  if (!is.null(attr(model,"temporal unit")))
    attr(modelNew,"temporal unit") <- attr(model,"temporal unit")
  
  covSurfSumMetric(modelNew, dist_grid) 
}

##################
## metric model ##
##################

vgmMetric <- function(model, dist_grid) {
  h = sqrt(dist_grid$spacelag^2 + (model$stAni * as.numeric(dist_grid$timelag))^2)

  cbind(dist_grid, "gamma" = variogramLine(model$joint, dist_vector=h)$gamma)
}

covMetric <- function(x, y, model) {
  stopifnot(inherits(x, c("STF", "STS", "STI")) & inherits(y, c("STF", "STS", "STI")))
  
  # the STF case
  if (inherits(x, "STF") & inherits(y, "STF")) {
    # calculate all spatial and temporal distances
    ds = spDists(x@sp, y@sp)
    dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    debug_time_unit(units(dt))
    dt <- as(dt, "matrix")
    
    # compose the cov-matrix
    h  = sqrt((matrix(1,nrow(dt),ncol(dt)) %x% ds)^2
              + (model$stAni * dt %x% matrix(1,nrow(ds),ncol(ds)))^2)
    Mm = variogramLine(model$joint, covariance = TRUE, dist_vector = h)
    
    return(Mm)
  } 
  
  # the STI case
  if(inherits(x, "STI") | inherits(y, "STI")) {
    # make sure that now both are of type STI
    x <- as(x, "STI")
    y <- as(y, "STI")
    
    # calculate all spatial and temporal distances
    ds = spDists(x@sp, y@sp)
    dt = abs(outer(index(x@time), index(y@time), "-"))
    if(!is.null(attr(model,"temporal unit")))
      units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
    debug_time_unit(units(dt))
    dt <- as(dt, "matrix")
    
    # compose the cov-matrix
    h  = sqrt(ds^2 + (model$stAni * dt)^2)
    Mm = variogramLine(model$joint, covariance = TRUE, dist_vector = h)
    
    return(Mm)
  }
  
  # the remaining cases, none of x and y is STI nor are both STF
  # make sure both are of type STS
  x <- as(x, "STS")
  y <- as(y, "STS")
  
  # calculate all spatial and temporal distances
  ds = spDists(x@sp, y@sp)
  dt = abs(outer(index(x@time), index(y@time), "-"))
  if(!is.null(attr(model,"temporal unit")))
    units(dt) <- attr(model, "temporal unit") # ensure the same temporal metric as in the variogram definition
  debug_time_unit(units(dt))
  dt <- as(dt, "matrix")
  
  # re-arrange the spatial and temporal distances
  sMat <- matrix(NA, nrow(x@index), nrow(y@index))
  tMat <- matrix(NA, nrow(x@index), nrow(y@index))
  for(r in 1:nrow(x@index)) {
    sMat[r,] <- ds[x@index[r,1], y@index[,1]]
    tMat[r,] <- dt[x@index[r,2], y@index[,2]]
  }
  
  # compose the cov-matrix
  h  = sqrt(sMat^2 + (model$stAni * tMat)^2)
  Mm = variogramLine(model$joint, covariance = TRUE, dist_vector = h)
  
  return(Mm)
}

# covariance for the circulant embedding in ST
covSurfMetric <- function(model, dist_grid) {
  h  = sqrt(dist_grid$spacelag^2 + (model$stAni * dist_grid$timelag)^2)
  
  cbind(dist_grid, "gamma" = variogramLine(model$joint, covariance = TRUE, dist_vector = h)$gamma)
}

###########################
## fitting ST variograms ##
###########################

fit.StVariogram <- function(object, model, ..., method = "L-BFGS-B", lower, upper, fit.method = 6, 
                            stAni=NA, wles) {
  if (!inherits(object, "StVariogram"))
    stop("\"object\" must be of class \"StVariogram\"")
  if (!inherits(model, "StVariogramModel"))
    stop("\"model\" must be of class \"StVariogramModel\".")
  
  sunit <- attr(object$spacelag, "units")
  tunit <- attr(object$timelag, "units")
  tu.obj = attr(model, "temporal unit")
  if (!is.null(tu.obj))
    stopifnot(identical(tunit, tu.obj))
  
  object$timelag = as.numeric(object$timelag) # needed for R 4.1
  object <- na.omit(object)
  
  ret <- model
  
  if(!missing(wles)) {
    if (wles)
      fit.method = 1
    else
      fit.method = 6
  }
  
  if (fit.method == 0) {
    attr(ret,"optim.output") <- "no fit"
    attr(ret, "MSE") <- mean((object$gamma - variogramSurface(model,
                                                              data.frame(spacelag=object$dist, timelag=object$timelag))$gamma)^2)
    attr(ret, "spatial unit")  <- sunit
    attr(ret, "temporal unit") <- tunit
    
    return(ret)
  }
  
  if ((fit.method == 7 | fit.method == 11) & is.null(model$stAni) & is.na(stAni)) {
    message("[An uninformed spatio-temporal anisotropy value of '1 (spatial unit)/(temporal unit)' is automatically selected. Consider providing a sensible estimate for stAni or using a different fit.method.]")
    stAni <- 1
  }
  
  weightingFun <- switch(fit.method,
                         function(obj, ...) obj$np, # 1
                         function(obj, gamma, ...) obj$np/gamma^2, # 2
                         function(obj, ...) obj$np, # 3
                         function(obj, gamma, ...) obj$np/gamma^2, # 4
                         function(obj, ...) stop("fit.method = 5 (REML), is not yet implemented"), # 5
                         function(obj, ...) 1, # 6
                         function(obj, curStAni, ...) 
                           if(is.na(stAni))
                             obj$np/(obj$dist^2+(curStAni*obj$timelag)^2)
                         else
                           obj$np/(obj$dist^2+(stAni*obj$timelag)^2), # 7
                         function(obj, ...) {
                           dist <- obj$dist
                           dist[dist == 0] <- min(dist[dist != 0], na.rm = TRUE)
                           obj$np/dist^2 # 8, pure space, 0 dist = min (dist > 0)
                         },
                         function(obj, ...) {
                           dist <- obj$timelag
                           dist[dist == 0] <- min(dist[dist != 0], na.rm = TRUE)
                           obj$np/dist^2
                         }, # 9, pure time
                         function(obj, gamma, ...) 1/gamma^2, # 10
                         function(obj, curStAni, ...) {
                           if(is.na(stAni))
                             1/(obj$dist^2+(curStAni*obj$timelag)^2)
                           else
                             1/(obj$dist^2+(stAni*obj$timelag)^2)
                         }, # 11
                         function(obj, ...) {
                           dist <- obj$dist
                           dist[dist == 0] <- min(dist[dist != 0], na.rm = TRUE)
                           1/(obj$dist^2) # 12, pure space
                         },
                         function(obj, ...) {
                           dist <- obj$timelag
                           dist[dist == 0] <- min(dist[dist != 0], na.rm = TRUE)
                           1/(obj$timelag^2)
                         }) # 13, pure time
  
  if(is.null(weightingFun))
    stop(paste("fit.method =", fit.method, "is not implementend"))
  
  fitFun = function(par, trace = FALSE, ...) {
    curModel <- insertPar(par, model)
    gammaMod <- variogramSurface(curModel,
                                 data.frame(spacelag=object$dist,
                                            timelag=object$timelag))$gamma
    resSq <- (object$gamma - gammaMod)^2
    resSq <- resSq * weightingFun(object, gamma=gammaMod, curStAni=curModel$stAni)
    if (trace)
      print(c(par, MSE = mean(resSq)))
    mean(resSq) # seems numerically more well behaved
  }
  
  if(missing(lower)) {
    min.s <- min(object$dist[object$dist>0])*0.05 # 5 % of the minimum distance larger 0
    min.t <- min(object$dist[object$timelag>0])*0.05 # 5 % of the minimum time lag 0),
    pos <- sqrt(.Machine$double.eps) # at least positive
    lower <- switch(strsplit(model$stModel, "_")[[1]][1],
                    separable=c(min.s, 0, min.t, 0, 0),
                    productSum=c(0, min.s, 0, 
                                 0, min.t, 0,
                                 pos),
                    productSumOld=c(0, min.s, 0, 
                                    0, min.t, 0, 0),
                    sumMetric=c(0, min.s, 0, 
                                0, min.t, 0,
                                0, pos, 0, pos),
                    simpleSumMetric=c(0, min.s,
                                      0, min.t,
                                      0, pos, 0, 0, pos),
                    metric=c(0, pos, 0, pos),
                    stop("Only \"separable\", \"productSum\", \"sumMetric\", \"simpleSumMetric\" and \"metric\" are implemented."))
  }
  if(missing(upper))
    upper <- switch(strsplit(model$stModel, "_")[[1]][1],
                    separable=c(Inf, 1, Inf, 1, Inf),
                    productSum=Inf,
                    productSumOld=Inf,
                    sumMetric=Inf,
                    simpleSumMetric=Inf,
                    metric=Inf,
                    stop("Only \"separable\", \"productSum\", \"sumMetric\", \"simpleSumMetric\" and \"metric\" are implemented."))
  
  
  
  pars.fit <- optim(extractPar(model), fitFun, ..., method = method, lower = lower, upper = upper)
  
  ret <- insertPar(pars.fit$par, model)
  attr(ret,"optim.output") <- pars.fit
  attr(ret, "MSE") <- mean((object$gamma - variogramSurface(insertPar(pars.fit$par, model),
                                                            data.frame(spacelag=object$dist, timelag=object$timelag))$gamma)^2)
  attr(ret, "spatial unit")  <- sunit
  attr(ret, "temporal unit") <- tunit
  
  return(ret)
}

###########
## tools ##
###########

# insert parameters into models
insertPar <- function(par, model) {
  switch(strsplit(model$stModel, "_")[[1]][1],
         separable=insertParSeparable(par, model),
         productSum=insertParProdSum(par, model),
         productSumOld=insertParProdSumOld(par, model),
         sumMetric=insertParSumMetric(par, model),
         simpleSumMetric=insertParSimpleSumMetric(par,model),
         metric=insertParMetric(par,model),
         stop("Only \"separable\", \"productSum\", \"sumMetric\", \"simpleSumMetric\" and \"metric\" are implemented."))
}

# extract parameters from models
extractPar <- function(model) {
  switch(strsplit(model$stModel, "_")[[1]][1],
         separable=c(range.s=model$space$range[2], nugget.s=model$space$psill[1],
                     range.t=model$time$range[2],  nugget.t=model$time$psill[1],
                     sill= model$sill[[1]]),
         productSumOld=c(sill.s = rev(model$space$psill)[1], range.s = rev(model$space$range)[1],
                         sill.t = rev(model$time$psill)[1],  range.t = rev(model$time$range)[1], 
                         sill=model$sill[[1]], nugget=model$nugget[[1]]),
         productSum=c(sill.s = model$space$psill[2], range.s = model$space$range[2], nugget.s = model$space$psill[1],
                      sill.t = model$time$psill[2],  range.t = model$time$range[2],  nugget.t = model$time$psill[1],
                      k=model$k),
         sumMetric=c(sill.s = model$space$psill[2], range.s = model$space$range[2], nugget.s = model$space$psill[1], 
                     sill.t = model$time$psill[2], range.t = model$time$range[2], nugget.t = model$time$psill[1],
                     sill.st = model$joint$psill[2], range.st = model$joint$range[2], nugget.st = model$joint$psill[1],
                     anis = model$stAni[[1]]),
         # simplified sumMetric model
         simpleSumMetric=c(sill.s = rev(model$space$psill)[1], range.s = rev(model$space$range)[1], 
                           sill.t = rev(model$time$psill)[1], range.t = rev(model$time$range)[1],
                           sill.st = rev(model$joint$psill)[1], range.st = rev(model$joint$range)[1], 
                           nugget = model$nugget[[1]], anis = model$stAni[[1]]),
         metric=c(sill = model$joint$psill[2], range = model$joint$range[2], nugget = model$joint$psill[1],
                  anis = model$stAni[[1]]),
         stop("Only \"separable\", \"productSum\", \"sumMetric\", \"simpleSumMetric\" and \"metric\" are implemented."))
}

# extract names
extractParNames <- function(model) {
  names(extractPar(model))
}

## dedicated insertion functions
################################

# separable model
insertParSeparable <- function(par, model) {
  vgmST("separable",
        space=vgm(1-par[2],as.character(model$space$model[2]),par[1],par[2],
                  kappa=model$space$kappa[2]),
        time= vgm(1-par[4],as.character(model$time$model[2]),par[3],par[4],
                  kappa=model$time$kappa[2]),
        sill=par[5])
}

# product sum model
insertParProdSumOld <- function(par, model) {
  vgmST("productSumOld",
        space=vgm(par[1],as.character(rev(model$space$model)[1]),par[2],
                  kappa=rev(model$space$kappa)[1]),
        time= vgm(par[3],as.character(rev(model$time$model)[1]),par[4],
                  kappa=rev(model$time$kappa)[1]),
        sill=par[5], nugget=par[6])
}

insertParProdSum <- function(par, model) {
  vgmST("productSum",
        space=vgm(par[1],as.character(model$space$model[2]),par[2],par[3],
                  kappa=model$space$kappa[2]),
        time= vgm(par[4],as.character(model$time$model[2]),par[5], par[6],
                  kappa=model$time$kappa[2]),
        k=par[7])
}

# sum metric model
insertParSumMetric <- function(par, model) {
  vgmST("sumMetric",
        space=vgm(par[1],as.character(model$space$model[2]),par[2],par[3],
                  kappa=model$space$kappa[2]),
        time= vgm(par[4],as.character(model$time$model[2]),par[5],par[6],
                  kappa=model$time$kappa[2]),
        joint=vgm(par[7],as.character(model$joint$model[2]),par[8],par[9],
                  kappa=model$joint$kappa[2]),
        stAni=par[10])
}

# simplified sum metric model
insertParSimpleSumMetric <- function(par, model) {
  vgmST("simpleSumMetric",
        space=vgm(par[1],as.character(rev(model$space$model)[1]),par[2],
                  kappa=rev(model$space$kappa)[1]),
        time= vgm(par[3],as.character(rev(model$time$model)[1]),par[4],
                  kappa=rev(model$time$kappa)[1]),
        joint=vgm(par[5],as.character(rev(model$joint$model)[1]),par[6],
                  kappa=rev(model$joint$kappa)[1]),
        nugget=par[7], stAni=par[8])
}

# metric model
insertParMetric <- function(par, model) {
  vgmST("metric",
        joint=vgm(par[1], as.character(model$joint$model[2]), par[2], par[3],
                  kappa=model$joint$kappa[2]),
        stAni=par[4])
}

Try the gstat package in your browser

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

gstat documentation built on Oct. 6, 2021, 5:06 p.m.