Dear mailing list,
As I try to test representativeness of study sites through variograms I am
calculating several anisotropic variograms in R 2.14, from which I only
take a subset of point pairs. My desired output is a range value, which
tells me up to which distance there is a statistical dependance between
study site and neighboring landscape. These subsets all include the center
point of a polygon. I split the variograms by direction, the result is 8
distinct variograms, so that autofitVariogram can handle them (actually
autofitVariogram cannot handle anisotropic variograms. By splitting them I
avoid this). Thus, the directions of the variogram are: N, NE, E, SE, S,
SW, W, NW and N, always regarded to the center point. I fit each of them
automatically witha modification of autofitVariogram. Unfortunately, the
result is not satisfying, as one direction is completely different from the
others and it is in 52 of 52 cases the same direction which acts different.
Most of the other directions are reliable, which is even stranger than the
result itself. So I guess I have some wrong assumptions in the
autofitVariogram or variogram function.
The code is:
## needed functions
as_radians = function(theta = 0) {
return(theta * pi/180)
}
## modification of autofitVariogram, it is now possible to input already
calculated variograms
autofit = function (experimental_variogram, diagonal, boundaries,
model = c("Sph", "Exp", "Gau","Ste"),
kappa = c(0.05, seq(0.2, 2, 0.1), 5, 10),
fix.values = c(NA, NA, NA), verbose = FALSE, GLS.model
= NA,
start_vals = c(NA, NA, NA), miscFitOptions = list(),
...)
{
if ("alpha" %in% names(list(...)))
warning("Anisotropic variogram model fitting not supported, see the
documentation of autofitVariogram for more details.")
miscFitOptionsDefaults = list(merge.small.bins = TRUE, min.np.bin = 5)
miscFitOptions = modifyList(miscFitOptionsDefaults, miscFitOptions)
diagonal = diagonal
if (miscFitOptions[["merge.small.bins"]]) {
if (verbose)
cat("Checking if any bins have less than 5 points, merging bins
when necessary...\n\n")
while (TRUE) {
if (length(experimental_variogram$np[experimental_variogram$np
<
miscFitOptions[["min.np.bin"]]]) == 0 | length(boundaries)
==
1)
break
boundaries = boundaries[2:length(boundaries)]
}
}
if (is.na(start_vals[1])) {
initial_nugget = min(experimental_variogram$gamma)
}
else {
initial_nugget = start_vals[1]
}
if (is.na(start_vals[2])) {
initial_range = 0.1 * diagonal
}
else {
initial_range = start_vals[2]
}
if (is.na(start_vals[3])) {
initial_sill = mean(c(max(experimental_variogram$gamma),
median(experimental_variogram$gamma)))
}
else {
initial_sill = start_vals[3]
}
if (!is.na(fix.values[1])) {
fit_nugget = FALSE
initial_nugget = fix.values[1]
}
else fit_nugget = TRUE
if (!is.na(fix.values[2])) {
fit_range = FALSE
initial_range = fix.values[2]
}
else fit_range = TRUE
if (!is.na(fix.values[3])) {
fit_sill = FALSE
initial_sill = fix.values[3]
}
else fit_sill = TRUE
getModel = function(psill, model, range, kappa, nugget, fit_range,
fit_sill, fit_nugget, verbose) {
if (verbose)
debug.level = 1
else debug.level = 0
if (model == "Pow") {
warning("Using the power model is at your own risk, read the
docs of autofitVariogram for more details.")
if (is.na(start_vals[1]))
nugget = 0
if (is.na(start_vals[2]))
range = 1
if (is.na(start_vals[3]))
sill = 1
}
obj = try(fit.variogram(experimental_variogram, model = vgm(psill =
psill,
model = model, range = range, nugget = nugget, kappa = kappa),
fit.ranges = c(fit_range), fit.sills = c(fit_nugget,
fit_sill), debug.level = 0), TRUE)
if ("try-error" %in% class(obj)) {
warning("An error has occured during variogram fitting.
Used:\n",
"\tnugget:\t", nugget, "\n\tmodel:\t", model,
"\n\tpsill:\t", psill, "\n\trange:\t", range,
"\n\tkappa:\t", ifelse(kappa == 0, NA, kappa),
"\n as initial guess. This particular variogram fit is not
taken into account. \nGstat error:\n",
obj)
return(NULL)
}
else return(obj)
}
test_models = model
SSerr_list = c()
vgm_list = list()
counter = 1
for (m in test_models) {
if (m != "Mat" && m != "Ste") {
model_fit = getModel(initial_sill - initial_nugget,
m, initial_range, kappa = 0, initial_nugget,
fit_range, fit_sill, fit_nugget, verbose = verbose)
if (!is.null(model_fit)) {
vgm_list[[counter]] = model_fit
SSerr_list = c(SSerr_list, attr(model_fit, "SSErr"))
}
counter = counter + 1
}
else {
for (k in kappa) {
model_fit = getModel(initial_sill - initial_nugget,
m, initial_range, k, initial_nugget, fit_range,
fit_sill, fit_nugget, verbose = verbose)
if (!is.null(model_fit)) {
vgm_list[[counter]] = model_fit
SSerr_list = c(SSerr_list, attr(model_fit,
"SSErr"))
}
counter = counter + 1
}
}
}
strange_entries = sapply(vgm_list, function(v) any(c(v$psill,
v$range) < 0) | is.null(v))
if (any(strange_entries)) {
if (verbose) {
print(vgm_list[strange_entries])
cat("^^^ ABOVE MODELS WERE REMOVED ^^^\n\n")
}
warning("Some models where removed for being either NULL or having
a negative sill/range/nugget, \n\tset verbose == TRUE for more information")
SSerr_list = SSerr_list[!strange_entries]
vgm_list = vgm_list[!strange_entries]
}
if (verbose) {
cat("Selected:\n")
print(vgm_list[[which.min(SSerr_list)]])
cat("\nTested models, best first:\n")
tested = data.frame(`Tested models` = sapply(vgm_list,
function(x) as.character(x[2, 1])), kappa = sapply(vgm_list,
function(x) as.character(x[2, 4])), SSerror = SSerr_list)
tested = tested[order(tested$SSerror), ]
print(tested)
}
result = list(exp_var = experimental_variogram, var_model =
vgm_list[[which.min(SSerr_list)]],
sserr = min(SSerr_list))
class(result) = c("autofitVariogram", "list")
return(result)
}
####################################################################################
##needed packages
library(gstat)
library(automap)
library(plotrix)
## extract imaginary center point
y = meuse[34,]
## extract data, which are north and south of specific point
north = meuse[which(meuse[,3] >= y[,3]), ]
south = meuse[which(meuse[,3] <= y[,3]), ]
## convert to SpatialPointsDataFrame
coordinates(meuse) = ~x+y
coordinates(north) = ~x+y
coordinates(south) = ~x+y
## calculate variogram for east and west direction
v = variogram(zinc ~ 1, data = meuse, width = 90, cutoff = 4789.86784786386,
cloud = T, alpha = c(45, 90, 135), tol.hor = 90/3)
## extract point pairs where y is contained. If y is left part of point
pair (%%),
## pair direction is east. If y is right part of point pair (%/%), pair
direction
## is west
vwest <- v[((v$np%/% attr(v, ".BigInt"))+1)==as.numeric(row.names(y)), ]
veast <- v[((v$np%% attr(v, ".BigInt"))+1)==as.numeric(row.names(y)), ]
## variogram calculation for north direction
v <- variogram(zinc ~ 1, data = north, cutoff = 4789.86784786386, width =
90,
cloud = T, alpha = c(0), tol.hor = 90/3)
## extract point pairs, which include center point of polygon
vN <- v[((v$np%/% attr(v, ".BigInt"))+1)==as.numeric(row.names(y)) |
((v$np%% attr(v, ".BigInt"))+1)==as.numeric(row.names(y)), ]
## variogram calculation for south direction
v <- variogram(zinc ~ 1, data = south, cutoff = 4789.86784786386, width =
90,
cloud = T, alpha = c(0), tol.hor = 90/3)
## extract point pairs, which include center point of polygon
vS <- v[((v$np%/% attr(v, ".BigInt"))+1)==as.numeric(row.names(y)) |
((v$np%% attr(v, ".BigInt"))+1)==as.numeric(row.names(y)), ]
## convert variogramCloud to gstatVariogram and change np-entry to plot the
## variograms
class(vwest) = c("gstatVariogram", "data.frame")
class(veast) = c("gstatVariogram", "data.frame")
class(vN) = c("gstatVariogram", "data.frame")
class(vS) = c("gstatVariogram", "data.frame")
vwest$np = rep(1, nrow(vwest))
veast$np = rep(1, nrow(veast))
vN$np = rep(1, nrow(vN))
vS$np = rep(1, nrow(vS))
## split anisotropic variograms by direction
pwest = split(vwest, vwest$dir.hor)
peast = split(veast, veast$dir.hor)
## create list of directions. Order: N, NE, E, SE, S, SW, W, NW
v = list(vN, peast[[1]], peast[[2]], peast[[3]], vS, pwest[[1]],
pwest[[2]],
pwest[[3]])
## fit variograms
auto = lapply(seq(v), function(i) {
autofit(experimental_variogram = v[[i]], diagonal =
4789.86784786386,
boundaries = seq(1, 4789.86784786386, 90), tol.hor = 90/3)
})
## extract row, which includes the sill and range values
val <-lapply(seq(auto), function(i) {
auto[[i]]$var_model[2,]
})
## assign directions to rows
dir <- c(0, 45, 90, 135, 180, 225, 270, 315)
srml = do.call("rbind", val)
srml = cbind(dir, srml)
srmr = as_radians(srml$dir)
## plot range values in radial plot, see the error in NE-direction
radial.plot(as.vector(srml$range), radial.pos = srmr, rp.type = "p",
start = pi/2, radial.lim = c(0, max(srml$range)),
labels = c("N", "NE", "E", "SE", "S", "SW", "W", "NW"),
clockwise = T, line.col="red",
grid.unit = "m", show.grid.labels = 3, cex.lab=0.5)
Don't mind the too less point pairs in some directions: in my original code
with my original data basis there are at least 600 point pairs in each
variogram, so I doubt that the fault is caused by data basis.
Does anyone knows how to fix my code, so that the results are reliable?
If you find any mistakes in the code which cause this strange calculation
of autofitVariogram, please tell me. Maybe it is not allowed to split
variograms by direction and then fit separately?
Many thanks in advance for any answer!
Laura
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
