Skip to content

Commit

Permalink
Merge pull request #3 from sigmafelix/action-test
Browse files Browse the repository at this point in the history
styling
  • Loading branch information
Insang Song authored May 1, 2024
2 parents cb491a3 + df98b69 commit 08dc895
Show file tree
Hide file tree
Showing 11 changed files with 1,001 additions and 939 deletions.
File renamed without changes.
157 changes: 82 additions & 75 deletions R/autoKrige.r
Original file line number Diff line number Diff line change
@@ -1,88 +1,95 @@
autoKrige = 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,
start_vals = c(NA,NA,NA), miscFitOptions = list(), ...)
autoKrige <- 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,
start_vals = c(NA, NA, NA), miscFitOptions = list(), ...)
# This function performs an automatic Kriging on the data in input_data
{
if(inherits(formula, "SpatialPointsDataFrame"))
# Is someone just passes a spatialpointsdataframe, assume he/she wants to interpolate the first column with Ordinary Kriging
{
input_data = formula
formula = as.formula(paste(names(input_data)[1], "~ 1"))
}

# Check if inpu_data and data_variogram are SpatialPointsDataFrame
if(!inherits(input_data,"SpatialPointsDataFrame") | !inherits(data_variogram,"SpatialPointsDataFrame"))
if (inherits(formula, "SpatialPointsDataFrame"))
# Is someone just passes a spatialpointsdataframe, assume he/she wants to interpolate the first column with Ordinary Kriging
{
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=''))
input_data <- formula
formula <- as.formula(paste(names(input_data)[1], "~ 1"))
}

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'.")

# Check if there are points or gridcells on the exact same coordinate and provide a more informative error message.
# Points on the same spot causes the interpolation to crash.
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], ]

}

# Check if inpu_data and data_variogram are SpatialPointsDataFrame
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'.")

# Check if there are points or gridcells on the exact same coordinate and provide a more informative error message.
# Points on the same spot causes the interpolation to crash.
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], ]
}
}

# If all the values return an informative error
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 all the values return an informative error
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)

## Perform some checks on the projection systems of input_data and new_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"))
}
if (missing(new_data)) new_data <- create_new_data(input_data)

# Fit the variogram model, first check which model is used
variogram_object = autofitVariogram(formula,
data_variogram,
#autoselect.model = autoselect.model,
model = model,
kappa = kappa,
fix.values = fix.values,
verbose = verbose,
GLS.model = GLS.model,
start_vals = start_vals,
miscFitOptions = miscFitOptions)
## Perform some checks on the projection systems of input_data and new_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"
))
}
}

## Perform the interpolation
krige_result = krige(formula,
input_data,
new_data,
variogram_object$var_model,
block = block,
...)
# Fit the variogram model, first check which model is used
variogram_object <- autofitVariogram(formula,
data_variogram,
# autoselect.model = autoselect.model,
model = model,
kappa = kappa,
fix.values = fix.values,
verbose = verbose,
GLS.model = GLS.model,
start_vals = start_vals,
miscFitOptions = miscFitOptions
)

krige_result$var1.stdev = sqrt(krige_result$var1.var)
## Perform the interpolation
krige_result <- krige(formula,
input_data,
new_data,
variogram_object$var_model,
block = block,
...
)

# Aggregate the results into an autoKrige object
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")
krige_result$var1.stdev <- sqrt(krige_result$var1.var)

return(result)
# Aggregate the results into an autoKrige object
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)
}
Loading

0 comments on commit 08dc895

Please sign in to comment.