GeoKrigloc: Spatial (bivariate) and spatio temporal optimal linear local...

In GeoModels: Procedures for Gaussian and Non Gaussian Geostatistical (Large) Data Analysis

Spatial (bivariate) and spatio temporal optimal linear local prediction for Gaussian and non-Gaussian random fields.

Description

For a given set of spatial location sites (and temporal instants), the function computes optimal local linear prediction and the associated mean squared error for the Gaussian and non-Gaussian case using a spatial (temporal) neighborhood computed using the function GeoNeighborhood

Usage

GeoKrigloc(estobj=NULL,data, coordx=NULL, coordy=NULL, coordz=NULL,coordt=NULL,
 coordx_dyn=NULL, corrmodel, distance="Eucl", grid=FALSE, 
 loc, neighb=NULL, maxdist=NULL, 
 maxtime=NULL, method="cholesky",
 model="Gaussian", n=1,nloc=NULL, mse=FALSE, 
 param, anisopars=NULL,radius=1,
 sparse=FALSE, time=NULL, type="Standard",
 type_krig="Simple",weigthed=TRUE, 
 which=1, copula=NULL,X=NULL,Xloc=NULL,Mloc=NULL,varcov=NULL,
 spobj=NULL,spdata=NULL,parallel=FALSE,ncores=6,progress=TRUE,
 check.duplicates=FALSE)

Arguments

estobj	An object of class Geofit that includes information about data, model and estimates.
data	A d-dimensional vector (a single spatial realisation) or a (d \times d)-matrix (a single spatial realisation on regular grid) or a (t \times d)-matrix (a single spatio-temporal realisation) or an (d \times d \times t \times n )-array (a single spatio-temporal realisation on regular grid) giving the data used for prediction.
coordx	A numeric d \times 2 or d \times 3 coordinate matrix. It may be omitted when coordx_dyn or spobj is supplied. Coordinates on a sphere for a fixed radius radius are passed in lon/lat format expressed in decimal degrees.
coordy	A numeric vector giving 1-dimension of spatial coordinates; Optional argument, the default is NULL.
coordz	A numeric vector giving 1-dimension of spatial coordinates; Optional argument, the default is NULL.
coordt	A numeric vector giving the temporal coordinates of the observations. The default is NULL, in which case a spatial random field is expected. Temporal coordinates may be irregularly spaced; temporal lags are computed from the supplied coordinate values.
coordx_dyn	For dynamic observation locations, a list with one two- or three-column coordinate matrix per element of coordt. The rows of coordx_dyn[[t]] correspond to data[[t]] and to the same temporal block of X. Dynamic observation locations do not alter the ordering of prediction tasks. See GeoModels-spacetime-ordering.
corrmodel	String; the name of a correlation model, for the see GeoCovmatrix for the list of implemented correlation models.
distance	String; the name of the spatial distance. The default is Eucl, the Euclidean distance. See GeoFit for details.
grid	Logical; if FALSE (the default) the data used for prediction are interpreted as spatial or spatio-temporal realisations on a set of non-equispaced spatial sites (irregular grid).
loc	A numeric (n \times 2)-matrix (where n is the number of spatial sites) giving 2-dimensions of spatial coordinates to be predicted.
neighb	Numeric; an optional positive integer indicating the order of the neighborhood.
maxdist	Numeric; an optional positive value indicating the distance in the spatial neighborhood.
maxtime	Numeric; an optional non-negative maximum temporal-distance threshold, expressed in the same units as coordt, used for the local temporal neighborhood.
method	String; matrix decomposition used to solve each local kriging system. The choices are cholesky (default) and svd. Only cholesky is available when sparse=TRUE.
n	Positive integer size parameter. For Binomial it may be scalar or contain one value per observation. For direct Negative Binomial it is the single common number r of successes; Geometric corresponds to r=1. Default is 1.
nloc	Positive integer size at prediction tasks. For Binomial, site-specific observation sizes require explicit prediction sizes. For direct Negative Binomial, nloc, if supplied, must equal the common r; for Geometric it must equal 1. No prediction size is inferred by averaging observation-side values.
mse	Logical; if TRUE (the default) MSE of the kriging predictor is computed. For model="Wrapped", the returned MSE is the sum of the sine- and cosine-component linear-prediction MSEs rather than an ordinary squared angular error.
model	String; the type of RF and therefore the densities associated to the likelihood objects. Gaussian is the default, see the Section Details.
param	A list of parameter values required for the correlation model. See Details for the accepted options.
anisopars	A list of two elements: "angle" and "ratio" i.e. the anisotropy angle and the anisotropy ratio, respectively.
radius	Numeric: the radius of the sphere if coordinates are passed in lon/lat format;Default value is 1.
sparse	Logical; if TRUE kriging is computed with sparse matrices algorithms using spam package. Default is FALSE. It should be used with compactly supported covariances.
time	A numeric vector giving the temporal instants to be predicted. Values need not be equally spaced and are interpreted on the same numeric time scale as coordt. The default is NULL, in which case only spatial prediction is performed.
type	String; currently only Standard is computed. Other legacy values produce a warning and are treated as Standard.
type_krig	String; Simple (default) treats fitted mean coefficients as known. Universal uses the same coefficients and adds their uncertainty from varcov to each local prediction MSE.
weigthed	Logical legacy argument retained for backward compatibility. It has no effect for type="Standard".
which	Numeric; In the case of bivariate (tapered) cokriging it indicates which variable to predict. It can be 1 or 2
copula	String; optional copula specification. Local linear prediction is implemented for "Gaussian", "Clayton", and "SkewGaussian".
X	Numeric design matrix at the observations. For fixed-location space-time data, rows follow c(t(data)), i.e. time then site. For dynamic observations, use a stacked matrix ordered by temporal blocks, or a list aligned with coordx_dyn.
Xloc	Numeric design matrix at prediction tasks. For space-time prediction its rows are location-major: all requested times for loc[1, ], then all requested times for loc[2, ], and so on.
Mloc	Numeric vector giving the known marginal location predictor at prediction tasks, in the same location-major order as Xloc. For Gaussian additive margins this is the mean itself; for models with link-scale parametrizations (for example Gamma/Weibull/LogGaussian or Beta2) it is on the same scale as param$mean. Use either Mloc or Xloc, not both.
varcov	Covariance matrix of the estimated parameters, required for type_krig="Universal" when mse=TRUE. With a composite fit it should be the sandwich covariance (inverse Godambe information). It is normally inherited from estobj.
spobj	An object of class sp or spacetime. For space-time objects, the current sp2Geo() conversion uses sequential temporal indices and does not preserve irregular original time spacing; to retain irregular temporal distances, use the explicit-coordinate interface with numeric coordt.
spdata	Character:The name of data in the sp or spacetime object
parallel	Logical; default FALSE. If TRUE, the independent local GeoKrig calls are evaluated with future.apply.
ncores	Numeric or NULL; default 6. With parallel=TRUE, a numeric value requests that many worker processes, capped only by the number of prediction tasks and detected cores. Thus the default parallel execution uses up to six workers. Set ncores=NULL to request the automatic/safe mode, in which the package-wide automatic cap of six workers and the internal RAM safety check may further reduce the worker count.
progress	If TRUE then a progress bar is shown.
check.duplicates	Logical. If TRUE, perform a fast exact scan for duplicated observation locations at this user entry point. The default is FALSE, so no duplicate-location scan is imposed. Internal bootstrap, cross-validation, and refitting calls do not repeat the scan.

Details

The optional pair cache is bounded by getOption("GeoModels.local_pair_cache_max_bytes", 256 * 1024^2). This is a per-process budget in bytes for native cache workspace and CSR output, not a limit on total R memory. The peak during hash growth is included. A zero budget disables cache construction. If the budget or integer-index limit would be exceeded, the uncached calculation is used automatically. Each parallel worker may hold its own R objects.

The univariate mean convention is the same as in GeoKrig: \mu=X\beta with coefficients named mean, mean1, and so on, or an external observation mean vector in param$mean. The intercept-only case is equivalent to a one-column matrix of ones. At prediction tasks, use either Xloc for X_{loc}\beta or Mloc for a directly supplied local mean. The rows of Xloc and the entries of Mloc follow the local-task order; for fixed-location space-time kriging this is location first and prediction time second.

For type_krig="Universal", each local prediction uses the mean coefficients already supplied or estimated by GeoFit; they are not re-estimated by GLS inside each neighborhood. If mse=TRUE, the local MSE includes the uncertainty of those coefficients using the corresponding block of varcov. For composite likelihood this is the sandwich covariance (inverse Godambe information). Mean coefficients missing from varcov are treated as fixed. External mean vectors are already known, so a requested universal prediction is treated as simple kriging with a warning.

For model="Wrapped", each neighborhood delegates to the dedicated circular predictor in GeoKrig: sine and cosine residual components are predicted linearly with their exact wrapped-Gaussian covariances and the final direction is reconstructed with atan2. A requested type_krig="Universal" is therefore replaced by componentwise Simple prediction.

For copula="Gaussian", copula="Clayton", and copula="SkewGaussian", each local prediction uses the same observed-scale copula covariance, marginal-mean calculations, solver, and MSE implementation as GeoKrig. The local predictor is an optimal linear predictor, not in general the full conditional mean.

GeoKrigloc uses a single computational path. First, GeoNeighborhood constructs the requested neighborhood for each prediction task. Then GeoKrig is called on each local data set. This keeps covariance, copula, mean, numerical-solver, and MSE logic centralized in GeoKrig rather than duplicating those calculations inside GeoKrigloc.

With parallel=TRUE, independent local GeoKrig calls are submitted through future.apply. The default ncores=6 requests up to six workers, subject only to the number of prediction tasks and detected cores. Any other numeric ncores value is treated in the same way and is not reduced by the RAM safety heuristic. Set ncores=NULL explicitly to select the automatic/safe mode, where the package-wide worker resolver and the internal RAM safety check may reduce the worker count. On platforms where future reports that forked workers are safe, the automatic backend uses multicore; otherwise it uses multisession. Parallel and serial execution use the same local neighborhoods and the same GeoKrig calculations.

When geometric anisotropy is supplied, observation coordinates, prediction locations, and neighborhood selection are transformed by the same metric. Thus the selected local neighbors are the nearest points under the covariance metric, not under the original isotropic distance.

This function uses GeoKrig with a spatial or spatio-temporal neighborhood computed using GeoNeighborhood. Each local prediction task is evaluated through the same GeoKrig implementation used by the non-local prediction interface. The neighborhood is specified with neighb, maxdist, and maxtime. Regular-grid data are internally vectorized in the same order as expand.grid; fixed-location space-time observations remain in time-major order.

Value

Returns an object of class Kg. An object of class Kg is a list containing at most the following components:

bivariate	TRUE if spatial bivariate cokriging is performed, otherwise FALSE;
coordx	A d-dimensional vector of spatial coordinates used for prediction;
coordy	A d-dimensional vector of spatial coordinates used for prediction;
coordz	A d-dimensional vector of spatial coordinates used for prediction;
coordt	A t-dimensional vector of temporal coordinates used for prediction;
corrmodel	String: the correlation model;
covmatrix	The covariance matrix if type is Standard. An object of class spam if type is Tapering
data	The vector or matrix or array of data used for prediction
distance	String: the type of spatial distance;
grid	TRUE if the spatial data used for prediction are observed in a regular grid, otherwise FALSE;
loc	A (n \times 2)-matrix of spatial locations to be predicted.
n	The Binomial number of trials, or the common Negative-Binomial number r of successes.
nozero	In the case of tapered simple kriging the percentage of non zero values in the covariance matrix. Otherwise is NULL.
numcoord	Numeric:he number d of spatial coordinates used for prediction;
numloc	Numeric: the number n of spatial coordinates to be predicted;
numtime	Numeric: the number d of the temporal instants used for prediction;
numt	Numeric: the number m of the temporal instants to be predicted;
model	The response model used for prediction after canonicalizing any misspecified-Gaussian fitting name.
fit_model	The model name supplied by the caller or stored in the input GeoFit object before prediction canonicalization.
param	The parameter list used for prediction;
pred	For spatio-temporal prediction, a m \times n matrix with rows corresponding to time and columns corresponding to rows of loc; for spatial prediction, a numeric vector.
radius	Numeric: the radius of the sphere if coordinates are pssed in lon/lat format;
spacetime	TRUE if spatio-temporal kriging and FALSE if spatial kriging;
tapmod	String: the taper model if type is Tapering. Otherwise is NULL.
time	A m-dimensional vector of temporal coordinates to be predicted;
type	String: the type of kriging (Standard or Tapering).
type_krig	String: the type of kriging: Simple or Universal
mse	When mse=TRUE, the prediction MSE; otherwise NULL. For model="Wrapped", this is the sum of the sine- and cosine-component MSEs.
wrapped_mse_sin	For model="Wrapped" and mse=TRUE, the local sine-component MSE. This component is present only for model="Wrapped".
wrapped_mse_cos	For model="Wrapped" and mse=TRUE, the local cosine-component MSE. This component is present only for model="Wrapped".
wrapped_resultant	For model="Wrapped", the norm of the predicted first trigonometric-moment vector, truncated at one. This component is present only for model="Wrapped".
wrapped_prediction_method	For model="Wrapped", a short description of the componentwise circular prediction method. This component is present only for model="Wrapped".

Spatio-temporal ordering

Observed fixed-location data use time-major order: a T \times N matrix is vectorized as c(t(data)). Dynamic observations are supplied as aligned lists data[[t]] and coordx_dyn[[t]], concatenated by time. The rows of X and any known observation mean use the same observation order.

Prediction tasks use the different, location-major order

loc[1, ] at time[1], ..., loc[1, ] at time[Tloc],
loc[2, ] at time[1], ..., loc[2, ] at time[Tloc], ...

Rows of Xloc and elements of Mloc must follow this order. The returned pred and mse objects are Tloc \times Nloc matrices with prediction times in rows and prediction locations in columns. See GeoModels-spacetime-ordering.

Author(s)

Moreno Bevilacqua, [e-mail],https://sites.google.com/view/moreno-bevilacqua/home, Víctor Morales Oñate, [e-mail], https://sites.google.com/site/moralesonatevictor/, Christian Caamaño-Carrillo, [e-mail],https://www.researchgate.net/profile/Christian-Caamano

References

Gaetan, C. and Guyon, X. (2010) Spatial Statistics and Modelling. Springer-Verlag, New York. Furrer R., Genton, M.G. and Nychka D. (2006). Covariance Tapering for Interpolation of Large Spatial Datasets. Journal of Computational and Graphical Statistics, 15-3, 502–523.

See Also

GeoCovmatrix

Examples

################################################################
############### Examples of Spatial local kriging #############
################################################################
require(GeoModels)
####
model="Gaussian"

# Define the spatial-coordinates of the points:
set.seed(759)
x = runif(1000, 0, 1)
y = runif(1000, 0, 1)
coords=cbind(x,y)
# Set the exponential cov parameters:
corrmodel = "GenWend"
mean=0; sill=1
nugget=0; scale=0.2
param=list(mean=mean,sill=sill,nugget=nugget,smooth=0,
scale=scale,power2=4)

# Simulation of the spatial Gaussian random field:
data = GeoSim(coordx=coords, corrmodel=corrmodel,
 param=param)$data

# Maximum pairwise likelihood fitting of the space time random field:

start=list(scale=scale,sill=sill,mean=mean)
fixed=list(power2=4,smooth=0,nugget=0)
fit = GeoFit(data, coordx=coords, corrmodel=corrmodel, 
 start=start,fixed=fixed,
 likelihood='Conditional', type='Pairwise',
 neighb=3)

# locations to predict
loc_to_pred=matrix(runif(8),4,2)
################################################################
###
### Example 1. Comparing spatial kriging with local kriging for
### a Gaussian random field with GenWend correlation.
### 
###############################################################
param=append(fit$param,fit$fixed)
pr=GeoKrig(fit,loc=loc_to_pred,mse=TRUE)

pr_loc=GeoKrigloc(fit,loc=loc_to_pred,neighb=100,mse=TRUE)

pr$pred;
pr_loc$pred


############################################################
#### Example: spatio temporal Gaussian local kriging ######
############################################################


require(GeoModels)
set.seed(78)
coords=cbind(runif(100),runif(100))
coordt=seq(0,5,0.25)
corrmodel="Matern_Matern"
param=list(nugget=0,mean=0,scale_s=0.2/3,scale_t=0.25/3,sill=2,
 smooth_s=0.5,smooth_t=0.5)

data = GeoSim(coordx=coords, coordt=coordt,
 corrmodel=corrmodel, param=param)$data


# Maximum pairwise likelihood fitting of the space time random field:
start = list(scale_s=0.2/3,scale_t=0.25,sill=2,mean=0)
fixed = list(smooth_s=0.5,smooth_t=0.5,nugget=0)
I=Inf
lower=list(scale_s=0,scale_t=0,sill=0,mean=-I)
upper=list(scale_s=I,scale_t=I,sill=I,mean=I)
fit = GeoFit(data, coordx=coords, coordt=coordt, model=model, corrmodel=corrmodel, 
 likelihood='Conditional', type='Pairwise',start=start,fixed=fixed,
 optimizer="nlminb",lower=lower,upper=upper,
 neighb=3,maxtime=1)

## four location to predict
loc_to_pred=matrix(runif(8),4,2)
## three temporal instants to predict
time=c(0.5,1.5,3.5)


pr=GeoKrig(fit,loc=loc_to_pred,time=time,mse=TRUE)
pr_loc=GeoKrigloc(fit,loc=loc_to_pred,time=time,
 neigh=25,maxtime=1, mse=TRUE)

## full and local prediction 
pr$pred
pr_loc$pred



############################################################
#### Example: spatio bivariate Gaussian local cokriging ######
############################################################
#set.seed(6)
#NN=1500 # number of spatial locations 
#x = runif(NN, 0, 1); 
#y = runif(NN, 0, 1) 
#coords=cbind(x,y) 

## setting parameters
#mean_1 = 2; mean_2= -1
#nugget_1 =0;nugget_2=0
#sill_1 =0.5; sill_2 =1; 

### correlation parameters
#CorrParam("Bi_Matern")
#scale_1=0.2/3; scale_2=0.15/3; scale_12=0.5*(scale_2+scale_1) 
#smooth_1=smooth_2=smooth_12=0.5
#pcol = -0.4
#param= list(nugget_1=nugget_1,nugget_2=nugget_2,
# sill_1=sill_1,sill_2=sill_2,
# mean_1=mean_1,mean_2=mean_2,
# smooth_1=smooth_1, smooth_2=smooth_2,smooth_12=smooth_12,
# scale_1=scale_1, scale_2=scale_2,scale_12=scale_12,
# pcol=pcol)

## simulation
#data = GeoSim(coordx=coords, corrmodel="Bi_Matern",model=model,param=param)$data

#fixed=list(mean_1=mean_1,mean_2=mean_2, nugget_1=nugget_1,nugget_2=nugget_2, 
# smooth_1=smooth_1, smooth_2=smooth_2,smooth_12=smooth_12)

#start=list( sill_1=sill_1,sill_2=sill_2,
# scale_1=scale_1,scale_2=scale_2,scale_12=scale_12, pcol=pcol)

## estimation with maximum likelihood 
#fit = GeoFit(data=data,coordx=coords, corrmodel="Bi_Matern",
 # likelihood="Marginal",type="Pairwise",optimizer="BFGS",neighb=5,
 #start=start,fixed=fixed)

###### co-kriging for the fist component ##############
#xx=seq(0,1,0.022)
#loc_to_pred=as.matrix(expand.grid(xx,xx))
#pr1 = GeoKrigloc(fit,which=1,mse=TRUE,loc=loc_to_pred,neighb=100)

#opar=par(no.readonly = TRUE)
#par(mfrow=c(1,2))
#zlim=c(-2.5,2.5)
#colour = rainbow(100)
#fields::quilt.plot(coords,data[1,] ,col=colour,main = paste(" Fist component")) 
#fields::quilt.plot(loc_to_pred,pr1$pred,col=colour,
# main = paste(" Kriging first component"),ylab="")
#par(opar)

GeoModels documentation built on Sept. 23, 2026, 5:07 p.m.