2

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?

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

3
  • You can, in fact, pass a raster to the autocor function. The window (matrix) defining the size of the weights matrix (Wij) is passed to the "w" argument. Per the help, the "x" argument can either be a vector or SpatRast object. Rasters are always numeric, even if dummy coded for nominal data. Commented Jul 20, 2023 at 15:10
  • @JeffreyEvans Well, if you pass just the SpatRaster to terra::autocor() you'll get the error: "Error in match.arg(tolower(method), c("moran", "geary")) : 'arg' should be one of “moran”, “geary”".
    – Corbjn
    Commented Jul 20, 2023 at 15:15
  • Unless the developer has mad substantive changes to this function, this is a bug in the current version. It is just not accepting anything other than "moran" or "geary" (global options) to the method argument. I will report it on GitHub. Commented Jul 20, 2023 at 15:46

2 Answers 2

2

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

1
  • While this works I used spdep in the end.
    – Corbjn
    Commented Jul 31, 2023 at 9:14
0

Working solution using spdep:

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


r <- terra::rast(ncol = 500, nrow = 500, val = rnorm(500 * 500))
nb <- spdep::cell2nb(
    terra::nrow(r),
    terra::ncol(r),
    type = "queen",
    torus = TRUE
    )
r_val <- terra::values(r)
spatcor <- spdep::localmoran_perm(
    as.vector(r_val),
    spdep::nb2listw(nb, style = "B"),
    nsim = 999
    )
r_spatcor <- terra::rast(
    nrows = terra::nrow(r),
    ncols = terra::ncol(r),
    extent = terra::ext(r),
    vals = as.vector(spatcor[, "Ii"])
)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.