Local Moran’s I using terra::autocor()

rterra

I have a SpatRaster with multiple layers and I want to calculate the Local Moran's I using R's terra package for each layer. Following terra's documentation I found the function terra::autocor(). It says:

If x is numeric or SpatRaster: "moran" for Moran’s I and "geary" for
Geary’s C. If x is numeric also: "Gi", "Gi*" (the Getis-Ord statistics), locmor
(local Moran’s I) and "mean" (local mean)

What I understand from that is that I can't just apply "locmor" to the SpatRaster but instead I have to convert it to numeric. I thought terra::values() would do it but I get the following error message:

Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘autocor’ for signature ‘"matrix"’

What am I doing wrong here? Thanks in advance!
Reproducible example:

library("terra")
library("spdep")

w <- matrix(c(0, 1, 0, 1, 0, 1, 0, 1, 0), nrow = 3)
r <- terra::rast(ncol = 500, nrow = 500, val = rnorm(500 * 500))
r_auto <- terra::autocor(terra::values(r), w, method = "locmor")

System Information:
Ubuntu 22.04.2 (64 bit), terra 1.7.39, spdep 1.2-8

Best Answer

Looked at the source on GitHub and not an error. Looks like the current implementation only takes numeric with an associated binary weights matrix indicating contingency (based on vector data). For local statistics, there is no method for SpatRaster and it looks like the function that was in the raster package is no longer exported. This may be by design because a moving window LISA has always been suspect as it does not integrate to the global I statistic. However, for exploratory analysis, it has been useful. Here is the original source from raster.

MoranLocal <- function(x, w=matrix(c(1,1,1,1,0,1,1,1,1),3,3)) { 
  z  <- x - cellStats(x, mean) 
    if (sum(! unique(w) %in% 0:1) > 0) {
      zz <- calc(z, fun=function(x) ifelse(is.na(x), NA ,1))
      W  <- focal( zz, w=w, na.rm=TRUE, pad=TRUE)       
    } else {
      w2 <- w
      w2[w2==0] <- NA
      W  <- focal( z, w=w2, fun=function(x, ...){ sum(!is.na(x)) }, na.rm=TRUE, pad=TRUE)
    }
    lz <- focal(z, w=w, na.rm=TRUE, pad=TRUE) / W   
    n <- ncell(x) - cellStats(x, 'countNA')
    s2 <- cellStats(x, 'sd')^2 
    (z / s2) * lz
}

You could coerce to a raster object and pass to this function but, this is a short term solution.

library(raster)
library(terra)

r <- rast(nrows=100, ncols=100, xmin=0)
  r[] <- runif(ncell(r))
    r <- focal(r, 9, mean)

lI <- MoranLocal(raster(r))
  r <-c(r,rast(lI))
    names(r) <- c("random Gaussian","Local Autocorrelation")
plot(r)

local autocorr

Related Question