Nothing
interpolate = function(observations, predictionLocations,
outputWhat = list(mean = TRUE, variance = TRUE), obsChar = NA,
methodName = "automatic", maximumTime = 30, optList = list(), cv = FALSE, ...) {
startTime = Sys.time()
# save.image("debug.img")
nPred = ifelse(is.numeric(predictionLocations), predictionLocations,
if (inherits(predictionLocations, "Spatial")) nrow(coordinates(predictionLocations)) else dim(predictionLocations)[1])
msg = paste("R", startTime, "interpolating", nrow(observations), "observations,",
nPred, "prediction locations\n")
# npred, "prediction locations", ifelse(cv, ", crossvalidating", " "), "\n")
cat(msg)
obsIsSf = FALSE
if (inherits(observations, "sf")) {observations = as(observations, "Spatial"); obsIsSf = TRUE}
if (inherits(predictionLocations, "sf")) predictionLocations = as(predictionLocations, "Spatial")
if ("formulaString" %in% names(optList)) {
formulaString = as.formula(optList$formulaString)
} else {
formulaString = as.formula("value~1")
}
localParams = list(confProj = TRUE,debug.level = 0)
localParams = modifyList(localParams,optList)
debug.level = localParams$debug.level
if (debug.level >=1) {
print("sending debug info to processDescription")
debugOutput = 1
rm(debugOutput)
debugOut = textConnection("debugOutput","w")
sink(debugOut,split=TRUE)
}
if (methodName == "automatic" & !(sum(names(optList) %in%
c("methodParameters", "variogramModel", "copulaParams"))) >= 1) {
methodName = chooseMethod(observations, predictionLocations,formulaString,
obsChar, maximumTime, outputWhat = outputWhat)
} else if (is.finite(maximumTime)) {
# predTime = try(predictTime("spatialPredict",nObs= dim(observations)[1],
# nPred = dim(coordinates(predictionLocations))[1], formula = value~1, class = methodName, outputWhat = outputWhat))
predTime = predictTime(nObs= dim(observations)[1],
nPred = nPred, formulaString = formulaString,
class = methodName, outputWhat = outputWhat, FUN = "spatialPredict")
if (is.na(predTime)) {
warning(paste("was not able to estimate prediction time for methodName",methodName))
} else if (maximumTime < predTime) {
warning(paste("the choosen method",methodName,
"is not able to give predictions within the maximum time of",maximumTime,"seconds"))
}
}
# set up intamap object:
krigingObject = createIntamapObject(
observations = observations,
predictionLocations = predictionLocations,
formulaString = formulaString,
outputWhat = outputWhat,
obsChar = obsChar,
params = getIntamapParams(localParams),
class = methodName )
# check:
debug.level = krigingObject$params$debug.level
checkSetup(krigingObject)
# do interpolation steps:
krigingObject = preProcess(krigingObject)
if (is.null(krigingObject$variogramModel) && is.null(krigingObject$copulaParams)
&& is.null(krigingObject$inverseDistancePower)) {
krigingObject = estimateParameters(krigingObject, ...)
}
krigingObjectMp = try(methodParameters(krigingObject))
if (!is(krigingObjectMp,"try-error")) krigingObject = krigingObjectMp
nsim = ifelse("nsim" %in% names(outputWhat), outputWhat$nsim, 0)
# Methods able to create simulations
if (cv) {
kObj = krigingObject
predictions = krigingObject$observations
depVar = as.character(krigingObject$formulaString[[2]])
predictions@data = data.frame(var1.pred = NA, var1.var = NA,
observed = observations@data[,depVar], residual = NA, zscore = NA)
for (i in 1:dim(krigingObject$observations)[1]) {
kObj$predictionLocations = krigingObject$observations[i,]
kObj$observations = krigingObject$observations[-i,]
if (debug.level == 0) {
tmp = capture.output(kObj <- spatialPredict(kObj, ...))
} else kObj <- spatialPredict(kObj)
if ("var1.pred" %in% names(kObj$predictions) & "var1.var" %in% names(kObj$predictions)) {
predictions@data[i,1:2] = kObj$predictions@data[,c("var1.pred", "var1.var")]
} else predictions@data[i,1:2] = kObj$predictions@data[,c("mean", "variance")]
}
predictions$residual = predictions$observed - predictions$var1.pred
predictions$zscore = predictions$residual/sqrt(predictions$var1.var)
krigingObject$predictions = predictions
} else krigingObject = spatialPredict(krigingObject, nsim = nsim)
krigingObject = postProcess(krigingObject)
# Add plot if wanted
if (!is.null(krigingObject$returnPlot) && krigingObject$returnPlot)
krigingObject$processPlot = ""#createPlot(krigingObject)
else
krigingObject$processPlot = ""
# Add to process description which method we used
krigingObject$processDescription =
paste("Spatial prediction using the method ",class(krigingObject) )
# Add debug information if wanted
if (debug.level >=1) {
sink()
close(debugOut)
krigingObject$processDescription = c(krigingObject$processDescription,debugOutput)
}
# Create a table easier to handle for the WPS
krigingObject$outputTable = toJava(krigingObject$outputTable)
attr(krigingObject$outputTable, "transposed") = TRUE
if (obsIsSf) {
krigingObject$observations = as(krigingObject$observations, "sf")
krigingObject$predictionLocations = as(krigingObject$predictionLocations, "sf")
}
return(krigingObject)
}
interpolateBlock = function(observations, predictionLocations, outputWhat, blockWhat = "none",
obsChar = NA, methodName = "automatic", maximumTime = 30, optList = list(), ...) {
obsIsSf = FALSE
if (inherits(observations, "sf")) {
obsIsSf = TRUE
observations = as(observations, "Sptial")
predictionLocations = as(predictionLocations, "Sptial")
}
startTime = Sys.time()
msg = paste("R", startTime, "interpolating", nrow(coordinates(observations)), "observations,",
nrow(coordinates(predictionLocations)), "prediction locations\n")
cat(msg)
if ("formulaString" %in% names(optList)) {
formulaString = as.formula(optList$formulaString)
} else {
formulaString = as.formula("value~1")
}
localParams = list(confProj = TRUE,debug.level = 0)
localParams = modifyList(localParams,optList)
debug.level = localParams$debug.level
if (debug.level >=1) {
debugOutput = ""
rm(debugOutput)
debugOut = textConnection("debugOutput","w")
sink(debugOut, split=TRUE)
}
if (inherits(observations, "Spatial")) {
nObs = dim(coordinates(observations))[1]
nPred = dim(coordinates(predictionLocations))[1]
} else if (inherits(observations, "sf")) {
nObs = dim(observations)[1]
nPred = dim(predictionLocations)[1]
}
if (methodName == "automatic") {
methodName = chooseMethod(observations,predictionLocations,formulaString,obsChar,maximumTime, outputWhat)
} else {
predTime = predictTime(nObs = nObs,
nPred = nPred, formulaString = formulaString,
class = methodName, outputWhat = outputWhat, FUN = "spatialPredict")
if (is.na(predTime)) {
warning(paste("was not able to estimate prediction time for methodName",methodName))
} else if (maximumTime < predTime) {
warning(paste("the choosen method",methodName,
"is not able to give predictions within the maximum time of",maximumTime,"seconds"))
}
}
if (inherits(predictionLocations, "Spatial")) {
if (!is(predictionLocations,"SpatialPolygons")) {
if (!is(predictionLocations,"SpatialGrid")) gridded(predictionLocations) = TRUE
cellsize = predictionLocations@grid@cellsize
block = cellsize
ptext = paste("blocks of size = ",block)
} else {
block = numeric(0)
ptext = "SpatialPolygons"
}
} else {
}
localParams$block = block
krigingObject = createIntamapObject(
observations = observations,
predictionLocations = predictionLocations,
formulaString = formulaString,
outputWhat = outputWhat,
blockWhat = blockWhat,
obsChar = obsChar,
params = getIntamapParams(localParams),
class = methodName)
# blockWhat is only necessary for prediction types that only exist for blocks
# e.g. fraction above threshold or max within block
# check:
debug.level = krigingObject$params$debug.level
checkSetup(krigingObject)
# do interpolation steps:
krigingObject = preProcess(krigingObject)
if (is.null(krigingObject$variogramModel) && is.null(krigingObject$copulaParams)
&& is.null(krigingObject$inverseDistancePower)) {
krigingObject = estimateParameters(krigingObject, ...)
}
krigingObjectMp = try(methodParameters(krigingObject))
if (!is(krigingObjectMp,"try-error")) krigingObject = krigingObjectMp
krigingObject = spatialPredict.block(krigingObject, ...)
krigingObject = postProcess(krigingObject)
# Add plot if wanted
if (!is.null(krigingObject$returnPlot) && krigingObject$returnPlot)
krigingObject$processPlot = createPlot(krigingObject)
else
krigingObject$processPlot = ""
# Add to process description which method we used
krigingObject$processDescription =
c(paste("Spatial block prediction using the method ",class(krigingObject)),
paste("Method applied on",ptext))
# Add debug information if wanted
if (debug.level >=1) {
sink()
close(debugOut)
krigingObject$processDescription = c(krigingObject$processDescription,debugOutput)
}
# Create a table easier to handle for the WPS
krigingObject$outputTable = toJava(krigingObject$outputTable)
if (obsIsSf) {
krigingObject$observations = as(krigingObject$observations, "sf")
krigingObject$predictionLocations = as(krigingObject$predictionLocations, "sf")
}
attr(krigingObject$outputTable, "transposed") = TRUE
return(krigingObject)
}
toJava = function(obj) {
obj = data.frame(obj)
obj = as.matrix(obj)
obj = t(obj)
return(obj)
}
createPlot = function(krigingObject) {
# if (sessionInfo()$R.version$os == "mingw32")
# if (sessionInfo()$R.version$os != "linux-gnu")
# return("empty")
if (!is.null(krigingObject$variogramModel)
&& !is.null(krigingObject$sampleVariogram)) {
fn = tempfile()
# Using pdf for the moment, more platform independent
pdf(fn)
# svg(fn)
print(plot(krigingObject$sampleVariogram, krigingObject$variogramModel,
main = "Sample variogram and fitted model"))
dev.off()
ret = readLines(file(fn, "r"))
unlink(fn)
return(ret)
} else
return("")
}
chooseMethod = function(observations, predictionLocations, formulaString,
obsChar, maximumTime,outputWhat) {
methodNames = c("copula","automap")
nPred = ifelse(is.numeric(predictionLocations), predictionLocations,
nrow(coordinates(predictionLocations)))
if (length(obsChar) > 0 && !is.na(obsChar) ) { #} && requireNamespace("psgp")) {
pTime = predictTime(nObs= dim(observations)[1],
nPred = nPred, formulaString = formulaString,
class="psgp", outputWhat = outputWhat)
if (pTime < maximumTime) return("psgp") else stop("no method is able to give predictions within the maximum time given")
} else {
for (i in 1:length(methodNames)) {
pTime = predictTime(nObs= dim(observations)[1],
nPred = nPred, formulaString = formulaString,
class = methodNames[i], outputWhat = outputWhat)
print(paste("estimated time for ",methodNames[i],pTime))
if (methodNames[i] == "copula") {
dataObs = observations[[as.character(formulaString[[2]]) ]]
test = isNonGauss(dataObs)
if (test) {
if (!is.na(pTime) && pTime < maximumTime && !("nsim" %in% names(outputWhat))) {
return("copula")
} else if (min(dataObs) <=0) {
t1 = "Data are non-gaussian with zero-values and/or negative values"
t2 = "time not sufficient for copulas (predicted time:"
t3 = paste(round(pTime,2), "maximumTime:",maximumTime,").")
t4 = "Will use automap."
warning(paste(t1,t2,t3,t4))
return("automap")
} else return("transGaussian")
}
}
if (pTime < maximumTime) return(methodNames[i])
}
stop("no method is able to give predictions within the maximum time given")
}
}
isNonGauss = function(dataObs, which = FALSE) {
test = logical(4)
md = min(dataObs)
if (md <= 0)
dataObs = dataObs + abs(md) + sd(dataObs)
fn = fivenum(dataObs)
iqr = IQR(dataObs)
test[1] = length(boxplot.stats(dataObs)$out)/length(dataObs) > 0.1
test[2] = fn[3] - fn[2] < iqr/3
test[3] = fn[4] - fn[3] < iqr/3
g=boxcox(dataObs ~ 1,lambda=seq(-2.5,2.5,len=101),plotit=FALSE)$y
test[4] = g[71] < sort(g)[91]
if (which)
return(list(newObs = dataObs, test=test))
else
return(any(test))
}
#chooseMethod(meuse,predictionLocations,as.formula(value~1),obsChar,Inf)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.