####################### ### Point processes ### ####################### ###################### # library 'spatstat' # ###################### library(spatstat) # Realization of a point process (a point pattern) is represented # as an object of the class 'ppp'. # It is necessary to specify the x- and y-coordinates of the points # and the observation window (by default the unit square). # It is also possible to add marks to the points. x <- runif(100) # vector of 100 uniformly distributed random numbers in [0,1] y <- runif(100) plot(x,y) # plotting the data with the standard 'plot' function (for numeric data) ### binomial point process with 100 points in [0,1]^2 bpp <- ppp(x,y) # point pattern with point coordinates x and y, defaultne observation window plot(bpp) # now the 'plot' function form spatstat is used (plot.ppp) plot(bpp,pch=19) # parameter 'pch' changes the style of the individual points # various parameters for the plotting can be set, e.g.: # chars ... appearance of the point (the value can be a number or a symbol), # - instead of 'chars' one can sometimes use 'pch' # cols ... color of the points (number or name of the color) # main ... title (string) ### adding marks mks <- rexp(100,2) # exponencially distributed marks mbpp <- ppp(x,y,marks=mks) # independently marked binomial point process plot(mbpp) ### the vectors x, y and the marks can be specified manually or loaded from a file x0 <- c(0.32,1.75,1.46,0.84,0.06,1.88) y0 <- c(0.15,2.21,0.92,2.78,1.29,1.52) pp0 <- ppp(x0,y0,window=owin(c(0,2),c(0,3))) # rectangular window (0,2) x (0,3) plot(pp0) ### or one can create a point pattern interactively P <- clickppp(10) # 10 points in a unit square ### this method can be also used to add points to an existing pattern plot(pp0) # left button to add, ESC to stop pp0add <- clickppp(add=T,win=owin(c(0,2),c(0,3))) # Task: try "generating" a Poisson process, i.e. no interactions between points, # complete spatial randomness. We will use the generated pattern later. pois.pattern <- clickppp(win=square(1)) ### list of available datasets data(package="spatstat") ### loading a dataset with the 'data' function (often not necessary in the newer versions) data(swedishpines) plot(swedishpines,chars=8,cols=2) summary(swedishpines) # other examples: bei, cells, gordon, japanesepines, nztrees, osteo (3D), ponderosa, redwood data(amacrine) # two-dimensional point pattern plot(amacrine,chars=c('A','B'),cols=c(3,6)) summary(amacrine) # other examples: ants, chorley, ganglia, hamster, mucosa, paracou, urkiola # multidimensional point patterns: hyytiala, lansing data(spruces) # marked point pattern plot(spruces) points(unmark(spruces),pch=19,col="blue") summary(spruces) # other examples: finpines, longleaf, waka # Object of the type 'ppp' contains: x (vector of x-coordinates) # y (vector of y-coordinates) # n (number of points) # window (observation window) # marks (marks, if present) # these can be accessed using $, e.g.: P$x # we created this one manually amacrine$y pp0add$n swedishpines$window spruces$marks ################################ # simulation of point patterns # ################################ # Another way of obtaining a point pattern is using simulation. ### 1. binomial point process # (homogeneous) binomial point process with 90 points in [0,3]^2 pp1a <- runifpoint(90,square(3)) plot(pp1a,pch=19) # (homogeneous) binomial point process with 50 points in a circle with radius 3 plot(pp1b <- runifdisc(50,radius=3),pch=19) f <- function(x,y) { abs(sin(x+y)) } # probability density function (unnormalized version is OK) xgrid <- seq(0,pi,len=101) # discretization used for plotting ygrid <- seq(0,2,len=101) z <- outer(xgrid,ygrid,"+") image(xgrid,ygrid,abs(sin(z))) # 50 independent points placed in a rectangular window (0,pi) x (0,2) according to f pp1c <- rpoint(50,f,win=owin(c(0,pi),c(0,2))) points(pp1c,pch=19) ### 2. Poisson point process lambda <- 10 pp2a <- rpoispp(lambda,win=pp1a$window) pp2a$n par(mfrow=c(1,2)) # for plotting two graphs next to each other # is there a difference? plot(pp1a,pch=19,main="binomial") plot(pp2a,pch=19,main="Poisson") # is there a difference? plot(pp2a,pch=19,main="Poisson") plot(pois.pattern, pch=19,main="hand-made") par(mfrow=c(1,1)) lambda <- function(x,y,a,b) { a*exp(-b*x) } # intensity function depending on a, b plot(pp2b <- rpoispp(lambda,a=200,b=3)) # we specify a, b ### 3. Matern cluster process a <- 20 b <- 5 pp3 <- rMatClust(kappa=a,r=0.15,mu=b,win=letterR) # The parameters are: # kappa ... intensity of the parent process # r ... radius of the circle for placing (uniformly) offspring points # mu ... mean number of points in a cluster # win ... observation window plot(pp3,pch=19) pp3$n # actual number of points in the realization a*b*area.owin(letterR) # expected number of points ### 4. Thomas process pp4 <- rThomas(5,0.05,10,owin(poly=list(x=c(0,2,0),y=c(0,0,3)))) # The parameters are: intensity of the parent process; standard deviation of the Gaussian distribution; # mean number of points in a cluster; observation window. plot(pp4,pch=19) ### 5. Matern hard-core process type I lambda <- 500 # intensity of a stationary Poisson process r <- 0.05 # hard-core distance plot(pp5 <- rMaternI(lambda,r,amacrine$window)) # observation window taken from another pattern # compare the mean number of points with the actual number of points: c(lambda*exp(-lambda*pi*r^2)*area.owin(pp5$window),pp5$n) ### 6. Matern hard-core process type II plot(pp6 <- rMaternII(lambda,r,pp5$window)) # compare the mean number of points with the actual number of points: c((1-exp(-lambda*pi*r^2))*area.owin(pp6$window)/(pi*r^2),pp6$n) # Task: generate several realizations of a Poisson process, # cluster process, regular process; plot them together # and describe the differences. ################################# # Point process characteristics # ################################# # estimate of the intensity function plot(density.ppp(pp2b)) points(pp2b, pch=20) # estimate of the pair-correlation function data(japanesepines) data(redwood) data(cells) plot(pcf(pp2a)) plot(pcf(japanesepines)) plot(pcf(redwood)) plot(pcf(cells)) # estimate of the K-function plot(Kest(pp2a)) par(mfrow=c(1,2)) plot(Kest(pois.pattern), legend=FALSE) plot(Kest(pp2a), legend=FALSE) par(mfrow=c(1,1)) # estimate of the nearest-neighbour distribution function plot(Gest(pp2a)) # estimate of the contact distribution function plot(Fest(pp2a)) # estimate of the J-function plot(Jest(pp2a)) # estimates of the four functions (F, G, J, K) using 'allstats' plot(allstats(pp2a), legend=FALSE) plot(allstats(japanesepines), legend=FALSE) plot(allstats(redwood), legend=FALSE) plot(allstats(cells), legend=FALSE) ################################ # operations on point patterns # ################################ # thinning plot(pp2a) # keep points with specified indices: points(pp2a[1:10],col="red",pch=19) # keep points fulfilling a condition: points(pp2a[pp2a$x+pp2a$y<=2.5],col="darkblue",pch=19) # independent random thinning pp2a_thin <- rthin(pp2a,0.7) plot(pp2a) # keep each point with probability 0.7 points(pp2a_thin,col="green",pch=19) f <- function(x,y) { 1-(x^2+y^2)/18 } pp2a_thin2 <- rthin(pp2a,f) plot(pp2a) # keep each point with probability depending on the distance from the bottom-left corner points(pp2a_thin2,pch=19) # restriction to a subwindow subwindow <- owin(c(1.2,2.5),c(0.4,1.7)) # [1.2,2.5] x [0.4,1.7] par(mfrow=c(1,2)) plot(pp2a,main="Window [0,3] x [0,3]") plot(subwindow,add=TRUE) plot(pp2a[subwindow],main="Window [1.2,2.5] x [0.4,1.7]") # superposition: ppsuperposition <- superimpose(pp1a,pp2a) par(mfrow=c(1,3)) plot(pp1a) plot(pp2a) plot(ppsuperposition) par(mfrow=c(1,1))