GeoVariogram: Empirical semivariogram estimation

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

Empirical semivariogram estimation

Description

Computes an empirical estimate of the semivariogram for spatial, spatio-temporal, and bivariate random fields.

Usage

GeoVariogram(data, coordx, coordy=NULL, coordz=NULL, coordt=NULL,
 coordx_dyn=NULL, cloud=FALSE, distance="Eucl",
 grid=FALSE, maxdist=NULL, neighb=NULL,
 maxtime=NULL, numbins=NULL,
 radius=1, type='variogram', bivariate=FALSE,
 subsample=1, subsample_t=1, numbins_t=NULL, directed=FALSE)

Arguments

data	A numeric vector of length d (a single spatial realisation), or a d \times d matrix (a single realisation on a regular grid), or a t \times d matrix (a single spatio-temporal realisation), or a d \times d \times t array (a single spatio-temporal realisation on a regular grid). See GeoFit for details.
coordx	Spatial coordinates. Either a numeric vector giving the first coordinate, or a d \times 2 (or d \times 3) matrix of coordinates. If distance refers to great-circle distances, coordinates must be provided in lon/lat format (decimal degrees) and the sphere radius is set by radius.
coordy	A numeric vector giving the second spatial coordinate. Optional, default is NULL.
coordz	A numeric vector giving the third spatial coordinate (if needed). Optional, default is NULL.
coordt	A numeric vector of temporal coordinates. If NULL (default), a purely spatial random field is assumed. Temporal coordinates may be irregularly spaced; temporal lags are computed from the supplied coordinate values.
coordx_dyn	For dynamic locations, a list with one two- or three-column coordinate matrix per temporal instant. The list length must match length(coordt); data[[t]] must have one value per row of coordx_dyn[[t]], in the same order.
cloud	Logical; if TRUE the semivariogram cloud is computed for a univariate purely spatial field. If FALSE (default), a binned empirical semivariogram is returned. For large data sets, combine cloud=TRUE with maxdist, neighb, or subsample to limit the number of returned pairs.
distance	String specifying the spatial distance. Default is "Eucl" (Euclidean distance). See GeoFit for details.
grid	Logical; if FALSE (default) data are interpreted as observations on irregularly spaced locations. If TRUE, data are interpreted as observations on a regular grid.
maxdist	A positive finite numeric maximum spatial distance. In the bivariate case a scalar or a length-3 vector (first marginal, cross, second marginal) is accepted. Use NULL to let the function use the available spatial range. See Details.
neighb	Numeric; an optional positive integer indicating the order of neighborhood (useful for large datasets). Neighborhood pair selection is available for spatial and spatio-temporal semivariograms. In the bivariate case a length-3 vector can be used for the first marginal, cross, and second marginal pair sets. See Details.
maxtime	A positive finite maximum temporal lag, expressed in the same units as coordt, to be considered for spatio-temporal semivariograms. See Details.
numbins	Historical GeoModels argument controlling the spatial bin grid. It is the number of spatial bin boundaries, so the number of empirical spatial classes is numbins - 1. The default is 13 boundaries (12 classes).
numbins_t	Optional positive integer giving the number of temporal classes for irregularly spaced coordt. If NULL, numbins - 1 classes are used. For regularly spaced times the attainable temporal lags are generated directly from the common spacing and this argument is not used.
directed	Logical used only when neighb is supplied. The default FALSE removes reciprocal nearest-neighbour duplicates so that the empirical semivariogram is based on unordered pairs. Set TRUE to retain the directed nearest-neighbour graph used by pairwise composite-likelihood code.
radius	Numeric; radius of the sphere when using great-circle distances. Default is 1.
type	String; type of semivariogram. Currently available: "variogram".
bivariate	Logical; if FALSE (default) data are interpreted as univariate spatial/spatio-temporal realisations. If TRUE, data is interpreted as a realisation from a bivariate field and (cross-)semivariograms are computed.
subsample	Numeric in (0,1]. Proportion of spatial locations to be used to compute the semivariogram (useful for large datasets). Default is 1 (use all locations).
subsample_t	Numeric in (0,1]. Proportion of time points to be used in spatio-temporal settings (when coordt is provided). Default is 1 (use all time points).

Details

We report the definition of the semivariogram in the spatial case; extensions to spatio-temporal and bivariate settings are based on the same principles.

For a spatial random field Z(\cdot), the (classical) binned semivariogram estimator is defined as

\hat{\gamma}(h) = \frac{1}{2 |N(h)|}\sum_{(x_i,x_j)\in N(h)} \{Z(x_i)-Z(x_j)\}^2,

where N(h) is the set of all sample pairs whose spatial distance falls within a tolerance region around lag h (equally spaced intervals are used when cloud=FALSE).

The historical numbins argument sets the number of spatial bin boundaries; hence numbins - 1 empirical spatial classes are formed when cloud=FALSE.

The maxdist argument sets a strictly positive finite maximum spatial distance. If no pair falls below the requested cutoff, a valid empty semivariogram object is returned instead of constructing decreasing bins.

The maxdist option can be combined with neighb to reduce the number of pairs when handling large datasets, by restricting computations to local neighborhoods. By default reciprocal directed nearest-neighbour edges are deduplicated before binning; directed=TRUE retains the directed graph. For chordal and geodesic distances, neighbour ordering is obtained from three-dimensional unit-sphere coordinates and the reported lags are then evaluated in the requested metric.

Spatial and temporal bins are left-closed and right-open, except for the last bin which is closed on the right. Thus a pair exactly at the maximum retained lag is not discarded.

The maxtime argument sets the maximum temporal lag considered for spatio-temporal semivariograms. For regularly spaced times, attainable temporal lags are generated in linear memory/time from the common spacing; no T \times T matrix of all time differences is formed. For irregular times a controlled grid of numbins_t temporal classes is used. The returned spatio-temporal surface is rectangular on centers by centert; cells with no valid pairs are returned as NA with count zero rather than being removed.

For dynamic sites, the temporal marginal is the empirical \gamma(0,u) and therefore uses only spatially collocated locations (up to numerical tolerance). Nearby but non-collocated locations contribute to the positive- distance space-time surface, not to the temporal marginal. Consequently, in a fully dynamic design with no spatial locations repeated across temporal instants, variogramt is NA in the affected temporal classes. This is an expected property of the sampling design, not a failure of the space-time variogram: the positive-distance surface \gamma(h,u) remains empirically estimable. In this case plot.GeoVariogram labels the panel “Temporal marginal unavailable (no repeated spatial locations)”. Repeated spatial locations across times are needed only when the empirical temporal marginal itself is required.

In the bivariate case the two marginal semivariograms and the cross-semivariogram are evaluated on the same spatial bins. The cross-semivariogram uses the classical increment-product estimator on unordered positive-distance pairs. The trivial collocated contribution at lag zero is not mixed into the first positive spatial class. If two dynamic/support coordinate sets are supplied, they must be aligned to define this estimator unambiguously.

The subsample and subsample_t arguments provide additional control for large datasets by using only a proportion of spatial locations and/or time points. With grid=TRUE, the grid is first normalized to explicit spatial locations, so full-grid and subsampled calculations follow the same path.

Missing values NA/NaN are allowed and pairs involving them are skipped consistently. Infinite observations are rejected.

Value

Returns an object of class "GeoVariogram". The list contains, as applicable:

bins	Spatial bin boundaries when cloud=FALSE; spatial pair distances when cloud=TRUE.
bint	Temporal lag representatives for a spatio-temporal variogram.
bivariate	Logical indicating a bivariate empirical variogram.
cloud	Logical indicating a variogram cloud.
centers	Spatial bin centers.
centert	Temporal lag representatives used by the rectangular spatio-temporal surface.
distance	Spatial distance type used to construct the empirical variogram.
radius	Sphere radius used for chordal or geodesic distances.
grid	Logical recording whether the original input was supplied as a grid.
neighb	Neighborhood order used for pair selection, or NULL.
directed	Whether reciprocal nearest-neighbour pairs were retained.
lenbins	Numbers of pairs in the spatial bins. In the bivariate case this is a two-row matrix.
lenbinst	Numbers of pairs in the cross/spatio-temporal bins. For space-time objects this follows the same row-major ordering as variogramst.
lenbint	Numbers of pairs in the temporal bins.
maxdist	Maximum spatial distance requested by the user.
maxtime	Maximum temporal lag requested by the user.
numbins_t	Requested number of temporal classes for irregular times, or NULL.
regular	Logical indicating regularly spaced temporal coordinates for a space-time object.
time.breaks	Internal temporal bin boundaries used for space-time binning.
spacetime_dyn	Logical indicating dynamic spatial coordinates.
temporal.margin	For space-time objects, a label indicating that the temporal margin is based on same-site/collocated pairs.
subsample	Spatial subsampling proportion.
subsample_t	Temporal subsampling proportion.
variograms	Empirical spatial semivariogram; a two-row matrix in the bivariate case.
variogramst	Empirical cross-semivariogram in the bivariate case, or the rectangular spatio-temporal surface stored in row-major order. Empty cells are NA.
variogramt	Empirical temporal marginal semivariogram.
type	Type of empirical variogram.

Spatio-temporal ordering

For fixed sites, data is a T \times N matrix with times in rows and sites in columns, and is internally read in the order c(t(data)). For dynamic sites, data and coordx_dyn are aligned lists; the function concatenates complete temporal blocks in list order. 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

Cressie, N. A. C. (1993) Statistics for Spatial Data. New York: Wiley.

Gaetan, C. and Guyon, X. (2010) Spatial Statistics and Modeling. Springer-Verlag, New York.

See Also

GeoFit for model fitting, GeoWLS for weighted least-squares estimation from empirical variograms, GeoCovariogram for fitted covariance and variogram values, plot.GeoVariogram for plotting empirical variograms.

Examples

library(GeoModels)

################################################################
### Example 1. Empirical semivariogram from a spatial Gaussian
### random field with Matérn correlation.
################################################################
set.seed(514)
x = runif(200, 0, 1)
y = runif(200, 0, 1)
coords = cbind(x,y)

corrmodel = "Matern"
mean = 0
sill = 1
nugget = 0
scale = 0.3/3
smooth = 0.5

data = GeoSim(coordx=coords, corrmodel=corrmodel,
 param=list(mean=mean, smooth=smooth, sill=sill,
 nugget=nugget, scale=scale))$data

vario = GeoVariogram(coordx=coords, data=data, maxdist=0.6)
plot(vario, pch=20, ylim=c(0,1), ylab="Semivariogram", xlab="Distance")

################################################################
### Example 2. Empirical semivariogram for a spatio-temporal
### Gaussian random field with Gneiting correlation.
################################################################
set.seed(331)
x = runif(200, 0, 1)
y = runif(200, 0, 1)
coords = cbind(x,y)
times = seq(1,10,1)

data = GeoSim(coordx=coords, coordt=times, corrmodel="gneiting",
 param=list(mean=0, scale_s=0.08, scale_t=0.4, sill=1,
 nugget=0, power_s=1, power_t=1, sep=0.5))$data

vario_st = GeoVariogram(data=data, coordx=coords, coordt=times,
 maxtime=7, maxdist=0.5)
plot(vario_st, pch=20)

################################################################
### Example 3. Empirical (cross-)semivariograms for a bivariate
### Gaussian random field with Bi-Matérn covariance.
################################################################
set.seed(293)
x = runif(400, 0, 1)
y = runif(400, 0, 1)
coords = cbind(x,y)

param = list(mean_1=0, mean_2=0,
 scale_1=0.1/3, scale_2=0.15/3, scale_12=0.15/3,
 sill_1=1, sill_2=1,
 nugget_1=0, nugget_2=0,
 smooth_1=0.5, smooth_12=0.5, smooth_2=0.5,
 pcol=0.3)

data = GeoSim(coordx=coords, corrmodel="Bi_matern", param=param)$data
biv_vario = GeoVariogram(data, coordx=coords, bivariate=TRUE, maxdist=0.5)
plot(biv_vario, pch=20)

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