In AdamWilsonLab/rasterAutocorr: Quickly calculate spatial autocorrelation on 2D rasters.

Load necessary libraries

library(raster)
library(fields)
library(rasterAutocorr)
library(dplyr)
library(knitr)
library(ggplot2)

Create a random raster with some spatial structure

nx=301  # size of raster
ny=round(nx*1.2)

r=raster(nrows=ny, ncols=nx,vals=1,xmn=-nx/2, xmx=nx/2, ymn=-ny/2, ymx=ny/2)
names(r)="z"

Simulate a Gaussian random field with an exponential covariance function

Theta is the scale of an exponential decay function. This controls degree of autocorrelation, values close to 1 are close to random while values near nx/4 have high autocorrelation

theta=50
grid=list(x=seq(xmin(r),xmax(r)-1,by=res(r)[1]),y=seq(ymin(r),ymax(r)-1,res(r)[2]))

obj<-Exp.image.cov(grid=grid, theta=theta, setup=TRUE)
look<- sim.rf( obj)

values(r)=t(look)

Plot the simulated raster

plot(r,ylab="Y",xlab="X",main="Original (Simulated) Raster")

Make the correlogram using fft()

Fit the complete spatial autocorrelation using fft()

a1=acorr(r)

Extract directions for each shift to facilitate plotting of the correlogram

d2=acorr_dir(r)

plot the autocorrlation and distance

  plot(a1[[1]],ylab="Shift in Y",xlab="Shift in X",main="Autocorrelation")
  plot(d2,ylab="Shift in Y",xlab="Shift in X",main="Direction from center (degrees)")

Combine the acorr values with distance

ftd=data.frame(cor=values(a1))#,dir=values(d2))

## round to  faciliate binning
ftd$dist=round(ftd$cor.dist,3)#,c(0:50,seq(51,1000,by=10)))
## take mean by km
ftd2 <- group_by(ftd, dist)
ftd2 <- summarise(ftd2,
                  min = min(cor.acor, na.rm = TRUE),
                  max = max(cor.acor, na.rm = TRUE),
                  sd = sd(cor.acor, na.rm = TRUE),
                  mean = mean(cor.acor, na.rm = TRUE)
)

This table now has the correlation values and distance for each pixel

kable(head(ftd2))

plot the correlogram

ggplot(data=ftd2,aes(y=mean, x=dist))+
  geom_linerange(aes(ymax=max,ymin=min),col="grey")+
  geom_line()+
  ylab("Correlation")+xlab("Distance") +
  geom_hline(aes(yintercept = 0))

Matlab/Octave script

x1=[3 6 5 ; 7 2 2 ; 4 NaN 0]
[n,p]=size(x1); 
ncols=2*p-1;
nrows=2*n-1;
nr2=5
nc2=5
x1id=~isnan(x1);
x1(~x1id)=zeros(sum(sum(~x1id),1);
fx1=fft2(x1,nr2,nc2);
fx1id=fft2(x1id,nr2,nc2);
nh11=round(real(ifft2(conj(fx1id) .*fx1id)));
m1=real(ifft2(conj(fx1) .*fx1id)) ./max(nh11,1);
m2=real(ifft2(conj(fx1id) .*fx1)) ./max(nh11,1)
gh11=real(ifft2(conj(fx1) .*fx1)); gh11=gh11./max(nh11,1)-m1.*m2;

AdamWilsonLab/rasterAutocorr documentation built on March 19, 2020, 12:01 a.m.