automap::autoKrige fit.method 의존성 문제
R2020. 7. 18. 23:30
공간 보간 기법 중 하나인 Kriging 을 R에서는 사용하기 편리하기 위해서 automap패키지에서 autoKrige로 지원하고 있다. 하지만 최근 fitting 과정에서 fit.method가 gstats 패키지에 있어 옵션이 수정되지 않는 현상을 겪었다.
따라서 아래와같이 function을 올려보겠다.
autoKrige2=function (formula, input_data, new_data, data_variogram = input_data,
block = 0, model = c("Sph", "Exp", "Gau", "Ste"), kappa = c(0.05,
seq(0.2, 2, 0.1), 5, 10), fix.values = c(NA, NA, NA),
remove_duplicates = TRUE, verbose = FALSE, GLS.model = NA, fit.method=fit.method,
start_vals = c(NA, NA, NA), miscFitOptions = list(), ...)
{
if (inherits(formula, "SpatialPointsDataFrame")) {
input_data = formula
formula = as.formula(paste(names(input_data)[1], "~ 1"))
}
if (!inherits(input_data, "SpatialPointsDataFrame") | !inherits(data_variogram,
"SpatialPointsDataFrame")) {
stop(paste("\nInvalid input objects: input_data or data_variogram not of class 'SpatialPointsDataFrame'.\n\tClass input_data: '",
class(input_data), "'", "\n\tClass data_variogram: '",
class(data_variogram), "'", sep = ""))
}
if (as.character(formula)[3] != 1 & missing(new_data))
stop("If you want to use Universal Kriging, new_data needs to be specified \n because the predictors are also required on the prediction locations.")
if ("newdata" %in% names(list(...)))
stop("The argument name for the prediction object is not 'newdata', but 'new_data'.")
if (remove_duplicates) {
zd = zerodist(input_data)
if (length(zd) != 0) {
warning("Removed ", length(zd)/2, " duplicate observation(s) in input_data:",
immediate. = TRUE)
print(input_data[c(zd), ])
input_data = input_data[-zd[, 2], ]
}
}
col_name = as.character(formula)[2]
if (length(unique(input_data[[col_name]])) == 1)
stop(sprintf("All data in attribute '%s' is identical and equal to %s\n Can not interpolate this data",
col_name, unique(input_data[[col_name]])[1]))
if (missing(new_data))
new_data = create_new_data(input_data)
p4s_obj1 = proj4string(input_data)
p4s_obj2 = proj4string(new_data)
if (!all(is.na(c(p4s_obj1, p4s_obj2)))) {
if (is.na(p4s_obj1) & !is.na(p4s_obj2))
proj4string(input_data) = proj4string(new_data)
if (!is.na(p4s_obj1) & is.na(p4s_obj2))
proj4string(new_data) = proj4string(input_data)
if (any(!c(is.projected(input_data), is.projected(new_data))))
stop(paste("Either input_data or new_data is in LongLat, please reproject.\n",
" input_data: ", p4s_obj1, "\n", " new_data: ",
p4s_obj2, "\n"))
if (proj4string(input_data) != proj4string(new_data))
stop(paste("Projections of input_data and new_data do not match:\n",
" input_data: ", p4s_obj1, "\n", " new_data: ",
p4s_obj2, "\n"))
}
variogram_object = autofitVariogram2(formula, data_variogram,
model = model, kappa = kappa, fix.values = fix.values,
verbose = verbose, GLS.model = GLS.model, start_vals = start_vals,
miscFitOptions = miscFitOptions, fit.method =fit.method)
krige_result = krige(formula, input_data, new_data, variogram_object$var_model,
block = block, ...)
krige_result$var1.stdev = sqrt(krige_result$var1.var)
result = list(krige_output = krige_result, exp_var = variogram_object$exp_var,
var_model = variogram_object$var_model, sserr = variogram_object$sserr)
class(result) = c("autoKrige", "list")
return(result)
}
autofitVariogram2=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(),fit.method = 1, ...)
{
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]
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(gstat::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,fit.method = fit.method), 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)
}
'R' 카테고리의 다른 글
딥로또 R 버전 (0) | 2020.02.10 |
---|---|
R 스케줄링 (0) | 2020.01.19 |
[ggplot2] 시각화 정리 (0) | 2019.12.08 |
[R][SQL] RMariaDB 외부 접속 설정하기 (0) | 2019.11.29 |
[공간정보오픈플렛폼]위경도 변환, 주소 변환 (0) | 2019.11.27 |