rgstatautomapgeostatistics
autofitVariogram in automap: cutoff option doesn't work
I'd like to create a spatial separation distance up to which point pairs are included in semivariance estimates (cutoff function in variogram {gstat}), but using autofitVariogram in automap package. Despite the use of miscFitOptions function nothing happened (error or expected output). In my example below I'd like to cutoff at 1000m the meuse data set:
# Packages
library(automap)
library(gstat)
# Classical meuse dataset example
data(meuse)
coordinates(meuse) = ~x+y
# Funcion autofitVariogram
autoZn=autofitVariogram(log(zinc)~1, meuse)
summary(autoZn)
# Plot variogram
plot(autoZn, pch=19, col="black")
# Now with 1000 meters cutoff
autoZn_cut=autofitVariogram(log(zinc)~1, meuse, cutoff=1000)
summary(autoZn_cut)
plot(autoZn_cut, pch=19, col="black")
# or
autoZn_cut=autofitVariogram(log(zinc)~1, meuse, miscFitOptions = list(cutoff=1000))
summary(autoZn_cut)
plot(autoZn_cut, pch=19, col="black")
But in the three plots do not change anything and I don't have any error?
Please, any help with it?
Solution
I modified the autofitVariogram function (I call it my_autofitVariogram) adding a boundaries option.
my_autofitVariogram <- function (formula, input_data, 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(), boundaries=NULL, ...) {
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)
longlat = !is.projected(input_data)
if (is.na(longlat))
longlat = FALSE
diagonal = spDists(t(bbox(input_data)), longlat = longlat)[1, 2]
if (is.null(boundaries)) {
boundaries = c(2, 4, 6, 9, 12, 15, 25, 35, 50, 65, 80, 100) * diagonal * 0.35/100
}
if (!is(GLS.model, "variogramModel")) {
experimental_variogram = variogram(formula, input_data, boundaries = boundaries, ...)
}
else {
if (verbose)
cat("Calculating GLS sample variogram\n")
g = gstat(NULL, "bla", formula, input_data, model = GLS.model,
set = list(gls = 1))
experimental_variogram = variogram(g, boundaries = boundaries,
...)
}
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(GLS.model, "variogramModel")) {
experimental_variogram = variogram(formula, input_data,
boundaries = boundaries, ...)
}
else {
experimental_variogram = variogram(g, boundaries = 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)
}
You can copy my_autofitVariogram in a file (named my_autofitVariogram.r), and put it in your working directory. Then, run this sample code:
library(automap)
library(gstat)
source("my_autofitVariogram.r")
data(meuse)
coordinates(meuse) = ~x+y
# Function my_autofitVariogram
autoZn <- my_autofitVariogram(log(zinc)~1, meuse,
boundaries=c(2, 4, 6, 9, 12, 15, 25, 35, 50, 65, 80, 100)*10
)
summary(autoZn)
plot(autoZn, pch=19, col="black")
- Is there a way to pass a string as a variable/column name to my function and use in a call to mutate?
- Connecting R Leaflet to esri vector tile service
- Add a secundary reverse x-axis ggplot2
- how sort in reverse order xts object
- ifelse() stripping POSIXct attribute from vector of timestamps?
- How to mutate a year quarter str into date?
- Change text of full_screen button ("Expand") in bslib card
- Returning a matrix after removing duplicate combinations of elements
- How can I get fine control over the position of a ggplot legend, now that a numeric 'legend.position' argument is deprecated?
- Adding percentage instead of observed counts on y axis on R figure
- patchwork - inset_element(): size and position
- R optimize function not reaching optimal value
- Why is E(g) not identical to itself
- control x position using ggrepel
- YAML parser error when rendering Quarto document
- Filter list of lists with purrr::keep
- How can I deal with the Lavaan error: "syntax error in lavaan model syntax"?
- Strange behavior of pickerInput, how to resolve it?
- R diagrammeR using html for formatting while reading text from r variable
- How to get a warning on "shiny app will not work if the same output is used twice"
- Condition facetting of ggplot2 plot
- Can we retrieve entries from a tibble using index matrix?
- Adding conditions to loop to generate class intervals using ClassInt in R for areas with fewer n classes
- ggsurvplot - reposition strata color indicators of risk table
- How to rearrange value boxes at a certain screen size?
- Remove last word from string
- Are doubly-linked structures allowed?
- Mermaid diagram line break
- Using gsub to extract only capital letters of a certain length
- How do I style a shinyWidgets::multiInput()