gwerrors <- function (x, vars, fp, adapt = NULL, bw, kernel, longlat = NULL, distMatrix = NULL, random=FALSE)
{
`%+%` <- function(x,y){ paste0(x,y)}
p4s <- as.character(NA)
Polys <- NULL
fp.missing <- missing(fp)
if (is(x, "SpatialPolygonsDataFrame")) {
Polys <- as(x, "SpatialPolygons")
gridded <- gridded(x)
dp <- coordinates(x)
p4s <- proj4string(x)
data <- as(x, "data.frame")
if (is.null(longlat) || !is.logical(longlat)) {
if (!is.na(is.projected(x)) && !is.projected(x)) {
longlat <- TRUE
}
else {
longlat <- FALSE
}
}
}
else if (is(x, "SpatialPointsDataFrame")) {
gridded <- gridded(x)
dp <- coordinates(x)
p4s <- proj4string(x)
data <- as(x, "data.frame")
if (is.null(longlat) || !is.logical(longlat)) {
if (!is.na(is.projected(data)) && !is.projected(data)) {
longlat <- TRUE
}
else {
longlat <- FALSE
}
}
}
else stop("x must be a Spatial Polygons or Points DataFrame")
if (is.null(longlat) || !is.logical(longlat))
longlat <- FALSE
if (is.integer(vars))
vars <- names(x)[vars]
stopifnot(is.character(vars))
x <- as.matrix(cbind(data[, vars], as.vector(data[, vars[1]] - data[, vars[2]])))
if (any(is.na(x)))
stop("x contains NAs")
nc <- ncol(x)
if (is.null(colnames(x)))
colnames(x) <- "V" %+% 1:nc
cn <- colnames(x)
n1 <- nrow(dp)
if (n1 != nrow(x))
stop("Differing row numbers between x and dp")
if (missing(fp)) {
fp <- dp
colnames(fp) <- colnames(dp)
}
else {
stopifnot(is(fp, "Spatial"))
gridded <- gridded(fp)
}
##NA of pixels in background grid are excluded
fp.mask <- !is.na(fp@data[,1])
fp.na <- is.na(fp@data[,1])
fp.coord <- coordinates(fp)
fp.mask.coord <- fp.coord[fp.mask, ]
fp.na.coord <- fp.coord[fp.na, ]
#coordiniations in all pixels (not masked yet)
fp <- coordinates(fp)
n2 <- nrow(fp.mask.coord)
if (is.null(adapt)) {
if (!missing(bw))
bw <- rep(bw, nrow(fp.mask.coord))
else stop("Bandwidth must be given for non-adaptive weights")
}
else {bw <- gw.adapt(dp = dp, fp = fp, quant = adapt, longlat = longlat)
###mapping bw###
###grid <- SpatialPointsDataFrame(coords=fp, data=data.frame(bw))
###spplot(grid)
bw <- bw[fp.mask]
}
###change here from for loop to apply func.
if(random==FALSE){
if (!missing(distMatrix)) {
dxs.all <- distMatrix
} else {
dxs.all <- spDists(dp, fp.mask.coord, longlat = longlat)
dxs.all[!is.finite(dxs.all)] <- 0
}
###mapping distance (from point 1)###
###grid <- SpatialPointsDataFrame(coords=fp.mask.coord, data=data.frame(dxs.all[1, ]))
###spplot(grid)
#weight
if(kernel == "gaussian")
wts.all <- apply(dxs.all, 1, function(x){ifelse(x>bw,0, exp(-0.5*(x/bw)^2))})
else if(kernel == "exponential")
wts.all <- apply(dxs.all, 1, function(x){ifelse(x>bw,0, exp(-x/bw))})
else if(kernel == "bisquare")
wts.all <- apply(dxs.all, 1, function(x){ifelse(x>bw,0, (1-(x/bw)^2)^2)})
else if(kernel == "tricube")
wts.all <- apply(dxs.all, 1, function(x){ifelse(x>bw,0, (1-(x/bw)^3)^3)})
else if(kernel == "boxcar")
wts.all <- apply(dxs.all, 1, function(x){ifelse(x>bw,0, 1)})
else stop("ivalid kernel type")
###mapping distance (from point 1)###
###grid <- SpatialPointsDataFrame(coords=fp.mask.coord, data=data.frame(wts.all[,1]))
###spplot(grid)
} else {
mcs <- sample(n1)
dxs.all <- spDists(dp[mcs, ], fp.mask.coord, longlat = longlat)
dxs.all[!is.finite(dxs.all)] <- 0
#weight
if(kernel == "gaussian")
wts.all <- apply(dxs.all[mcs, ], 1, function(x){ifelse(x>bw,0, exp(-0.5*(x/bw)^2))})
else if(kernel == "exponential")
wts.all <- apply(dxs.all[mcs, ], 1, function(x){ifelse(x>bw,0, exp(-x/bw))})
else if(kernel == "bisquare")
wts.all <- apply(dxs.all[mcs, ], 1, function(x){ifelse(x>bw,0, (1-(x/bw)^2)^2)})
else if(kernel == "tricube")
wts.all <- apply(dxs.all[mcs, ], 1, function(x){ifelse(x>bw,0, (1-(x/bw)^3)^3)})
else if(kernel == "boxcar")
wts.all <- apply(dxs.all[mcs, ], 1, function(x){ifelse(x>bw,0, 1)})
else stop("ivalid kernel type")
}
if( any(wts.all[which(wts.all < 0) | is.na(wts.all) | sum(wts.all)==0])){
stop(print("Invalid weights"))
}
cov.wt2_run <- function(y){
out <- tryCatch(error_diagnostics(as.matrix(x), y) , error =function(e){
ans <- c(rmse = NA, mae = NA, difference = NA, calc.pointn = NA, cor.test = NA, cor.pval = NA
#,mse = mse, bias2 = bias2, variance = v
)
return(ans)}
)}
##cov.wt2 includes the calc of rmse & difference
res.error <- apply(wts.all, 1, function(y)(cov.wt2_run(y)))
#colname of res.error
rnm <- c( "rmse_" %+% cn[2] %+% "_" %+% cn[1],
"mae_" %+% cn[2] %+% "_" %+% cn[1],
"difference_" %+% cn[2] %+% "_" %+% cn[1],
"calcpointn_" %+% cn[2] %+% "_" %+% cn[1],
"cor_" %+% cn[2] %+% "_" %+% cn[1],
"cor.pvalue_" %+% cn[2] %+% "_" %+% cn[1]
)
rownames(res.error) <- rnm
SDF_fill <- SpatialPointsDataFrame(coords = fp.mask.coord,
data = data.frame(t(res.error)),
proj4string = CRS(p4s))
if(dim(fp.na.coord)[1]>0){
SDF_NA <- SpatialPointsDataFrame(coords = fp.na.coord,
data = data.frame(matrix(NA,nrow=dim(fp.na.coord)[1],ncol=dim(t(res.error))[2] )),
proj4string = CRS(p4s))
names(SDF_NA) <- names(SDF_fill)
SDF <- rbind(SDF_fill, SDF_NA)
} else {
SDF <- SDF_fill
}
if (gridded)
gridded(SDF) <- TRUE
else if (!is.null(Polys) && fp.missing) {
df <- data.frame(SDF@data)
rownames(df) <- sapply(slot(Polys, "polygons"), function(i) slot(i, "ID"))
SDF <- SpatialPolygonsDataFrame(Sr = Polys, data = df)
}
res <- list(SDF = SDF, bandwidth = bw, adapt = adapt)
res
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.