fgeo.krig: Analyze Soils Data

#' Version of geoR::ksline that avoids a warning.
#'
#' See https://github.com/forestgeo/krig/issues/9.
#' @seealso [geoR::ksline()]
#' @keywords internal
#' @noRd
krig_ksline <- function(geodata, coords = geodata$coords, data = geodata$data,
                        locations, borders = NULL, cov.model = "matern", cov.pars = stop("covariance parameters (sigmasq and phi) needed"),
                        kappa = 0.5, nugget = 0, micro.scale = 0, lambda = 1, m0 = "ok",
                        nwin = "full", n.samples.backtransform = 500, trend = 1,
                        d = 2, ktedata = NULL, ktelocations = NULL, aniso.pars = NULL,
                        signal = FALSE, dist.epsilon = 1e-10, messages) {
  if (missing(messages)) {
    messages.screen <- as.logical(ifelse(is.null(getOption("geoR.messages")),
      TRUE, getOption("geoR.messages")
    ))
  } else {
    messages.screen <- messages
  }
  locations <- .check.locations(locations)
  if (length(unique(locations[, 1])) == 1 |
      length(unique(locations[, 2])) == 1) {
    krige1d <- TRUE
  } else {
    krige1d <- FALSE
  }
  if (!is.null(borders)) {
    locations <- locations.inside(locations, borders, bound = TRUE)
    if (messages.screen) {
      cat("ksline: results will be returned only for prediction locations inside the borders\n")
    }
  }
  call.fc <- match.call()
  cov.model <- match.arg(cov.model, choices = .geoR.cov.models)
  if (abs(lambda - 1) > 0.001) {
    if (messages.screen) {
      cat("ksline: Box-Cox data transformation performed.\n")
    }
    if (abs(lambda) < 0.001) {
      data <- log(data)
    } else {
      data <- ((data^lambda) - 1) / lambda
    }
  }
  coords <- as.matrix(coords)
  locations <- as.matrix(locations)
  dimnames(coords) <- list(NULL, NULL)
  dimnames(locations) <- list(NULL, NULL)
  if (!is.null(ktedata) & !is.null(ktelocations) & m0 != "kte") {
    cat("ksline: external variable (covariate) provided. Kriging ste to KTE\n")
    m0 <- "kte"
  }
  if (!is.null(aniso.pars)) {
    if (length(aniso.pars) != 2 | mode(aniso.pars) != "numeric") {
      stop("anisotropy parameters must be provided as a numeric vector with two elements: the rotation angle (in radians) and the anisotropy ratio (a number greater than 1)")
    }
    if (messages.screen) {
      cat("ksline: anisotropy correction performed\n")
    }
    coords.c <- coords.aniso(coords = coords, aniso.pars = aniso.pars)
    locations.c <- coords.aniso(coords = locations, aniso.pars = aniso.pars)
  }
  else {
    coords.c <- coords
    locations.c <- locations
  }
  if (m0 == "kte") {
    ktedata <- as.matrix(ktedata)
    ktelocations <- as.matrix(ktelocations)
    dimnames(ktedata) <- list(NULL, NULL)
    dimnames(ktelocations) <- list(NULL, NULL)
  }
  n <- length(data)
  ni <- length(locations[, 1])
  tausq <- nugget
  sigmasq <- cov.pars[1]
  phi <- cov.pars[2]
  if (nwin == "full") {
    est <- rep(0, ni)
    dif <- rep(0, ni)
    kvar <- rep(0, ni)
    sumw <- rep(0, ni)
    wofmean <- rep(0, ni)
    iv <- varcov.spatial(
      coords = coords.c, cov.model = cov.model,
      kappa = kappa, nugget = nugget, cov.pars = cov.pars,
      inv = TRUE, det = FALSE, func.inv = "cholesky"
    )$inverse
    av <- mean(data)
    sd <- sqrt(var(data))
    one <- rep(1, n)
    tone <- t(one)
    toneiv <- crossprod(one, iv)
    den <- solve(toneiv %*% one)
    ml <- den %*% toneiv %*% data
    kmsd <- sqrt(den)
    means <- c(average = av, stdev = sd, kmean = ml, kmsd = kmsd)
    if (m0 != "kt") {
      mktlocations <- "Constant trend"
      beta <- ml
    }
    else {
      mktlocations <- rep(0, ni)
      if (m0 == "kt" & trend == 1) {
        if (d == 1) {
          xmat <- cbind(rep(1, n), coords[, 2])
          xmati <- cbind(rep(1, ni), locations[, 2])
        }
        else {
          xmat <- cbind(rep(1, n), coords[, 1], coords
          [
            ,
            2
          ])
          xmati <- cbind(
            rep(1, ni), locations[, 1],
            locations[, 2]
          )
        }
        iviv <- solve(crossprod(xmat, iv) %*% xmat)
        txiv <- crossprod(xmat, iv)
        beta <- iviv %*% txiv %*% data
        mkt <- xmat %*% beta
      }
      if (m0 == "kt" & trend == 2) {
        if (d == 1) {
          xmat <- cbind(rep(1, n), coords[, 2], (coords
          [
            ,
            2
          ])^2)
          xmati <- cbind(
            rep(1, ni), locations[, 2],
            (locations[, 2])^2
          )
        }
        else {
          xmat <- cbind(rep(1, n), coords[, 1], coords
          [
            ,
            2
          ], (coords[, 1])^2, (coords[, 2])^2, coords
          [
            ,
            1
          ] * coords[, 2])
          xmati <- cbind(
            rep(1, ni), locations[, 1],
            locations[, 2], (locations[, 1])^2, (locations
            [
              ,
              2
            ])^2, locations[, 1] * locations[, 2]
          )
        }
        iviv <- solve(crossprod(xmat, iv) %*% xmat)
        txiv <- crossprod(xmat, iv)
        beta <- iviv %*% txiv %*% data
        mkt <- xmat %*% beta
      }
    }
    if (m0 != "kte") {
      mktelocations <- "No external trend"
    } else {
      if (m0 == "kte") {
        mktelocations <- rep(0, ni)
        xmat <- cbind(rep(1, n), ktedata)
        xmati <- cbind(rep(1, ni), ktelocations)
        iviv <- solve(crossprod(xmat, iv) %*% xmat)
        txiv <- crossprod(xmat, iv)
        beta <- iviv %*% txiv %*% data
        mkte <- xmat %*% beta
      }
    }
    for (i in 1:ni) {
      if (messages.screen) {
        if (ni < 11) {
          cat(paste(
            "ksline: kriging location: ", i,
            "out of", ni, "\n"
          ))
        } else {
          if (ni < 101 & (i %% 10 == 1)) {
            cat(paste(
              "ksline: kriging location: ",
              i, "out of", ni, "\n"
            ))
          }
          if (ni > 100 & i %% 100 == 1) {
            cat(paste(
              "ksline: kriging location: ",
              i, "out of", ni, "\n"
            ))
          }
          if (i == ni) {
            cat(paste(
              "ksline: kriging location: ",
              i, "out of", ni, "\n"
            ))
          }
        }
      }
      coords0 <- cbind((coords.c[, 1] - locations.c[
        i,
        1
      ]), (coords.c[, 2] - locations.c[i, 2]))
      dm0 <- sqrt(coords0[, 1]^2 + coords0[, 2]^2)
      v0 <- cov.spatial(
        obj = dm0, cov.model = cov.model,
        kappa = kappa, cov.pars = cov.pars
      )
      v0[dm0 < dist.epsilon] <- micro.scale + sigmasq
      tv0 <- t(v0)
      v0iv <- crossprod(v0, iv)
      v0ivv0 <- v0iv %*% v0
      skw <- crossprod(v0, iv)
      wofmean[i] <- 1 - sum(skw)
      if (mode(m0) == "numeric") {
        dif[i] <- skw %*% (data - m0)
        est[i] <- m0 + dif[i]
        if (signal == TRUE) {
          kvar[i] <- sigmasq - v0ivv0
        } else {
          kvar[i] <- tausq + sigmasq - v0ivv0
        }
        sumw[i] <- sum(skw)
      }
      if (m0 == "av") {
        dif[i] <- skw %*% (data - av)
        est[i] <- av + dif[i]
        if (signal == TRUE) {
          kvar[i] <- sigmasq - v0ivv0
        } else {
          kvar[i] <- tausq + sigmasq - v0ivv0
        }
        sumw[i] <- sum(((tone / n) + skw - ((skw %*% one %*%
          tone) / n)))
      }
      if (m0 == "ok") {
        dif[i] <- skw %*% (data - as.vector(ml))
        est[i] <- ml + dif[i]
        redu <- as.vector(1 - toneiv %*% v0)
        if (signal == TRUE) {
          kvar[i] <- sigmasq - v0ivv0 + (redu %*% den %*%
            redu)
        } else {
          kvar[i] <- tausq + sigmasq - v0ivv0 + (redu %*%
            den %*% redu)
        }
        sumw[i] <- sum((den %*% one + tv0 - v0iv %*%
          one %*% den %*% tone) %*% iv)
      }
      if (m0 == "kt") {
        dif[i] <- skw %*% (data - mkt)
        est[i] <- xmati[i, ] %*% beta + dif[i]
        redu <- as.vector(xmati[i, ]) - as.vector(txiv %*%
          v0)
        if (signal == TRUE) {
          kvar[i] <- sigmasq - v0ivv0 + (redu %*% iviv %*%
            redu)
        } else {
          kvar[i] <- tausq + sigmasq - v0ivv0 + (redu %*%
            iviv %*% redu)
        }
        sumw[i] <- sum(skw + xmati[i, ] %*% iviv %*%
          txiv - skw %*% xmat %*% iviv %*% txiv)
        mktlocations[i] <- xmati[i, ] %*% beta
      }
      if (m0 == "kte") {
        dif[i] <- skw %*% (data - mkte)
        est[i] <- xmati[i, ] %*% beta + dif[i]
        redu <- as.vector(xmati[i, ]) - as.vector(txiv %*%
          v0)
        if (signal == TRUE) {
          kvar[i] <- sigmasq - v0ivv0 - (redu %*% iviv %*%
            redu)
        } else {
          kvar[i] <- tausq + sigmasq - v0ivv0 + (redu %*%
            iviv %*% redu)
        }
        sumw[i] <- sum(skw + xmati[i, ] %*% iviv %*%
          txiv - skw %*% xmat %*% iviv %*% txiv)
        mktelocations[i] <- xmati[i, ] %*% beta
      }
      NULL
    }
    message <- "Kriging performed using global neighbourhood"
    if (messages.screen) {
      cat(paste(message, "\n"))
    }
    results <- list(
      predict = est, krige.var = kvar, dif = dif,
      summary = means, ktrend = mktlocations, ktetrend = mktelocations,
      beta = beta, wofmean = wofmean
    )
  }
  else {
    nwin <- min(n, nwin)
    avwin <- rep(0, ni)
    sdwin <- rep(0, ni)
    mlwin <- rep(0, ni)
    kmsdwin <- rep(0, ni)
    estwin <- rep(0, ni)
    difwin <- rep(0, ni)
    kvarwin <- rep(0, ni)
    sumwwin <- rep(0, ni)
    wofmean <- rep(0, ni)
    if (m0 != "kt") {
      mkt <- "Constant position trend"
    } else {
      mkt <- rep(0, ni)
    }
    if (m0 != "kte") {
      mkte <- "No external trend"
    } else {
      mkte <- rep(0, ni)
    }
    if (m0 != "kt" & m0 != "kte") {
      betawin <- "No polynomial or external trend"
    }
    if (m0 == "kt") {
      if (trend == 1) {
        if (d == 1) {
          xmati <- cbind(rep(1, ni), locations[, 2])
        } else {
          xmati <- cbind(rep(1, ni), locations
          [
            ,
            1
          ], locations[, 2])
        }
      }
      if (trend == 2) {
        if (d == 1) {
          xmati <- cbind(
            rep(1, ni), locations[, 2],
            locations[, 2]^2
          )
        } else {
          xmati <- cbind(rep(1, ni), locations
          [
            ,
            1
          ], locations[, 2], (locations[, 1])^2, (locations
          [
            ,
            2
          ])^2, locations[, 1] * locations[, 2])
        }
      }
      betawin <- matrix(0,
        nrow = (ncol(xmati) * ni),
        ncol = ncol(xmati)
      )
    }
    if (m0 == "kte") {
      xmati <- cbind(rep(1, ni), ktelocations)
      if (is.vector(ktedata) == TRUE) {
        betawin <- matrix(0, nrow = (2 * ni), ncol = 2)
      } else {
        betawin <- matrix(0, nrow = ((ncol(ktedata) +
          1) * ni), ncol = (ncol(ktedata) + 1))
      }
    }
    for (i in 1:ni) {
      temp.win <- .ksline.aux.1(
        coords = coords, coords.c = coords.c,
        data = data, n = n, locations = locations[i, ], locations.c = locations.c[i, ], cov.pars = cov.pars,
        nugget = nugget, cov.model = cov.model, kappa = kappa,
        m0 = m0, nwin = nwin, trend = trend, d = d,
        ktedata = ktedata, ktelocations = ktelocations,
        micro.scale = micro.scale, location.number = i,
        xmati = xmati[i, ], mkte = NULL, mkt = NULL,
        betawin = NULL, signal = signal, dist.epsilon = dist.epsilon
      )
      avwin[i] <- temp.win$avwin
      sdwin[i] <- temp.win$sdwin
      mlwin[i] <- temp.win$mlwin
      kmsdwin[i] <- temp.win$kmsdwin
      estwin[i] <- temp.win$estwin
      difwin[i] <- temp.win$difwin
      kvarwin[i] <- temp.win$kvarwin
      sumwwin[i] <- temp.win$sumwwin
      wofmean[i] <- temp.win$wofmean
      if (m0 == "kt") {
        mkt[i] <- temp.win$mkt
      }
      if (m0 == "kte") {
        mkte[i] <- temp.win$mkte
      }
      if (m0 == "kt" | m0 == "kte") {
        betawin[i, ] <- temp.win$betawin
      }
      if (messages.screen) {
        if (ni < 11) {
          cat(paste(
            "ksline: kriging location: ", i,
            "out of", ni, "\n"
          ))
        } else {
          if (ni < 101 & (i %% 10 == 1)) {
            cat(paste(
              "ksline: kriging location: ",
              i, "out of", ni, "\n"
            ))
          }
          if (ni > 100 & i %% 100 == 1) {
            cat(paste(
              "ksline: kriging location: ",
              i, "out of", ni, "\n"
            ))
          }
          if (i == ni) {
            cat(paste(
              "ksline: kriging location: ",
              i, "out of", ni, "\n"
            ))
          }
        }
      }
    }
    message <- "kriging performed in moving neighbourhood"
    if (messages.screen) {
      cat(paste(message, "\n"))
    }
    results <- list(
      predict = estwin, krige.var = kvarwin,
      dif = difwin, avtrend = avwin, sd = sdwin, oktrend = mlwin,
      oksd = kmsdwin, ktrend = mkt, ktetrend = mkte, beta = betawin,
      sumw = sumwwin, wofmean = wofmean
    )
  }
  if (abs(lambda - 1) > 0.001) {
    if (messages.screen) {
      cat("Back-transforming the predicted mean and variance.\n")
    }
    if (abs(lambda) < 0.001) {
      predict.transf <- results$predict
      results$predict <- exp(predict.transf) - 0.5 * results$krige.var
      results$krige.var <- (exp(2 * predict.transf - results$krige.var)) *
        (exp(results$krige.var) - 1)
    }
    if (lambda > 0.001) {
      if (messages.screen) {
        cat("Back-transformation by simulating from the normal predictive distribution\n")
      }
      temp.data <- suppressWarnings(matrix(rnorm(ni *
        n.samples.backtransform,
      mean = results$predict,
      sd = sqrt(results$krige.var)
      ), nrow = ni))
      temp.data[(results$krige.var == 0), ] <- results$predict[(results$krige.var ==
        0)]
      temp.data[temp.data < -1 / lambda] <- -1 / lambda
      temp.data <- ((temp.data * lambda) + 1)^(1 / lambda)
      results$predict <- as.vector(rowMeans(temp.data))
      results$krige.var <- as.vector(apply(
        temp.data,
        1, var
      ))
    }
    if (lambda < -0.001) {
      cat("Resulting distribution has no mean for lambda < 0 - back transformation not performed\n")
    }
  }
  results$locations <- locations
  results$message <- message
  results$call <- call.fc
  attr(results, "sp.dim") <- ifelse(krige1d, "1d", "2d")
  attr(results, "prediction.locations") <- call.fc$locations
  if (!is.null(call.fc$borders)) {
    attr(results, "borders") <- call.fc$borders
  }
  oldClass(results) <- c("kriging")
  return(invisible(results))
}