Spatial correlogram using a raster package

Question

Dear Crowd

Problem

I tried to calculate the spatial correlogram with the packages nfc, pgirmess, SpatialPack and spdep. However, it was difficult for me to determine the beginning and end of the distance. I'm only interested in spatial autocorrelation at smaller distances, but at smaller cells. In addition, since the raster is quite large (1.8 megapixels), I run into memory problems with these packages, but with SpatialPack.

So I tried to create my own code using the Moran function from the raster package. But I should have some error, since the result for the full data set is slightly different from the result from other packages. If there is no error in my code, this may at least help others with similar problems.

Question

I'm not sure if my focal matrix is ​​wrong. Could you tell me if a central pixel is needed? Using testdata, I cannot show the differences between the methods, but there are differences on my complete dataset, as shown in the figure below. However, the bunkers are not exactly the same (50 m versus 69 m), so this may explain some of the differences. However, in the first box, this explanation seems implausible to me. Or maybe the irregular shape of my raster and different ways of processing NA cause a difference?

Comparing your own method with one of SpatialPack

Runable example

Testdata STRONG>

The code for calculating test data is taken from http://www.petrkeil.com/?p=1050#comment-416317

# packages used for the data generation library(raster) library(vegan) # will be used for PCNM # empty matrix and spatial coordinates of its cells side=30 my.mat <- matrix(NA, nrow=side, ncol=side) x.coord <- rep(1:side, each=side)*5 y.coord <- rep(1:side, times=side)*5 xy <- data.frame(x.coord, y.coord) # all paiwise euclidean distances between the cells xy.dist <- dist(xy) # PCNM axes of the dist. matrix (from 'vegan' package) pcnm.axes <- pcnm(xy.dist)$vectors # using 8th PCNM axis as my atificial z variable z.value <- pcnm.axes[,8]*200 + rnorm(side*side, 0, 1) # plotting the artificial spatial data r <- rasterFromXYZ(xyz = cbind(xy,z.value)) plot(r, axes=F) 

Native code

 library(raster) sp.Corr <- matrix(nrow = 0,ncol = 2) formerBreak <- 0 #for the first run important for (i in c(seq(10,200,10))) #Calculate the Morans I for these bins { cat(paste0("..",i)) #print the bin, which is currently calculated w = focalWeight(r,d = i,type = 'circle') wTemp <- w #temporarily saves the weigtht matrix if (formerBreak>0) #if it is the second run { midpoint <- ceiling(ncol(w)/2) # get the midpoint w[(midpoint-formerBreak):(midpoint+formerBreak),(midpoint-formerBreak):(midpoint+formerBreak)] <- w[(midpoint-formerBreak):(midpoint+formerBreak),(midpoint-formerBreak):(midpoint+formerBreak)]*(wOld==0)#set the previous focal weights to 0 w <- w*(1/sum(w)) #normalizes the vector to sum the weights to 1 } wOld <- wTemp #save this weight matrix for the next run mor <- Moran(r,w = w) sp.Corr <- rbind(sp.Corr,c(Moran =mor,Distance = i)) formerBreak <- i/res(r)[1]#divides the breaks by the resolution of the raster to be able to translate them to the focal window } plot(x=sp.Corr[,2],y = sp.Corr[,1],type = "l",ylab = "Moran I",xlab="Upper bound of distance") 

Other methods for calculating spatial correlogram

 library(SpatialPack) sp.Corr <- summary(modified.ttest(z.value,z.value,coords = xy,nclass = 21)) plot(x=sp.Corr$coef[,1],y = data$coef[,4],type = "l",ylab = "Moran I",xlab="Upper bound of distance") library(ncf) ncf.cor <- correlog(x.coord, y.coord, z.value,increment=10, resamp=1) plot(ncf.cor) 

Answer 1

To compare the results of a correlogram, two things should be considered in your case. (i) your code only works for boxes proportional to the resolution of your raster. In this case, a small difference in the bins can make the inclusion or exclusion of an important number of pairs. (ii) An irregular raster shape has a strong effect on pairs that are believed to calculate correlations for a given distance interval. Therefore, your code must deal with both, allow any value for the length of the hopper, and consider the irregular shape of the raster. Below is a small modification of your code to solve these problems.

 # SpatialPack correlation library(SpatialPack) test <- modified.ttest(z.value,z.value,coords = xy,nclass = 21) # Own correlation bins <- test$upper.bounds library(raster) sp.Corr <- matrix(nrow = 0,ncol = 2) for (i in bins) { cat(paste0("..",i)) #print the bin, which is currently calculated w = focalWeight(r,d = i,type = 'circle') wTemp <- w #temporarily saves the weigtht matrix if (i > bins[1]) { midpoint <- ceiling(dim(w)/2) # get the midpoint half_range <- floor(dim(wOld)/2) w[(midpoint[1] - half_range[1]):(midpoint[1] + half_range[1]), (midpoint[2] - half_range[2]):(midpoint[2] + half_range[2])] <- w[(midpoint[1] - half_range[1]):(midpoint[1] + half_range[1]), (midpoint[2] - half_range[2]):(midpoint[2] + half_range[2])]*(wOld==0) w <- w * (1/sum(w)) #normalizes the vector to sum the weights to 1 } wOld <- wTemp #save this weight matrix for the next run mor <- Moran(r,w=w) sp.Corr <- rbind(sp.Corr,c(Moran =mor,Distance = i)) } # Comparing plot(x=test$upper.bounds, test$imoran[,1], col = 2,type = "b",ylab = "Moran I",xlab="Upper bound of distance", lwd = 2) lines(x=sp.Corr[,2],y = sp.Corr[,1], col = 3) points(x=sp.Corr[,2],y = sp.Corr[,1], col = 3) legend('topright', legend = c('SpatialPack', 'Own code'), col = 2:3, lty = 1, lwd = 2:1) 

The image shows that the results of using SpatialPack and native code are the same.