R/f.points.CKrige.R

In constrainedKriging: Constrained, Covariance-Matching Constrained and Universal Point or Block Kriging

f.points.CKrige<- function(
  formula,
  data,
  locations,
  object,
  method = 2,
  ex.out = FALSE
)
### purpose: calculate the constrained (C), covariance-matching constrained (CM)
###          and universal (U) kriging predictor for one block configuration
### ...........................................................................
### arguments:
###  pixconfig = list,
###                       1. element = coords of the lower left pixel
###                       2. element = vector with first number = number
###                          of pixel second number = number of polygons
###                       3. element = logical vector, TURE means the
###                          appropriate polyon area is smaller than the
###                          pixel area
###                       4. element = vector, numbers are the pixel centers
###                          in the appropiate polygon  e.g 0, 2, 1 means 0
###                          pixels in first poly  2 pixel in second poly
###                          and 1 pixel in
###                       5. element = matrix, coords of the centrois of the
###                          polygons
###                       6. element = list, each element is a number, the number
###                          describe the polygon in which the appropriate
###                          pixel center is, e.g 13th element hast number
###                          3 that means that the 13th pixel center is
###                            in the 3th polygon
###                       7. element = vector, numbers represent which polyon are
###                          considered e.g. c(6,3,10) means its a polygon
###                          configuration with the 6th, 3th and 10th polygon
### t.bb.cov = list, elements are block-block covariances matrices
### t.pred.designmat = n x number of predicton location matrix, model matrix of
###                 the prediction location
### model  = list, every element describ the autocorrelation
### t.cov.spline = Spline object of the  covariance between the sampling points
###                t.sample and the cells (rectangles)
### t.inv.Sigma = n x n inverse of the covariance matrix of the sampling points
###              t.sample
### t.list.beta.coef = list of two elements
###                        1. element = 1 x pi vector with the parameter of the
###                           general least square regression
###                        2. element = pi x pi matrix of the covariance between
###                           the general least square parameter
### t.Y =  n x 1 vector target value at the sampling points
### t.X = design matrix of the samples
### t.sample = n X 2 matrix with the coordinates of the smapling points,
###            first row are the x values and second row are the y values
### t.eps =  is the tuning parameter. It will replace those eigen values
###          when the positive definiteness of the Q or P matrix  is not
###          satisfied.
###
### ..............................................................................
# 2023-12-12 A. Papritz function returns gls residuals now as vector (and
#                       not asl one-column matrix)
{
  #
  t.support.designmat <- model.matrix(object = formula, data = data)
  #
  t.response  <- model.response( model.frame( formula, data ) )
  #
  locations <- terms( x= locations)
  attr(locations, "intercept") = 0
  locations <- model.matrix(object = locations, data = data)
  #
  if( dim( object@data )[1] == 0 && dim( object@data )[2] == 0 ) # ordinary case
  {
    t.pred.designmat <- matrix( rep( 1, dim( object@coords )[1] ) , ncol = 1 )
  }
  else
  {
    # t.terms.wr = terms without the response variable
    t.terms.wr <- delete.response( termobj = terms( x = formula ) )
    t.pred.designmat <- model.matrix( object = t.terms.wr, data = object@data)
    rm( t.terms.wr )
  }
  #
  ### covmodel of measurement free error process
  model.me.free <- object@model[unlist(lapply(1:length(object@model), function(i,m){m[[i]]$model != "mev"}, m = object@model))]
  ## Kovarianzmatrix der ST?tzpunkte
  t.support.covmat <- f.covmat.support( object@model, locations )

  # Inversere Choleskymatrix L^-1
  t.ichol <- forwardsolve(
    t.chol  <- t( chol(t.support.covmat ) ),
    diag( nrow( t.support.covmat ) )
  )

  # Inverse Kovarianzmatrix
  t.isupport.covmat <- crossprod( t.ichol, t.ichol )
  #orthogonalisierte Designmatrix und Messwerte der Datenpunkte
  t.orth.designmat <- t.ichol %*% t.support.designmat
  t.orth.data <- t.ichol %*% t.response

  # QR  <- Zerlegung der orthogonalisierten Designmatrix
  t.qr <- qr( as.data.frame( t.orth.designmat  )  )

  t.R <- qr.R( t.qr )
  t.iR <- backsolve( t.R, diag( nrow( t.R ) ) )
  # General least square Koeffizienten (beta GLS)
  t.beta.coef <- matrix( qr.coef( t.qr, as.data.frame( t.orth.data )), ncol = 1)

  # (Co)-Varianzen der beta Koeffizienten
  t.cov.beta.coef  <- crossprod( t( t.iR ), t( t.iR ) )
  # GLS-Residuen
  t.orth.resid <- qr.resid( t.qr, t.orth.data )
  t.residuals  <-  t.chol %*% t.orth.resid

  # lapply Laufvariable
  t.index  <- as.list( 1:dim( object@coords )[1] )

  t.krige.res <- lapply(t.index,
    function(i, coords, posindex, t.bb.covmat.list, t.pred.designmat,
      locations, model, t.ichol, t.isupport.covmat, t.beta.coef,
      t.cov.beta.coef , t.orth.resid, t.orth.designmat, t.iR, method)
    {

      t.indices <- posindex[[ i ]]

      coords <- matrix( coords[ t.indices, ], ncol= 2 )

      ## m X m covariance matrix of the block
      t.bb.covmat <- as.matrix( t.bb.covmat.list[[ t.indices[1] ]] )

      ## block design matrix
      t.pred.designmat <- matrix(t.pred.designmat[ t.indices, ],
        ncol = length( t.beta.coef ), nrow = length( t.indices ) )

      #cat("Parzelle: ", i, "\n")
      # Kovarianz zwischen den St?tzpunkten und den Zielbl?cken (n X m Matrix)
      t.dist <- f.row.dist( locations, coords )

      t.pred.covmat <- matrix(f.pp.cov( as.vector(t.dist), model), ncol = ncol(t.dist) )

      rm( t.dist )

      t.pred.covmat.ichol.trans <-  crossprod( t.pred.covmat, t( t.ichol) )

      # t(c) %*% Sigma^-1 = t(c) %*% t(L^-1) %*% L^1
      t.sk.weights  <- t.pred.covmat.ichol.trans %*% t.ichol


      ## least square prediction

      t.lin.trend.est <-  crossprod( t( t.pred.designmat), t.beta.coef )


      ## residuals weithed by the simple kriging weights

      t.weighted.resid <- crossprod( t( t.pred.covmat.ichol.trans ), t.orth.resid )

      ## UK Vorhersage
      t.uk.pred <- t.lin.trend.est  + t.weighted.resid

      # UK MSPE
      t.aux <- crossprod( t( t.pred.designmat - crossprod(t(t.pred.covmat.ichol.trans), t.orth.designmat) ) , t.iR)
      t.uk.mspe <- t.bb.covmat - tcrossprod( t.pred.covmat.ichol.trans, t.pred.covmat.ichol.trans ) + t.aux %*% t( t.aux )

      #print(t.bb.covmat - tcrossprod( t.pred.covmat.ichol.trans, t.pred.covmat.ichol.trans ))

      # Kriging Vorhersage und MSPE
      t.result <- f.kriging.method( t.lin.trend.est, t.cov.beta.coef, t.weighted.resid,
        t.uk.mspe, t.bb.covmat, t.pred.designmat,
        t.orth.designmat, t.pred.covmat.ichol.trans, method)

      if(method == 1)
      {
        return( list( prediction = t.result[[1]], sqrt.mspe = t.result[[2]], sk.weights = t.sk.weights ) )
      }
      if( method == 2 )
      {
        return( list( prediction = t.result[[1]], sqrt.mspe = t.result[[2]],
            P1 = t.result[[3]], Q1 =t.result[[4]], K = t.result[[5]], sk.weights = t.sk.weights) )
      }
      if( method == 3)
      {
        return( list( prediction = t.result[[1]], mspe = t.result[[2]],
            P1 = t.result[[3]], Q1 =t.result[[4]], K = t.result[[5]], sk.weights = t.sk.weights ) )
      }

    },
    coords = object@coords,
    posindex = object@posindex,
    t.bb.covmat.list = object@covmat,
    t.pred.designmat = t.pred.designmat,
    locations = locations,
    model = model.me.free,
    t.ichol = t.ichol,
    t.isupport.covmat = t.isupport.covmat,
    t.beta.coef = t.beta.coef,
    t.cov.beta.coef = t.cov.beta.coef,
    t.orth.resid = t.orth.resid,
    t.orth.designmat = t.orth.designmat,
    t.iR = t.iR,
    method = method
  ) ## end of lapply
  ##
  ## output
  ##
  # # #
  if( method == 1 ){
    krige.result <- data.frame( matrix( NaN, nrow = nrow( t.pred.designmat ), ncol = 2 ) )
    krige.result[,1] <- unlist( lapply( t.krige.res, function( poly ){ return( poly$prediction ) } ) )
    krige.result[,2] <- unlist( lapply( t.krige.res, function( poly ){ return( poly$sqrt.mspe ) } ) )
    colnames( krige.result) = c("prediction", "prediction.se")
  }
  #
  if( method == 2 )
  {
    krige.result <- data.frame( matrix( NaN, nrow = nrow(t.pred.designmat), ncol = 5 ) )
    krige.result[,1] <- unlist( lapply( t.krige.res, function( poly ){ return( poly$prediction ) } ) )
    krige.result[,2] <- unlist( lapply( t.krige.res, function( poly ){ return(  poly$sqrt.mspe  ) } ) )
    krige.result[,3] <- unlist( lapply( t.krige.res, function( poly ){ return( poly$P1 ) } ) )
    krige.result[,4] <- unlist( lapply( t.krige.res, function( poly ){ return( poly$Q1 ) } ) )
    krige.result[,5] <- unlist( lapply( t.krige.res, function( poly ){ return( poly$K ) } ) )
    colnames( krige.result) = c("prediction", "prediction.se", "P1", "Q1", "K")
  }
  #
  if( method == 3)
  {
    krige.result <- data.frame( matrix( NaN, nrow = nrow(t.pred.designmat), ncol = 5 ) )
    krige.result[,1] <- unlist( lapply( t.krige.res, function( poly ){ return( poly$prediction[1,1] ) } ) )
    krige.result[,2] <- unlist( lapply( t.krige.res, function( poly ){ return( sqrt( poly$mspe[1,1] ) ) } ) )
    krige.result[,3] <- unlist( lapply( t.krige.res, function( poly ){ return( poly$P1[1,1] ) } ) )
    krige.result[,4] <- unlist( lapply( t.krige.res, function( poly ){ return( poly$Q1[1,1] ) } ) )
    krige.result[,5] <- unlist( lapply( t.krige.res, function( poly ){ return( poly$K[1,1] ) } ) )
    colnames( krige.result) = c("prediction", "prediction.se","P1.11", "Q1.11", "K.11")
  }
  #
  if(  ex.out == TRUE)
  {
    if( method == 1 | method == 2 )
    {
      #sk.weights.matrix
      if( is.null( dim(t.krige.res[[1]]$sk.weights) ) )
      {
        sk.weights.matrix <- t(matrix( unlist( lapply( t.krige.res, function( points ){ return( points$sk.weights ) } ) ), ncol = nrow( data ), byrow = TRUE ))
      }else{
        sk.weights.matrix <- t(matrix( unlist( lapply( t.krige.res, function( points ){ return( points$sk.weights[1,] ) } ) ), ncol = nrow( data ), byrow = TRUE))
      }
      ## extended output is a list
      object <- SpatialPointsDataFrame( SpatialPoints( object@coords ), data = krige.result, match.ID = FALSE)
      res <- list(
        object = object,
        krig.method = method,
        parameter = list( beta.coef = t.beta.coef, cov.beta.coef = t.cov.beta.coef),
        sk.weights = sk.weights.matrix,
        inv.Sigma = t.isupport.covmat,
        residuals = c(t.residuals)
      )
      class( res )  <- "CKrige.exout.points"
      return( res )
    }
    else
    {
      object <- SpatialPointsDataFrame( SpatialPoints(object@coords ), data = krige.result, match.ID = FALSE)
      P1.list = lapply( t.krige.res, function( points ){ return(points$P1)})
      Q1.list = lapply( t.krige.res, function( points ){ return(points$Q1)})
      K.list = lapply( t.krige.res, function( points ){ return(points$K)})
      sk.weights.list <-lapply( t.krige.res, function( points ){ return(t(points$sk.weights))})
      res <- list(
        object = object,
        krig.method = method,
        CMCK.par = list(P1 = P1.list, Q1 =Q1.list,K = K.list),
        parameter = list( beta.coef = t.beta.coef, cov.beta.coef = t.cov.beta.coef),
        sk.weights = sk.weights.list,
        inv.Sigma = t.isupport.covmat,
        residuals = c(t.residuals)
      )
      class( res )  <- "CKrige.exout.points"
      return( res )
    }
  }
  else
  {

    row.names( krige.result ) <- row.names( object@coords )
    object <- SpatialPointsDataFrame( SpatialPoints( object@coords ), krige.result )
    return( object )
  }
  #
} # end of function