Statistically comparing two spatially congruent continuous rasters?

Question

First: Rasters are the same spatial resolution & extent and contain the same kind of information. Each pixel contains a probability value from 0-1.

Sounds simple to compare. Conventional wisdom might point to a neighborhood analysis of some size (along with a sensitivity analysis of Nhood size). I could always use a pearson linear correlation, a reflective correlation, a spearman rank correlation, etc. However, none of these analyses will quite give me what I want. Since these two rasters contain extremely similar values for most locations (pixels) what I'm really interested in knowing:

How do they differ in relation to isotropy?

I suspect there are small directional differences in the values of these two rasters. To be more direct, there are two possible comparisons to consider; (1) values, and (2) position. Most of the analyses listed above provide information about values and not position.

Is there a way to compare 2 semivariograms?

Answer 1

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.