Kriging_Ordinary: Ordinary Kriging interpolation
In InterpolateR: A Comprehensive Toolkit for Fast and Efficient Spatial Interpolation
View source: R/Kriging_Ordinary.R
Kriging_Ordinary R Documentation
Ordinary Kriging interpolation
Description
This function performs Ordinary Kriging, which is a spatial interpolation method that provides the best linear unbiased estimator
for unknown values based on the spatial correlation structure of the observed data. Unlike IDW, Kriging incorporates the spatial
autocorrelation through variogram modeling to determine optimal weights.
Usage
Kriging_Ordinary(
BD_Obs,
BD_Coord,
shapefile,
grid_resolution,
variogram_model = c("exponential", "spherical", "gaussian", "linear"),
max_dist = NULL,
n_lags = 15,
min_stations = 2,
n_round = NULL,
training = 1,
stat_validation = NULL,
Rain_threshold = NULL,
save_model = FALSE,
name_save = NULL
)
Arguments
BD_Obs
A data.table or data.frame containing observational data with the following structure:
The first column (Date): A Date object representing the observation date.
The remaining columns: Each column corresponds to a unique ground station, where the column name is the station identifier.
The dataset should be structured as follows:
> BD_Obs
# A data.table or data.frame with n rows (dates) and m+1 columns (stations + Date)
Date ST001 ST002 ST003 ST004 ...
<date> <dbl> <dbl> <dbl> <dbl> ...
1 2015-01-01 0 0 0 0 ...
2 2015-01-02 0 0 0 0.2 ...
3 2015-01-03 0.1 0 0 0.1 ...
Each station column contains numeric values representing observed measurements.
The column names (station identifiers) must be unique and match those in BD_Coord$Cod to ensure proper spatial referencing.
BD_Coord
A data.table or data.frame containing the metadata of the ground stations. It must include the following columns:
-
"Cod":
Unique identifier for each ground station.
-
"X":
Latitude of the station in UTM format.
-
"Y":
Longitude of the station in UTM format.
shapefile
A shapefile defining the study area, used to constrain the interpolation to the region of interest.
The shapefile must be of class SpatVector (from the terra package) and should have a UTM coordinate reference system.
grid_resolution
A numeric value indicating the resolution of the interpolation grid in kilometers (km).
variogram_model
Character string specifying the variogram model to fit. Options are:
-
"exponential" (default): Exponential variogram model
-
"spherical": Spherical variogram model
-
"gaussian": Gaussian variogram model
-
"linear": Linear variogram model
max_dist
Numeric value specifying the maximum distance (in meters) for variogram calculation and kriging prediction.
If NULL (default), it will be set to half the maximum distance between stations.
n_lags
Integer specifying the number of lag bins to use in the empirical variogram calculation. Default is 15.
min_stations
Integer specifying the minimum number of stations required for kriging interpolation.
If fewer stations are available (due to NAs), the function will return a constant field or use fallback methods.
Default is 2, which is more permissive for sparse precipitation data.
n_round
An integer specifying the number of decimal places to round the interpolated results.
If set to NULL, all decimal places will be preserved. The default value is 1.
training
Numerical value between 0 and 1 indicating the proportion of data used for model training. The remaining data are used for validation. Note that if you enter, for example, 0.8 it means that 80 % of the data will be used for training and 20 % for validation.
If you do not want to perform validation, set training = 1. (Default training = 1).
stat_validation
A character vector specifying the names of the stations to be used for validation.
This option should only be filled in when it is desired to manually enter the stations used for validation. If this parameter is NULL, and the formation is different from 1, a validation will be performed using random stations.
The vector must contain the names of the stations selected by the user for validation.
For example, stat_validation = c("ST001", "ST002"). (Default stat_validation = NULL).
Rain_threshold
List of numerical vectors defining precipitation thresholds to classify precipitation into different categories according to its intensity.
This parameter should be entered only when the validation is to include categorical metrics such as Critical Success Index (CSI), Probability of Detection (POD), False Alarm Rate (FAR), etc.
Each list item should represent a category, with the category name as the list item name and a numeric vector specifying the lower and upper bounds of that category.
Note: See the "Notes" section for additional details on how to define categories, use this parameter for validation, and example configurations.
save_model
Logical value indicating whether the interpolation file should be saved to disk. The default value is FALSE. indicating that the interpolated file should not be saved.
If set to TRUE, be sure to set the working directory beforehand using setwd(path) to specify where the files should be saved.
name_save
Character string indicating the name under which the interpolation raster file will be saved. By default the algorithm sets as output name: 'Model_Kriging'.
Value
The return value will depend on whether validation has been performed or not.
If validation is not performed, the function will return a SpatRaster object with the interpolated values.
If validation is performed, the function will return a list with two elements:
-
Ensamble: A SpatRaster object with the interpolated values.
-
Validation: A data.table with the validation results, including goodness-of-fit metrics and categorical metrics (if Rain_threshold is provided).
Details
Ordinary Kriging is a geostatistical interpolation technique that provides the best linear unbiased estimator (BLUE)
for spatial data. Unlike deterministic methods like IDW, Kriging incorporates the spatial correlation structure of the
data through variogram modeling.
The Ordinary Kriging estimator is defined as:
\hat{Z}(s_0) = \sum_{i=1}^{n} \lambda_i Z(s_i)
where:
\hat{Z}(s_0) is the estimated value at the unknown point.
Z(s_i) are the known values.
\lambda_i are the Kriging weights that sum to 1.
n is the total number of known points used.
The Kriging weights are obtained by solving the Kriging system:
\begin{bmatrix} \Gamma & \mathbf{1} \\ \mathbf{1}^T & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda} \\ \mu \end{bmatrix} = \begin{bmatrix} \boldsymbol{\gamma} \\ 1 \end{bmatrix}
where:
\Gamma is the variogram matrix between observation points.
\boldsymbol{\gamma} is the variogram vector between prediction and observation points.
\mu is the Lagrange multiplier ensuring unbiasedness.
The variogram models available are:
- Exponential
-
\gamma(h) = \sigma^2(1 - \exp(-3h/a)) for h > 0
- Spherical
-
\gamma(h) = \sigma^2(1.5h/a - 0.5(h/a)^3) for h \leq a, \sigma^2 for h > a
- Gaussian
-
\gamma(h) = \sigma^2(1 - \exp(-3(h/a)^2)) for h > 0
- Linear
-
\gamma(h) = \sigma^2 \cdot h/a for h \leq a, \sigma^2 for h > a
The Rain_threshold parameter is used to calculate categorical metrics such as the Critical Success Index (CSI),
Probability of Detection (POD), False Alarm Rate (FAR), success ratio (SR), Hit BIAS (HB),Heidke Skill Score (HSS);
Hanssen-Kuipers Discriminant (HK); Equal Threat Score (ETS) or Gilbert Skill Score.
The parameter should be entered as a named list, where each item represents a category and the name of the item is the category name.
The elements of each category must be a numeric vector with two values: the lower and upper limits of the category.
For example:
Rain_threshold = list(
no_rain = c(0, 1),
light_rain = c(1, 5),
moderate_rain = c(5, 20),
heavy_rain = c(20, 40),
violent_rain = c(40, Inf)
)
Author(s)
Marco Mogro marcov.mogro@ucuenca.edu.ec
References
Matheron, G. (1963). Principles of geostatistics. Economic Geology, 58(8), 1246-1266.
Cressie, N. (1993). Statistics for Spatial Data. John Wiley & Sons.
Examples
# Load data from on-site observations
data("BD_Obs", package = "InterpolateR")
data("BD_Coord", package = "InterpolateR")
# Load the study area where the interpolation is performed.
shp_path = system.file("extdata", "study_area.shp", package = "InterpolateR")
shapefile = terra::vect(shp_path)
# Perform the interpolation
Interpolated_data <- Kriging_Ordinary(BD_Obs, BD_Coord, shapefile,
grid_resolution = 5, variogram_model = "exponential",
max_dist = NULL, n_lags = 15, n_round = 1, training = 0.8,
Rain_threshold = NULL, stat_validation = NULL,
save_model = FALSE, name_save = NULL)
InterpolateR documentation built on Sept. 9, 2025, 5:59 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
View source: R/Kriging_Ordinary.R
| Kriging_Ordinary | R Documentation |
Ordinary Kriging interpolation
Description
This function performs Ordinary Kriging, which is a spatial interpolation method that provides the best linear unbiased estimator for unknown values based on the spatial correlation structure of the observed data. Unlike IDW, Kriging incorporates the spatial autocorrelation through variogram modeling to determine optimal weights.
Usage
Kriging_Ordinary(
BD_Obs,
BD_Coord,
shapefile,
grid_resolution,
variogram_model = c("exponential", "spherical", "gaussian", "linear"),
max_dist = NULL,
n_lags = 15,
min_stations = 2,
n_round = NULL,
training = 1,
stat_validation = NULL,
Rain_threshold = NULL,
save_model = FALSE,
name_save = NULL
)
Arguments
BD_Obs |
A
The dataset should be structured as follows: > BD_Obs # A data.table or data.frame with n rows (dates) and m+1 columns (stations + Date) Date ST001 ST002 ST003 ST004 ... <date> <dbl> <dbl> <dbl> <dbl> ... 1 2015-01-01 0 0 0 0 ... 2 2015-01-02 0 0 0 0.2 ... 3 2015-01-03 0.1 0 0 0.1 ...
|
BD_Coord |
A
|
shapefile |
A shapefile defining the study area, used to constrain the interpolation to the region of interest.
The shapefile must be of class |
grid_resolution |
A numeric value indicating the resolution of the interpolation grid in kilometers ( |
variogram_model |
Character string specifying the variogram model to fit. Options are:
|
max_dist |
Numeric value specifying the maximum distance (in meters) for variogram calculation and kriging prediction.
If |
n_lags |
Integer specifying the number of lag bins to use in the empirical variogram calculation. Default is 15. |
min_stations |
Integer specifying the minimum number of stations required for kriging interpolation. If fewer stations are available (due to NAs), the function will return a constant field or use fallback methods. Default is 2, which is more permissive for sparse precipitation data. |
n_round |
An integer specifying the number of decimal places to round the interpolated results.
If set to |
training |
Numerical value between 0 and 1 indicating the proportion of data used for model training. The remaining data are used for validation. Note that if you enter, for example, 0.8 it means that 80 % of the data will be used for training and 20 % for validation. If you do not want to perform validation, set training = 1. (Default training = 1). |
stat_validation |
A character vector specifying the names of the stations to be used for validation. This option should only be filled in when it is desired to manually enter the stations used for validation. If this parameter is NULL, and the formation is different from 1, a validation will be performed using random stations. The vector must contain the names of the stations selected by the user for validation. For example, stat_validation = c("ST001", "ST002"). (Default stat_validation = NULL). |
Rain_threshold |
List of numerical vectors defining precipitation thresholds to classify precipitation into different categories according to its intensity. This parameter should be entered only when the validation is to include categorical metrics such as Critical Success Index (CSI), Probability of Detection (POD), False Alarm Rate (FAR), etc. Each list item should represent a category, with the category name as the list item name and a numeric vector specifying the lower and upper bounds of that category. Note: See the "Notes" section for additional details on how to define categories, use this parameter for validation, and example configurations. |
save_model |
Logical value indicating whether the interpolation file should be saved to disk. The default value is |
name_save |
Character string indicating the name under which the interpolation raster file will be saved. By default the algorithm sets as output name: 'Model_Kriging'. |
Value
The return value will depend on whether validation has been performed or not.
If validation is not performed, the function will return a SpatRaster object with the interpolated values.
If validation is performed, the function will return a list with two elements:
-
Ensamble: ASpatRasterobject with the interpolated values. -
Validation: Adata.tablewith the validation results, including goodness-of-fit metrics and categorical metrics (ifRain_thresholdis provided).
Details
Ordinary Kriging is a geostatistical interpolation technique that provides the best linear unbiased estimator (BLUE) for spatial data. Unlike deterministic methods like IDW, Kriging incorporates the spatial correlation structure of the data through variogram modeling.
The Ordinary Kriging estimator is defined as:
\hat{Z}(s_0) = \sum_{i=1}^{n} \lambda_i Z(s_i)
where:
\hat{Z}(s_0)is the estimated value at the unknown point.
Z(s_i)are the known values.
\lambda_iare the Kriging weights that sum to 1.
nis the total number of known points used.
The Kriging weights are obtained by solving the Kriging system:
\begin{bmatrix} \Gamma & \mathbf{1} \\ \mathbf{1}^T & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda} \\ \mu \end{bmatrix} = \begin{bmatrix} \boldsymbol{\gamma} \\ 1 \end{bmatrix}
where:
\Gammais the variogram matrix between observation points.
\boldsymbol{\gamma}is the variogram vector between prediction and observation points.
\muis the Lagrange multiplier ensuring unbiasedness.
The variogram models available are:
- Exponential
-
\gamma(h) = \sigma^2(1 - \exp(-3h/a))forh > 0 - Spherical
-
\gamma(h) = \sigma^2(1.5h/a - 0.5(h/a)^3)forh \leq a,\sigma^2forh > a - Gaussian
-
\gamma(h) = \sigma^2(1 - \exp(-3(h/a)^2))forh > 0 - Linear
-
\gamma(h) = \sigma^2 \cdot h/aforh \leq a,\sigma^2forh > a
The Rain_threshold parameter is used to calculate categorical metrics such as the Critical Success Index (CSI),
Probability of Detection (POD), False Alarm Rate (FAR), success ratio (SR), Hit BIAS (HB),Heidke Skill Score (HSS);
Hanssen-Kuipers Discriminant (HK); Equal Threat Score (ETS) or Gilbert Skill Score.
The parameter should be entered as a named list, where each item represents a category and the name of the item is the category name.
The elements of each category must be a numeric vector with two values: the lower and upper limits of the category.
For example:
Rain_threshold = list(
no_rain = c(0, 1),
light_rain = c(1, 5),
moderate_rain = c(5, 20),
heavy_rain = c(20, 40),
violent_rain = c(40, Inf)
)
Author(s)
Marco Mogro marcov.mogro@ucuenca.edu.ec
References
Matheron, G. (1963). Principles of geostatistics. Economic Geology, 58(8), 1246-1266. Cressie, N. (1993). Statistics for Spatial Data. John Wiley & Sons.
Examples
# Load data from on-site observations
data("BD_Obs", package = "InterpolateR")
data("BD_Coord", package = "InterpolateR")
# Load the study area where the interpolation is performed.
shp_path = system.file("extdata", "study_area.shp", package = "InterpolateR")
shapefile = terra::vect(shp_path)
# Perform the interpolation
Interpolated_data <- Kriging_Ordinary(BD_Obs, BD_Coord, shapefile,
grid_resolution = 5, variogram_model = "exponential",
max_dist = NULL, n_lags = 15, n_round = 1, training = 0.8,
Rain_threshold = NULL, stat_validation = NULL,
save_model = FALSE, name_save = NULL)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.