Spatial process can be introduced by the underlying data that went into a probabilistic estimate, the model itself and the grain of the raster. You also have the addition of potentially of introduced nonstationarity in the random field. How do you plan on teasing out these confounding factors in understanding the anisotropy of the data? What exactly is the hypothesis that is being tested here?
There has been some interesting work done on using ordination techniques to evaluate the uncertainty in spatial estimates. Mostly it has been applied to univariate estimates but, one could expand it to include a bivariate (or even multivariate) assessment. Here is a worked example using some correlated dummy data.
References
Correia, P., L. Azevedo, R.F.M. Nunes, G. Schwedersky neto (2014) Multidimensional scaling for the evaluation of a geostatistical seismic elastic inversion methodology. Geophysics 79(1):M1-M10 DOI10.1190/geo2013-0037.1
Scheidt, C., and J. Caers (2009) Representing Spatial Uncertainty Using Distances and Kernels. Mathematical Geosciences 41:397.
Add required libraries and generate some data
library(sp)
library(spdep)
library(raster)
library(gstat)
library(sp)
data(meuse)
data(meuse.grid)
coordinates(meuse) <- ~x + y
coordinates(meuse.grid) <- ~x + y
gridded(meuse.grid) = TRUE
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
proj4string(meuse) <- CRS("+init=epsg:28992")
v1 <- variogram(log(copper) ~ 1, meuse)
x1 <- fit.variogram(v1, vgm(1, "Sph", 800, 1))
x <- krige(zinc ~ 1, meuse, meuse.grid, x1, nmax = 30)
gridded(x) <- TRUE
x@data = as.data.frame(x@data[,-2])
x@data[,1] <- x@data[,1] / max(x@data[,1])
v2 <- variogram(log(elev) ~ 1, meuse)
x2 <- fit.variogram(v2, vgm(.1, "Sph", 1000, .6))
y <- krige(elev ~ 1, meuse, meuse.grid, x2, nmax = 30)
gridded(y) <- TRUE
y@data <- as.data.frame(y@data[,-2])
y@data[,1] <- y@data[,1]
y@data[,1] <- y@data[,1] / max(y@data[,1])
x@data$y <- y@data[,1]
names(x)[1] <- "x"
Now, we can specify some needed parameters
cs=40 # cell size of raster array
d = sqrt(2*((cs*3)^2)) # neighborhood distance
x.idx = 1 # index of value x column in raster
y.idx = 2 # index of value y column in raster
Here we create a neighbors list using the d parameter
nb <- spdep::dnearneigh(sp::coordinates(x),0, d)
Now, we can loop through the neighbors list to subset local neighborhood values and calculate the local Multidimensional Scaling (MDS) for the center cell of the matrix.
xy.mds <- vector()
for(i in 1:length(nb)) {
if( ncol(x) > 1) {
x.var <- x@data[nb[[i]],][x.idx][,1]
y.var <- x@data[nb[[i]],][y.idx][,1]
} else {
x.var <- x@data[nb[[i]],]
y.var <- x@data[nb[[i]],]
}
cmd <- stats::cmdscale(stats::dist(cbind(x.var,y.var)), k=1)
xy.mds <- append(xy.mds, cmd[round((length(cmd)/2)+1,0)])
}
We now assign the MDS values and calculate a cell-wise absolute difference between x and y and display the results.
x@data$diff <- abs(x@data[,1] - x@data[,2])
x@data$mds <- xy.mds
x <- stack(x)
plot(x)
And, here are the results. Keep in mind that the MDS values are two-tailed so, this is not a cold (low values) to hot (high values) plot like the other three.
