##################### ### Random fields ### ##################### install.packages(c("spdep","maptools","RANN","RandomFields","geoR")) ############### # Ising model # ############### # Wolfram Demonstration: # http://msekce.karlin.mff.cuni.cz/~dvorak/teaching/2015_2016/teaching_2015_2016.html ################### # library 'spdep' # ################### library(spdep) # examples of regional data data(package="spdep") library(maptools) # crime data from 49 regions ('neighbourhoods') in Columbus (Ohio) # https://en.wikipedia.org/wiki/Columbus,_Ohio example(columbus) plot(columbus) coords <- coordinates(columbus) points(coords,pch=19) ### neigbours col.gal.nb # neigbourhood relation plot(col.gal.nb,coords,add=TRUE) summary(col.gal.nb) is.symmetric.nb(col.gal.nb) print(cards <- card(col.gal.nb)) # counts of neigbours ### weights for assessing the spatial autocorrelation col.b <- nb2listw(col.gal.nb,style="B") # binary weights col.b$weights col.w <- nb2listw(col.gal.nb,style="W") # normalized binary weights (sum up to 1 for each point) col.w$weights ### second-order neigbours col.lags <- nblag(col.gal.nb, 2) # neighbours of neighbours plot(col.lags[[2]], coords, add=TRUE, col="red", lty=2) summary(col.lags[[2]], coords) ### k nearest neighbours (in Euclidean space, not necessarily in the neigbourhood relation) library(RANN) col.k4.nb <- knn2nb(knearneigh(coords, 4)) # four nearest neighbours is.symmetric.nb(col.k4.nb) plot(columbus,border="grey") plot(col.k4.nb,coords,arrows=TRUE,add=TRUE) ### neighbourhood relation based on the Euclidean distance col.d.nb <- dnearneigh(coords,0,0.5) # neighbours if closer than 0.5 is.symmetric.nb(col.d.nb) plot(columbus,border="grey") plot(col.d.nb,coords,add=TRUE) ### continuous data # residential burglaries and vehicle thefts per thousand households columbus$CRIME summary(columbus$CRIME) # Tests for "no spatial autocorrelation" hypothesis: # We assume a constant mean value and constant variance! # Moran index for normalized weights, test based on the randomization assumption # (all permutations of assigning the observed values to the lattice points are equally probable) moran.test(columbus$CRIME,listw=col.w) # Moran index for normalized weights, test based on the normality assumption # (uncorrelated random variables with Gaussian distribution) moran.test(columbus$CRIME,listw=col.w,randomisation=FALSE) # Moran index for binary weights: moran.test(columbus$CRIME,listw=col.b) # Moran index for normalized weights (four nearest neighbours): moran.test(columbus$CRIME,listw=nb2listw(col.k4.nb)) # Geary index for normalized weights, test based on the randomization assumption: geary.test(columbus$CRIME,listw=col.w) ### binary data HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high")) plot(columbus) points(coords[HICRIME=="high",],col="red",pch=19) points(coords[HICRIME=="low",],col="blue",pch=19) joincount.test(HICRIME, listw=col.w) # BB statistic joincount.test(HICRIME, listw=col.b) joincount.test(HICRIME, listw=nb2listw(col.k4.nb)) # the tests are performed under the assumption of 'nonfree sampling', i.e. the hypergeometric sampling; # we have a fixed number of red and blue regions ### neigbours in a regular lattice nb7rt <- cell2nb(7, 7) nb7rt summary(nb7rt) xyc <- attr(nb7rt, "region.id") # locations of the sites (lattice points) xy <- matrix(as.integer(unlist(strsplit(xyc, ":"))), ncol=2, byrow=TRUE) plot(nb7rt, xy) print(card(nb7rt)) # counts of neighbours nb7rtt <- cell2nb(7, 7, torus=TRUE) # toroidal situation, no edge-effects summary(nb7rtt) # every point has now four neighbours print(card(nb7rtt)) ### uncorrelated data on the regular lattice data7rt <- runif(49,0,1) moran.test(data7rt,nb2listw(nb7rt)) geary.test(data7rt,nb2listw(nb7rt)) ########################## # library 'RandomFields' # ########################## library(RandomFields) ### R >= 3.0.2 ### simulation of a Gaussian random field # spherical autocovariance function: model1 <- RMspheric(var=2,scale=0.2) plot(model1,xlim=c(-0.5,0.5),ylim=c(0,2)) # two-dimensional regular lattice of points in [0,1]^2 x <- seq(0,1,by=0.02); y <- seq(0,1,by=0.02) # realization of a stationary Gaussian random field Z1 <- RFsimulate(model1,x,y,grid=TRUE) plot(Z1) # we can add some trend (non-constant mean value) and nugget (residual variance) # in the following we use the exponential autocovariance function model2 <- RMexp(var=1,scale=0.1)+RMnugget(var=0.2)+RMtrend(mean=0.5) plot(model2,ylim=c(0,1.2)) plot(model2,ylim=c(0,1.2), fct.type="Variogram") # semivariogram: 1/2 * var( Z(x+h) - Z(x) ) plot(model2,xlim=c(-0.3,0.3), fct.type="Variogram", dim=2) # semivariogram in 2D # simulation of a three-dimensional Gaussian random field z <- seq(0,1,by=0.02) Z2 <- RFsimulate(model2,x,y,z,grid=TRUE) plot(Z2,MARGIN=c(2,3),MARGIN.slices=1,n.slices=4,col=terrain.colors(64)) # Task: try different values of the parameters and # watch their effect on the resulting realizations # fractional Brownian motion, Hurst parameter H=alpha/2 (H=1/2 -> Brownian motion) model3 <- RMfbm(alpha=1) plot(model3,ylim=c(0,1)) # semivariogram is not a bounded function Z3 <- RFsimulate(model3,x=seq(0,10,by=0.02),grid=TRUE) plot(Z3) # the bigger the value of nu, the smoother the realization # (nu=0.5 corresponds to the exponential autocovariance function) model4 <- RMmatern(nu=2) Z4 <- RFsimulate(model4,x=cbind(x,y),grid=TRUE) plot(Z4) ### empirical variogram ev4 <- RFempiricalvariogram(data=Z4) plot(ev4) # shows also the number of pairs in the given distance plot(ev4, model=model4) ################## # library 'geoR' # ################## library(geoR) # library for a geostatistical analysis ### Kattegat basin: https://en.wikipedia.org/wiki/Kattegat data(kattegat) # Kattegat basin salinity data plot(c(550,770),c(6150,6420),type="n",xlab="X UTM",ylab="Y UTM") points(kattegat, pch=19, add=TRUE) lapply(kattegat$dk, lines, lwd=2) ### River Meuse: https://en.wikipedia.org/wiki/Meuse data(meuse.riv) # data from the library 'sp' plot(meuse.riv,type="l",asp=1) # river Meuse data(meuse) # topsoil heavy metal concentrations in a flood plain points(meuse) summary(meuse) plot(y~x,meuse) plot(meuse[,1:4]) # spatial coordinates, cadmium and copper meuse_zn <- as.geodata(meuse,data.col=6) # zinc plot(c(178000,182000),c(33000,33400),type="n") points(meuse_zn, pch=19) ### Estimation of the semivariogram # nonparametric estimate of the semivariogram for zinc: est_gamma <- variog(meuse_zn,max.dist=1500) plot(est_gamma,type="b") # kernel estimate: smooth <- variog(meuse_zn,option="smooth",max.dist=1500,n.points=100,kernel="normal",band=250) lines(smooth,type="l",lty=2) # interactive fitting of a parametric model (performed by the user!) eyefit(est_gamma) ### Anisotropy est_gamma_45 <- variog(meuse_zn,max.dist=1500,direction=pi/4) plot(est_gamma_45,type="b")