# Kriging using multiple spatial input variables and one spatial output (response) variable

I have a data table with the following columns in a CSV file.

``````X    Y     Zinc      L    F     C
---  ---   -----     --   ---  ----
``````

Of course, I can import the CSV file as a data frame

``````a <- read.csv('givendata.csv', header = T)
``````

And convert this data into spatial data by

``````coordinates(a) = ~X + Y
``````

Now, I wish to run kriging by using L, F and C as spatial input variables and Zinc as the output variable. How do I achieve this?

Edit: I am new to Kriging.

I am adding some more explanation to the question, based on a comment that I do not have an idea about linear modelling.

If I do the following

``````A1 <- read.csv('givendata.csv', header = T)
mymodel <- lm(Zinc ~ L + F + C, data = A1)
``````

it gives the output, with all the regression coefficients and the constant. But, these observations have spatial inputs (see the meuse dataset), which is why it is referred to as a spatial variable (this is GIS.SE!).

Now suppose I were to take the following steps:

``````x.range <- (range(A1@coords[,1]))
y.range <- (range(A1@coords[,2]))

A1.grd <- expand.grid(x=seq(from=x.range, to=x.range, by=30),
y=seq(from=y.range, to=y.range, by=30))

coordinates(A1.grd) <- ~x+y
gridded(A1.grd) <- TRUE

library(automap)
kZinc <- autoKrige(Zinc ~ L + F + C, A1, A1.grd)
``````

I get the following error:

``````Error in eval(expr, envir, enclos) : object 'L' not found
``````

which is funny because the `lm` was working and the Kriging does not work!

• Are L,F, and C explanatory variables? If this was a non-spatial model you'd just do `m = glm(Zinc ~ L+F+C)` yes? – Spacedman May 5 '17 at 19:45
• `Zinc ~ L + F + C` does not appear to work with kriging. – Indian May 6 '17 at 9:43
• I didn't say it did. I just asked if L, F, and C were explanatory variables. You call them "spatial input variables" which is not a recognised term. In linear modelling terms, is "Zinc" your response, and L, F, and C, your explanatory variables, or covariates? If you don't understand those terms you need to understand linear modelling before you attempt kriging. – Spacedman May 6 '17 at 14:09
• Your claim about me not understanding response, iid or covariates is blatant and non-chalant. If we could get back to how R addresses the issue, moving away the attention from me, it would be great. Thanks! – Indian May 7 '17 at 11:14
• Added sufficient explanation to the question, so that any doubt does not arise about my understanding. – Indian May 7 '17 at 11:36

As I stated in the comments, you need values of your explanatory variables at the prediction locations in order to make predictions. Here is a fully reproducible example:

``````library(automap)
``````

Make 100 random data points since we don't have your data file:

``````A1 = data.frame(x=runif(100),y=runif(100), Zinc=runif(100), L=runif(100), F=runif(100), C=runif(100))
coordinates(A1)=~x+y
``````

Predict over this grid (note I've used `len=20` instead of `by=30` from your code because my numbers are different):

``````x.range <- (range(A1@coords[,1]))
y.range <- (range(A1@coords[,2]))

A1.grd <- expand.grid(x=seq(from=x.range, to=x.range, len=20),
y=seq(from=y.range, to=y.range, len=20))
coordinates(A1.grd) <- ~x+y
gridded(A1.grd) <- TRUE
``````

Now let's reproduce what you did:

``````kZinc <- autoKrige(Zinc ~ L + F + C, A1, A1.grd)
``````

That reproduces your error. Now I can see what your problem is, something impossible to tell from your initial post. Thank you for editing.

The problem is the model wants to predict Zinc at the grid locations, but you need values of L, F, and C at those locations, in order to do predictions. Let's just make some up to illustrate. I have a 20x20 grid so that's 400 points:

``````A1.grd\$L=runif(400)
A1.grd\$F=runif(400)
A1.grd\$C=runif(400)
``````

Moment of truth:

kZinc <- autoKrige(Zinc ~ L + F + C, A1, A1.grd)

And that works.

Now we have no idea what your L, F, and C represent. Normally they will be one of two kinds of things:

1. Spatially continuous quantities that you can get the value of at any (x,y) location, such as altitude, or distance to nearest road. Compute them over your (x,y) grid and do kriging as above.

2. Non-continuous quantities only valid at the sample level. For example if these are the amounts of Zinc in human blood samples you might have the age of the human sampled. How do you then do kriging predictions at a point? You don't have sample-specific measurements there because you haven't any samples there.

In this case you can literally make up some values. To continue my example you might choose the average age of people in your study, do the kriging with that, and then show a kriged map of "Zinc amount for 27 year olds". With more variables you may end up with something like "Zinc amount for unvaccinated 27 year old men" if your variables are vaccination status, age, and sex.

You could also do separate kriged maps for extreme values of your covariates, and show them on two plots, something like "Range of Zinc values for ages 12 (left figure) to 52 (right figure)".

Or you might find that your covariates aren't significant at all and just drop them.

• I am using Satellite Data. Therefore multiple bands (or should I say channels?) should be in the form of multiple columns. – Indian May 8 '17 at 8:21

There are two widely used packages to achieve Kriging interpolation and variogram model fitting in `R`, these are gstat and geoR. Depending on the package you use, you can build a native `gstat` object or a `geoR` object. Also, you can manage to do all your work using the known `SpatialPointsDataFrame` object.

Some general steps prior to achieve Kriging interpolation would be:

1. Representing Variation: describe the variation of your data by meanings of Frequency distribution (Histograms and Box-plots), Centre (mean, median and mode) and Dispersion (range, interquartile range, standard deviation, variance, coefficient of variation, skewness and kurtosis) analysis.

For Frequency distribution analysis you can use `histogram` and `boxplot` functions or if you want more attractive plots you can use `ggplot2` library (`geom_histogram` and `geom_boxplot` functions). For Centre analysis you can use `mean`, `median` and this `Mode` function: `Mode <- function(x) {ux = unique(x); ux[which.max(tabulate(match(x, ux)))]}`. For Dispersion analysis you can use `range`, `IQR`, `sd` and `var` functions, calculate coefficient of variation as `(sd/mean)*100` and `skewness` and `kurtosis` functions from `moments` library. Also, try to find outliers at this points and evaluate if it's better to remove them or not.

1. Transformation: attempt to transform the measured values to a new scale on which the distribution is more nearly normal to overcome difficulties arising from departures from normality.

To test normality of your data you can perform a Shapiro-Wilk test with `shapiro.test` function. If your data is not normal you can try simple transformations (`log` and `sqrt`) or see Box-Cox transformations using `forecast` library (`BoxCox.lambda` and `BoxCox` functions). Also, you can perform a "visual" analysis of normality using Q-Q plots, you can use `ggQQ <- ggData + geom_qq(aes(sample = var)) + geom_abline(intercept = mean(var), slope = sd(var), col = "red")` from `ggplot2` library.

1. Trend analysis: discover trends between the response variable and co-variables using scatter plots with smoothed lines and/or correlations and/or correlograms.

To make a scatter plot you just can simply `plot(x = var1, y = var2)` and try a linear model `abline(lm(var2 ~ var1, data = data), col = "red")` or using `ggplot2` library `ggTrend_var1var2 <- ggData + geom_point(aes(y = var1, x = var2)) + geom_smooth(aes(y = var1, x = var2))`. For correlation you can use Pearson correlation test: `cor.test(x = var2, y = var1, method = "pearson")`. To perform correlograms you can use the function `corrgram` from library `corrgram`. After analysing possible trends by means of correlations and plots, you can evaluate some trend models e.g.: multiple linear regression model with vars 2, 3 and 4 (`model1 <- lm(formula = var1 ~ var2 + var3 + var4, data = data)`) or a simple linear regression model with only var 2 (`model2 <- lm(formula = var1 ~ var2, data = data)`) and after that choose the best model with Akaike criterion `AIC(model1, model2)`.

1. Find spatial outliers: try to find which points values are atypically different from their neighborhood point values, you can use Moran Index to quantify it and make a "posting" plot to explore your data. At this point you need to evaluate to remove or not (depending on your data) those spatial outliers points.

You can use `moran.plot` from `spdep` library to make a Moran Scatterplot and get influenced values. To perform a "posting" plot you can use the function `spplot(obj = SpatialPointsDataFrame, zcol = "var1", colorkey = TRUE)` from `sp` library or the great interactive map `mapview(SpatialPointsDataFrame, zcol = "var1")` from `mapview` library.

1. Empiric variogram: characterizing the spatial processes of your data.

Here you can use the function `variog` from `geoR` library or `variogram` function from `gstat` library. In `geoR` and `gstat` the trend can be included as a formula: `trend = trend.spatial(var1 ~ var2 + var3)` option in `geoR::variog` function and `object = var1 ~ var2 + var3` option in 'gstat::variogram'. At this point you should consider find possible anisotropy process in your data. To check that, you have to perform directional variograms (typically 0, 45, 90 and 135 degrees). See `direction`, `tolerance` and `unit.angle` options in `geoR::variog` function or `alpha` option in `gstat::variogram` function.

1. Adjust theoretical variogram: adjust a theoretical variogram model to the empirical variogram. The choose of the model will depend on the best fit of a series of theoretic variogram models tried or one choose based on your previously knowledge of how the spatial structure of your data is best described.

See `gstat::vgm()`, `geoR::likfit` and `geoR::variofit` functions.

1. Kriging interpolation: there are different types of Kriging interpolations, you should consider which one is best for you (or your data!). At this point, you need your variogram model (output from `vgm`, `likfit` or `variofit`) and know which points positions in space you are interested to predict (it can be a grid, or some irregular point data).

One example of Ordinary Kriging using `gstat` could be `gstat::krige(formula = var1 ~ var2 + var3, locations = data, newdata = pointsPositionsToPredict, model = variogramModel)`

• Note: `geoR` and `gstat` packages are not compatible, these libraries work with different objects. For example, you can't use the variogram object model from `geoR` in a `gstat::krige` function. Try use one package first, but you can use both packages building manually the objects using the parameters derived from one package function to the other.

• For further information you can check this book Geostatistics for Environmental Scientists, 2nd Edition (Webster & Oliver).

• A code work example IntroductionToGeostatistics (only spanish version)