NR 512 – Spatial Statistical Modeling
Exercise #7 – Modeling Geostatistical Data
Variogram Models and Ordinary Kriging
As discussed in class, a geostatistical dataset contains a series of sample points that characterizes some continuously varying phenomena occurring in a defined study region. These samples could represent anything that varies continuously and has a measurable location including soil samples, mineral deposits, vegetation cover, temperature, etc. The primary goal of geostatistics is to leverage spatial structure (as determined by the sample points) to infer the values of the continuously varying phenomena in unsampled areas (i.e. across the entire area of interest).
Geostatistical datasets have first-order properties, which are related to the mean value of the attribute across the study region, and second-order properties, which are related to the spatial dependence (i.e. spatial arrangement) of the continuously varying phenomena across the study area.
Objectives: The goal of today’s lab is to learn how to fit variogram models and apply these models in the form of kriging to interpolate the phenomena at unsampled locations.
The R Package: Today we will be using R packages geoR and gstat, which provide functions for geostatistical data analysis using the software R.
The Datasets: Today you be working with a dataset called mesue, which contains locations and topsoil heavy metal concentrations, along with a number of soil and landscape variables at the observation locations, collected in a flood plain of the river Meuse, near the village of Stein (NL). Heavy metal concentrations are from composite samples of an area of approximately 15 m x 15 m.
Exercise 1 – Fitting Variogram Models
In geostatistics the theoretical variogram is a function describing the degree of spatial dependence between samples of a spatial process. It is defined as the variance of the difference between field values at two locations across a study region.
In the last lab we explored geostatistical data by plotting empirical variograms. As discussed in class, in order to model geostatistical data via kriging we first need to fit an appropriate model to the empirical variogram. This process is called variogram fitting.
Today’s lab will be conducted using various R packages.
- Start R
- If you have not already, download the geoR and gstat packages and required dependencies from CRAN (ask the instructor for assistance).
- Load the geoR and gstat packages using:
> require(geoR)
> require(gstat)
- We may also need additional packages for today’s lab. Type the following to access them. Note if you get an error, you may need to download and install the packages.
> require(lattice)
> require(RColorBrewer)
> require(sp)
- We concluded the last lab by plotting an empirical directional variograms to the meuse data. As you will recall, we plotted directional variograms because the meuse data were clearly anisotropic.
Let’s create and plot the directional variograms for the zinc information along the following azimuth directions (0, 45, 90, 135) again.
> data(meuse)
> data(meuse.grid)
> vgm.aniso <- variogram(log(zinc)~1, ~x+y, meuse, alpha = c(0, 45, 90, 135),
+ cutoff = 2000, width = 100)
> vgm.aniso
> plot(vgm.aniso)
What are the cutoff and width parameters doing in the variogram command?
______________________________________________________________________________
- Typically, when there is a trend, you will find the sill appears to be different in different directions. This is known as “zonal anisotropy”. There are a few ways to deal with this, 1) Fit a nested semivariogram model with an infinite range in one direction, 2) Fit a semivariogram just along a direction with no trend, or 3) Fit a trend to the data, and then fit a semivariogram with a sill to the residual (however, optimal fitting of a trend requires a semivariogram).
We will start with 2). For the most part, the data appear to be stationary in the 45 degree direction, and the trend appears to be largest in the 135 degree direction. Let’s plot the variogram in the 45 degree direction. To create a new object with just the corresponding values for the 45 degree direction, type the following:
> dir.45 <- vgm.aniso$dir.hor == 45
> plot(vgm.aniso[dir.45,])
Based on the last 45 degree directional variogram plot visually estimate the following:
Sill: __________, Range: __________, Nugget: ___________
- Variogram models can be fit automatically and manually. We will experiment with both. Regardless of the methods, we have to decide which type model form (e.g., spherical, Gaussian, exponential, etc.) to fit.
When manually fitting to the empirical variogram we need to guess which model might fit
the data best. Based on plot you just made, which model seems most appropriate. Think if the data seems to have a convex or concave shape for the bins shorter than the range: ____________
- To fit the variogram manually we will need to provide the model form and estimates of the sill, range, and nugget to the variogram model fitting command, vgm.
The order of parameters to “vgm” is partial sill, model form, range, and nugget, where the partial sill is the sill minus nugget. Plug in your visually estimated values to the next line of code separated by commas, when doing so, put the model in as a spherical model which the command will want as “Sph”, including the quotations.
> eye.sph <- vgm()
> plot(vgm.aniso[dir.45,], eye.sph)
Expand the data frame that this created, what do you notice about the fit parameters compared to what you provided the vgm command?
_____________________________________________________________________________
_____________________________________________________________________________
- Overall, this probably looks pretty good. Now let’s fit the variogram automatically. This will fit spherical, exponential, matern, and gaussian models to the data. According to the online manual, the effective ranges for these models are a, 3a, sqrt(3a), and 4a for a Matern(kappa = 1) model, where a is the range.
> fit.sph <- fit.variogram(vgm.aniso[dir.45,], vgm(0.42, “Sph”, 1200, 0.05))
> fit.exp <- fit.variogram(vgm.aniso[dir.45,], vgm(0.42, “Exp”, 1200/3 , 0.05))
> fit.mat.025 <- fit.variogram(vgm.aniso[dir.45,],
+ vgm(0.42, “Mat”, 1200, 0.05, kappa = 0.25))
> fit.mat.100 <- fit.variogram(vgm.aniso[dir.45,],
+ vgm(0.42, “Mat”, 1200/4, 0.05, kappa = 1))
> fit.gau <- fit.variogram(vgm.aniso[dir.45,], vgm(0.42, “Gau”, 1200/sqrt(3), 0.05))
We can plot all of these for comparison using:
> plot(vgm.aniso[dir.45,], fit.sph)
> plot(vgm.aniso[dir.45,], fit.exp)
> plot(vgm.aniso[dir.45,], fit.mat.025)
> plot(vgm.aniso[dir.45,], fit.mat.100)
> plot(vgm.aniso[dir.45,], fit.gau)
Look through each of the model datasets created, what metric can we use to statistically evaluate which model best fits our data? _________________________________________________________
Which model does this point us to? ________________________________________________
- Either the other models just don’t fit the data, or the automatic procedure had trouble with them. If you have reasons to theoretically prefer one of the other semivariogram models, you would probably have to play with the parameters to find an acceptable fit.
Now, replot the series of directional empirical variograms so that we can use it to modify our best fitting model:
> plot(vgm.aniso)
Based on this, lets modify the best fitting spherical semivariogram model for anisotropy. By eyeballing the directionally sampled semivariograms, if the major direction is 45 degrees with a range of 1200, I will say that the minor axis (135 degrees) has a range of a little under 400. This is ~30% of 1200.
> eye.sph.aniso <- vgm(0.42, “Sph”, 1200, 0.05, anis = c(45, 0.3))
> plot(vgm.aniso, eye.sph.aniso)
We will now perform ordinary kriging based on the isotropic and anisotropic variogram models. To fit an isotropic model, use:
> x.iso <- krige(log(zinc)~1, ~x+y, meuse, meuse.grid, model = fit.sph)
> x2.iso <- x.iso
> gridded(x.iso) <- ~x+y
Challenge: Try rewriting this code to perform the ordinary kriging with the anisotropic variogram model and report out the min, mean, and max of the predicted values.
min = ____________; mean ____________; max ____________
We can plot the resulting predicted raster by using:
> spplot(x.iso[“var1.pred”], main = “ordinary kriging predictions”)
Print this map and turn in with your lab.
We can also plot the kriging variance by typing:
> spplot(x.iso[“var1.var”], main = “ordinary kriging variance”)
To make this a little more impactful, we can add the sample locations to the kriging variance map by combining a few commands:
> plot(x.iso[“var1.var”], main = “Isotropic Ordinary Kriging Variance”)
> points(meuse$x, meuse$y, pch = “+”, cex=0.5, col = “white”)
> box(which = “plot”, lty = “solid”)Print this map and turn in with your lab.
Print this map and turn it in with your lab.
What do you notice about the relationship between kriging variance and the locations of the of the sample locations? What does this tell you about the kriging variance function?
______________________________________________________________________________
______________________________________________________________________________
Challenge: Now try and create maps of the anisotropic kriging prediction and variance.
Print these maps and turn it in with your lab.
What do you notice by visually inspecting the difference between the isotropic and anisotropic maps?
_____________________________________________________________________________
_____________________________________________________________________________
- We can compare the two different kriging estimates by determining and plot the difference map in the maps by using the following code:
> x.diff <- x2.iso
> x.diff$var1.pred <- x2.iso$var1.pred – x2.aniso$var1.pred
> gridded(x.diff) <- ~x+y
> spplot(x.diff[“var1.pred”], panel = function(…) {
panel.gridplot(…)
panel.xyplot(meuse$x, meuse$y, pch=”+”)
},
main = “Difference between isotropy and anisotropy”)
- Comment on the differences you see: ___________________________________________
______________________________________________________________________________
______________________________________________________________________________
- Now fit a model with zonal anisotropy. Look again at the anisotropic model we already fit.
> plot(vgm.aniso, eye.sph.aniso)
I will say that I want the range in the 135 direction to be 800, and the sill to be 0.4 higher than it is now. I only want this model in the 135 degree direction, so I will make it disappear in the 45 degree direction. To do this, give it an effective range of infinity.
> eye.sph.aniso2 <- vgm(0.4, “Sph”, 800*10000, anis = c(45, 0.0001),
+ add.to = eye.sph.aniso)
> plot(vgm.aniso, eye.sph.aniso2)
Not exactly the best fit, but this is perhaps the best that can be done.
- Now krige with this new anisotropic plot.
> x.aniso2 <- krige(log(zinc)~1, ~x+y, meuse, meuse.grid, model = eye.sph.aniso2)
> x2.aniso2 <- x.aniso2
> gridded(x.aniso2) <- ~x+y
# Plot the kriging result using:
> spplot(x.aniso2[“var1.pred”], main = “ordinary kriging predictions”)
# Plot the kriging variance:
> spplot(x.aniso2[“var1.var”], main = “ordinary kriging variance”)
# Add sample points to the variance plot:
> spplot(x.aniso2[“var1.var”], panel = function(…) {
panel.gridplot(…)
panel.xyplot(meuse$x, meuse$y, pch = “+”)
},
main = “Ordinary kriging variance”)
# Compare the anisotropic predictions by creating a difference map.
> x.diff2 <- x2.aniso2
> x.diff2$var1.pred <- x2.aniso$var1.pred – x2.aniso2$var1.pred
> gridded(x.diff2) <- ~x + y
> spplot(x.diff2[“var1.pred”], panel = function(…) {
panel.gridplot(…)
panel.xyplot(meuse$x, meuse$y, pch = “+”)
},
main = “Differences between two anisotropy models”)
- Comment on the differences you see:
______________________________________________________________________________
______________________________________________________________________________
______________________________________________________________________________