# Faster way to compute Moran’s Index from a very large distance matrix between geographic points

I am computing Moran’s Index from a very large distance matrix between geographic points. The time to execute the code is too high. I used the `dism` function in the `geosphere` package.

Here a reproducible example in which the Moran's index is calculated for three covariates:

``````## Generate a dataset
N = 100000
dat <- data.frame(long = runif(N, min = -130, max = -60), lat = runif(N, min = 16, max = 60),
var1 = runif(N, min = 0, max = 40),
var2 = runif(N, min = 0, max = 40),
var3 = runif(N, min = 0, max = 40))

## Generate a distance matrix
mat_dist <- as.matrix(geosphere::distm(cbind(dat[,c("long")], dat[,c("lat")])))/1000 ## distance in km
## summary(mat_dist)

## Define the data frame containing the Moran's index for each covariate
moranI <- data.frame(matrix(NA, nrow = 1, ncol = 3))
names(moranI) <- paste0("var", 1:3)
## print(moranI)

## Compute the Moran's index for each covariate
ncol_moranI <- 1
while(ncol_moranI <= dim(moranI)[2]){

## Modify the distance matrix
mat_dist_moranI <- mat_dist

## Generate a matrix of inverse distance weights
mat_dist_moranI <- 1/mat_dist_moranI

## Replace -Inf with 0
diag(mat_dist_moranI) <- 0

## Calculate the Moran's index
moranI[ncol_moranI] <- abs(ape::Moran.I(dat[, c(names(moranI)[ncol_moranI])], mat_dist_moranI, scaled = TRUE)\$observed)

## Update the incremental variable
ncol_moranI <- ncol_moranI + 1
}
## print(moranI)
``````

I tested the `distance` function in the `terra` package, but the execution time is also too high.

Is there any way to speed up the code ? Any help would be greatly appreciated. Thanks so much.

• "The time to execute the code is too high." is a bit vague - how long does it currently take for a given N and what would be acceptable? Is it purely the distance calculation that takes all the time? Or is it the Moran computation? Apr 26, 2023 at 7:44
• Thanks Spacedman. It's the distance calculation that takes all the time. For N = 1000, the time execution is 1.36s and N = 10000, it's 137.63 s (a bit too long) Apr 26, 2023 at 15:46
• Is 10,000 the most you need? Because your example has `N = 100000` - ten times that... Apr 26, 2023 at 15:59
• I need to run the code using a dataframe with 100,000 rows. Apr 26, 2023 at 16:07
• If this works for you I think we can close this Q as a duplicate! gis.stackexchange.com/questions/291304/… from github.com/mcooper/moranfast although it doesn't do geospatial distance... Apr 26, 2023 at 20:30