使用栅格包的空间相关图

Dear Crowd

Problem

I tried to calculate a spatial correlogram with the packages nfc, pgirmess, SpatialPack and spdep. However, I was troubling to define the start and end-point of the distance. I'm only interested in the spatial autocorrelation at smaller distances, but there on smaller bins. Additionally, as the raster is quite large (1.8 Megapixels), I run into memory troubles with these packages but the SpatialPack.

我试图用包nfc,pgirmess,SpatialPack和spdep来计算空间相关图。但是,我很难定义距离的起点和终点。我只对较小距离的空间自相关感兴趣,但是在较小的箱子上。此外,由于光栅非常大(1.8百万像素),我遇到了这些软件包的内存问题,但是SpatialPack。

So I tried to produce my own code, using the function Moran from the package raster. But I must have some error, as the result for the complete dataset is somewhat different than the one from the other packages. If there is no error in my code, it might at least help others with similar problems.

所以我尝试使用包栅格中的Moran函数生成我自己的代码。但我必须有一些错误,因为完整数据集的结果与其他包中的结果略有不同。如果我的代码中没有错误,它至少可以帮助其他类似问题的人。

Question

I'm not sure, whether my focal matrix is erroneous. Could you please tell me whether the central pixel needs to be incorporated? Using the testdata I can't show the differences between the methods, but on my complete dataset, there are differences visible, as shown in the Image below. However, the bins are not exactly the same (50m vs. 69m), so this might explain parts of the differences. However, at the first bin, this explanation seems not to be plausible to me. Or might the irregular shape of my raster, and different ways to handle NA's cause the difference?

我不确定我的焦点矩阵是否错误。你能否告诉我中心像素是否需要合并?使用testdata我无法显示方法之间的差异,但在我的完整数据集上,可见差异,如下图所示。然而,箱子不完全相同(50米对69米),所以这可以解释部分差异。然而,在第一个箱子里,这个解释对我来说似乎不合理。或者我的光栅的不规则形状,以及处理NA的不同方式会导致差异吗?

Comparison of Own method with the one from SpatialPack

自己的方法与SpatialPack的比较

Runable Example

Testdata

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

计算测试数据的代码取自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)

Own 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's I",xlab="Upper bound of distance")

Other methods to calculate the 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's I",xlab="Upper bound of distance")

library(ncf)
ncf.cor <- correlog(x.coord, y.coord, z.value,increment=10, resamp=1)
plot(ncf.cor)

1 个解决方案

#1

In order to compare the results of the correlogram, in your case, two things should be considered. (i) your code only works for bins proportional to the resolution of your raster. In that case, a bit of difference in the bins could make to include or exclude an important amount of pairs. (ii) The irregular shape of the raster has a strong impact of the pairs that are considered to compute the correlation for certain distance interval. So your code should deal with both, allow any value for the length of bin and consider the irregular shape of the raster. A small modification of your code to tackle those problems are below.

为了比较相关图的结果,在您的情况下,应该考虑两件事。 (i)您的代码仅适用于与光栅分辨率成比例的分档。在这种情况下,箱中的一些差异可以包括或排除重要数量的对。 (ii)光栅的不规则形状具有对的强烈影响,被认为是计算某个距离间隔的相关性。因此,您的代码应该处理两者,允许bin的长度的任何值并考虑栅格的不规则形状。下面将对您的代码进行一些小修改以解决这些问题。

# 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's 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 the SpatialPack package and the own code are the same.

该图显示使用SpatialPack包和自己的代码的结果是相同的。