WS 2016/2017
Prerequisites
setwd, library, data
setwd("YOUR-PATH")
raw_data <- read.table(file = "data/FS_UTM.csv",
header = TRUE,
sep = ";",
stringsAsFactors = FALSE,
quote = "")
Create a spatial object
Libraries and creation of sp object
Version 1
library(sp)
library(rgdal)
sites <- SpatialPointsDataFrame(coords = cbind(raw_data$xUTM,raw_data$yUTM),
data = data.frame(raw_data[,2]),
proj4string = CRS("+init=epsg:32634")
)
class(sites)
## [1] "SpatialPointsDataFrame" ## attr(,"package") ## [1] "sp"
Libraries and creation of sp object
Version 2
sites <- raw_data
coordinates(sites) <- ~xUTM+yUTM
proj4string(sites) <- CRS("+init=epsg:32634")
Plot it 2.0
Let's plot with leaflet. Since OpenStreetMap, google maps, etc. are global, we need a geographic coordinate system like WGS84
is.projected(sites)
## [1] TRUE
sites.wgs84 <- spTransform(sites, "+init=epsg:4326") is.projected(sites.wgs84)
## [1] FALSE
What's this?!
Plot it 2.0
library(magrittr)
library(leaflet)
m1 <- leaflet() %>%
addProviderTiles("Thunderforest.Landscape") %>%
addMarkers(lng=sites.wgs84@coords[,1],
lat=sites.wgs84@coords[,2],
popup = sites.wgs84@data$Site.Name
)
m1
Export to Shapefile
writeOGR(obj = sites,
dsn = "./data",
layer = "A_Sites",
driver = "ESRI Shapefile",
overwrite_layer = TRUE
)
## Warning in writeOGR(obj = sites, dsn = "./data", layer = "A_Sites", driver ## = "ESRI Shapefile", : Field names abbreviated for ESRI Shapefile driver
list.files("./data")
## [1] "A_Sites.dbf" "A_Sites.prj" "A_Sites.shp" "A_Sites.shx" ## [5] "FS_UTM.csv" "FS_UTM.xlsx" "srtm_41_03.tif"
Spatial Point Patterns
Infos
THE references for this subject are: Baddeley et al. (2016), Baddeley (2008), Bivand et al. (2008), Diggle (2013)
And the best thing is: the authors are active developers of R. The package we will use in this topic very often is spatstat. A package so full of functions that the manual alone has over 1589 pages (as of 2016-02-03).
If you use the package in your research, do not forget to cite it (this holds true for every package used, including the base software R !):
library(spatstat) citation(package = "spatstat")
## ## To cite spatstat in publications use: ## ## Adrian Baddeley, Ege Rubak, Rolf Turner (2015). Spatial Point ## Patterns: Methodology and Applications with R. London: Chapman ## and Hall/CRC Press, 2015. URL ## http://www.crcpress.com/Spatial-Point-Patterns-Methodology-and-Applications-with-R/Baddeley-Rubak-Turner/9781482210200/ ## ## If you use hybrid models, please also cite: ## ## Adrian Baddeley, Rolf Turner, Jorge Mateu, Andrew Bevan (2013). ## Hybrids of Gibbs Point Process Models and Their Implementation. ## Journal of Statistical Software, 55(11), 1-43. URL ## http://www.jstatsoft.org/v55/i11/. ## ## In survey articles, please cite the original paper on spatstat: ## ## Adrian Baddeley, Rolf Turner (2005). spatstat: An R Package for ## Analyzing Spatial Point Patterns. Journal of Statistical ## Software 12(6), 1-42. URL http://www.jstatsoft.org/v12/i06/.
Infos 2
Some resources for spatstat package:
Create ppp object
library(spatstat)
sites.window <- owin(xrange = c(min(sites@coords[,1]),max(sites@coords[,1])),
yrange = c(min(sites@coords[,2]),max(sites@coords[,2]))
)
sites.pp <- ppp(x = sites@coords[,1],
y = sites@coords[,2],
window = sites.window
)
## Warning in ppp(x = sites@coords[, 1], y = sites@coords[, 2], window = ## sites.window): data contain duplicated points
unitname(sites.pp) <- c("meter", "meters")
inspect and plot
str(sites.pp)
## List of 5 ## $ window :List of 4 ## ..$ type : chr "rectangle" ## ..$ xrange: num [1:2] 443157 696380 ## ..$ yrange: num [1:2] 5012062 5145452 ## ..$ units :List of 3 ## .. ..$ singular : chr "meter" ## .. ..$ plural : chr "meters" ## .. ..$ multiplier: num 1 ## .. ..- attr(*, "class")= chr "units" ## ..- attr(*, "class")= chr "owin" ## $ n : int 508 ## $ x : num [1:508] 443157 459882 462115 466751 467930 ... ## $ y : num [1:508] 5105906 5095991 5088415 5089685 5109462 ... ## $ markformat: chr "none" ## - attr(*, "class")= chr "ppp"
plot(sites.pp)
Measures of Central tendency
Mean center (mc) and standard distance (stdist)
mc <- cbind(sum(sites@coords[,1])/length(sites@coords[,1]),
sum(sites@coords[,2])/length(sites@coords[,2])
)
stdist <- sqrt(sum((sites@coords[,1]-mean(sites@coords[,1]))^2 +
(sites@coords[,2]-mean(sites@coords[,2]))^2) /
length(sites@coords[,1])
)
Measures of Central tendency
plot(sites)
points(mc, col ="red")
library(plotrix)
draw.circle(x = mc[1],
y = mc[2],
radius = stdist,
border = "red"
)
First-order effects
Global intensity
Global intensity = Number of points per area
## A = a * b area.sqm <- diff(sites.pp$window$xrange) * diff(sites.pp$window$yrange) area.sqkm <- area.sqm*10^-6 # area <- area/1000000 area.sqkm
## [1] 33777.36
## calculate intensity intensity <- sites.pp$n/area.sqkm intensity
## [1] 0.01503966
Local intensity
qc.sites <- quadratcount(X = sites.pp) plot(qc.sites) points(sites.pp, pch = 20, cex = .5, col = rgb(.2,.2,.2,.5))
Local intensity
The question is: do the quadratcounts indicate CSR? To check we use a \(\chi^2\) test approach
(relation between observed (i.e. empirical) and expected (i.e. theoretical, here CSR) amounts of points in quadrants)
qt.sites <- quadrat.test(X = sites.pp) qt.sites
## ## Chi-squared test of CSR using quadrat counts ## Pearson X2 statistic ## ## data: sites.pp ## X2 = 610.11, df = 24, p-value < 2.2e-16 ## alternative hypothesis: two.sided ## ## Quadrats: 5 by 5 grid of tiles
Local intensity
observed (lop left) | expected (top right) | Pearson residual (bottom)
+/- 2 = unusual; > +/- 2 = gross departure from fitted model
plot(qt.sites)
Density calculation - 1st approach
Amount of events within pixel…in imprecise language: a histogram in space. This is a step by step example but the code is rather inefficient.
Step 1 - data preparation
cs <- 10000 # cellsize # enlarge the study area by ONE pixel (in E-W and N-S direction) xmin <- sites.pp$window$xrange[1] - cs/2 xmax <- sites.pp$window$xrange[2] + cs/2 ymin <- sites.pp$window$yrange[1] - cs/2 ymax <- sites.pp$window$yrange[2] + cs/2 ## calculate the number of rows/columns ## (add 1 just because a pixel might get "lost" through the rounding operation) rows <- round((ymax-ymin)/cs, 0) + 1 columns <- round((xmax-xmin)/cs, 0) + 1 z <- cbind(1:(columns*rows)) # create a vector with all grid cells df <- data.frame(z) # create a data.frame of it
Density calculation - 1st approach
Step 2 - create grid topology
# create a topological description of a grid
gt <- GridTopology(cellcentre.offset = c(sites.pp$window$xrange[1] - cs/2,
sites.pp$window$yrange[1] - cs/2),
cellsize = c(cs,cs),
cells.dim = c(columns,rows))
gt
## X1 X2 ## cellcentre.offset 438157 5007062 ## cellsize 10000 10000 ## cells.dim 27 15
Density calculation - 1st approach
Step 3 - create grid
sgdf <- SpatialGridDataFrame(grid = gt,
data = df,
proj4string = CRS("+init=epsg:32634")
)
plot(sgdf)
Density calculation - 1st approach
Step 4 - calculate density
for (i in seq(along=coordinates(gt)[,1])){
## because the coordinate is defined for the
## center of the cell, the half cellsize
## is substracted
x <- coordinates(gt)[i,1] - cs/2
y <- coordinates(gt)[i,2] - cs/2
## which events lie within the x/y direction?
xi <- which(sites.pp$x > x &
sites.pp$x < x + cs)
yi <- which(sites.pp$y > y &
sites.pp$y < y + cs)
## how many objects in x and y direction intersect?
pz <- length(intersect(xi,yi))
## divide the number of points by the area
sgdf@data$z[i]<- pz / (cs/1000)^2
}
plot(raster::raster(sgdf))
points(sites.pp$x, sites.pp$y, pch=16, cex=0.4)
Density calculation - 2nd approach
Kernel density estimation; again here performed step by step and rather inefficiently/partially incorrect coded. Use it to understand the concept.
- for all cells
- get cell center
- for all points
- calculate distance to cell center; Euclidean distance -> \(\sqrt{(x_i-x)^2+(y_i-y)^2}\)
- weight the distance according to kernel function
- write value in pixel
Density calculation - 2nd approach
sgdf_kde <- sgdf
## kernel bandwidth
sd <- 50000
for (i in seq(along=coordinates(gt)[,1])){
x <- coordinates(gt)[i,1]
y <- coordinates(gt)[i,2]
g2 <- 0
for (j in seq(along=sites.pp$x)){
distance <- sqrt((sites.pp$x[j] - x)^2 +
(sites.pp$y[j] - y)^2
)
g1 <- dnorm(distance, mean=0, sd=sd)
g2 <- g2 + g1
}
sgdf_kde@data$z[i]<- g2
}
library(raster)
plot(raster(sgdf_kde))
points(sites.pp$x, sites.pp$y, pch=16, cex=0.4)
Density calculation - 3rd approach
Kernel density estimation using the built in function of spatstat. Since the package is great, it is perfectly coded and fast.
cs <- 10000
sd <- 50000
sites.dens <- density(x = sites.pp,
bw = sd,
eps=cs,
edge=FALSE,
at="pixels"
)
plot(sites.dens, col = grey(seq(0.2,1,.1)))
contour(sites.dens, add=T)
points(sites.pp$x, sites.pp$y, pch=16, cex=0.4)
Density calculation - 3rd approach
Use another measure for sd - in this case three times the mean nearest neighbor distance.
sdev <- 3*mean(nndist(sites.pp))
sites.dens2 <- density(x = sites.pp,
bw = sdev,
eps=cs,
edge=TRUE,
at="pixels"
)
plot(sites.dens2, col = grey(seq(0.2,1,.1)))
contour(sites.dens2, add=T)
points(sites.pp$x, sites.pp$y, pch=16, cex=0.4)
Second-order effects
Nearest-neighbor distance
sites.nn <- nndist(X = sites.pp) str(sites.nn)
## num [1:508] 19443 7899 3519 132 3026 ...
mean(sites.nn)
## [1] 1623.151
hist(sites.nn) abline(v=mean(sites.nn)) abline(v=median(sites.nn), lty=2)
Nearest neighbor Distance – Clark and Evans’ R
\(R = \bar{d_{min}} / \frac{1}{2\sqrt{\lambda}}\)
"An \(R\) value of less than 1 indicates of a tendency toward clustering, since it shows that observed nearest neighbor distances are shorter than expected. An R value of more than 1 indicatives of a tendency toward evenly spaced events" (O'Sullivan & Unwin 2010, 144)
nnE <- 1/(2*sqrt((sites.pp$n/area.sqm))) nnE
## [1] 4077.097
R.sites <- mean(sites.nn)/nnE R.sites
## [1] 0.3981143
G function
First approach based on theoretical assumptions, i.e. CSR
sites.g <- Gest(sites.pp) plot(sites.g)
G function
Second approach based on simulations of a theortical process, i.e. CSR
sites.g.env <- envelope(sites.pp,
fun = Gest,
nsim = 10)
## Generating 10 simulations of CSR ... ## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. ## ## Done.
plot(sites.g.env)
F function
First approach based on theoretical assumptions, i.e. CSR
sites.f <- Fest(sites.pp) plot(sites.f)
F function
Second approach based on simulations of a theortical process, i.e. CSR
sites.f.env <- envelope(sites.pp,
fun = Fest,
nsim = 10)
## Generating 10 simulations of CSR ... ## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. ## ## Done.
plot(sites.f.env)
K function
First approach based on theoretical assumptions, i.e. CSR
sites.k <- Kest(sites.pp) plot(sites.k)
K function
Second approach based on simulations of a theortical process, i.e. CSR
sites.k.env <- envelope(sites.pp,
fun = Kest,
nsim = 10)
## Generating 10 simulations of CSR ... ## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. ## ## Done.
plot(sites.k.env)
L function
First approach based on theoretical assumptions, i.e. CSR
sites.l <- Lest(sites.pp) plot(sites.l)
L function
Second approach based on simulations of a theortical process, i.e. CSR
sites.l.env <- envelope(sites.pp,
fun = Lest,
nsim = 10)
## Generating 10 simulations of CSR ... ## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. ## ## Done.
plot(sites.l.env)
From global to local measure
G and moving window
Where are the points clustered/ordered/…
Procedure
- take a raster
- choose a moving window
- define sample radius for G-function
- check [cell by cell] characteristics of G
Preparation
sgdf.g <- sgdf # grid #sites.nn # nearest neighbors #sites.window # window radius <- 25000 # Radius of moving window r <- seq(0, radius, 50) # Disctance vector for G; it is not advised to use it!
G and moving window
for (i in seq(along=sgdf.g@data$z)) {
xr <- coordinates(sgdf.g)[i,1]
yr <- coordinates(sgdf.g)[i,2]
distances <- sqrt((sites.pp$x -xr)^2 + (sites.pp$y -yr)^2)
## Which points are within the moving window?
index <- which(distances<radius)
## at least three points need to be in the moving window
if (length(index) > 2 ) {
x <- sites.pp$x[index]
y <- sites.pp$y[index]
name <- sites$Site.Name[index]
sites2 <- SpatialPointsDataFrame(cbind(x, y),
as.data.frame(name),
proj4string= CRS("+init=epsg:32634")
)
## Point pattern from points in moving window
sites.ppt <- ppp(sites2@coords[,1],sites2@coords[,2],window=sites.window)
g <- Gest(sites.ppt, r=r)
## mean differnce "theoretical - empirical"; positiv = orered
meandiff <- mean(g$theo-g$km)
sgdf.g@data$z[i] <- meandiff
}
else {sgdf.g@data$z[i] <-NA}
}
## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points
G and moving window
library(raster) plot(raster(sgdf.g)) points(sites, cex = .2, pch = 19)
F and moving window
Where are the points clustered/ordered/…
Procedure
- take a raster
- choose a moving window
- define sample radius for G-function
- check [cell by cell] characteristics of F
Preparation
sgdf.f <- sgdf # grid #sites.nn # nearest neighbors #sites.window # window radius <- 25000 # Radius of moving window r <- seq(0, radius, 50) # Disctance vector for G; it is not advised to use it!
F and moving window
for (i in seq(along=sgdf.f@data$z)) {
xr <- coordinates(sgdf.f)[i,1]
yr <- coordinates(sgdf.f)[i,2]
distances <- sqrt((sites.pp$x -xr)^2 + (sites.pp$y -yr)^2)
## Which points are within the moving window?
index <- which(distances<radius)
## at least three points need to be in the moving window
if (length(index) > 2 ) {
x <- sites.pp$x[index]
y <- sites.pp$y[index]
name <- sites$Site.Name[index]
sites2 <- SpatialPointsDataFrame(cbind(x, y),
as.data.frame(name),
proj4string= CRS("+init=epsg:32634")
)
## Point pattern from points in moving window
sites.ppt <- ppp(sites2@coords[,1],sites2@coords[,2],window=sites.window)
f <- Fest(sites.ppt, r=r)
## mean differnce "theoretical - empirical"; positiv = orered
meandiff <- mean(f$theo-f$km)
sgdf.f@data$z[i] <- meandiff
}
else {sgdf.f@data$z[i] <-NA}
}
## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points ## Warning in ppp(sites2@coords[, 1], sites2@coords[, 2], window = ## sites.window): data contain duplicated points
F and moving window
library(raster) plot(raster(sgdf.f)) points(sites, cex = .2, pch = 19)
G and F moving window
References
Baddeley, A., 2008. Analysing spatial point patterns in R, CSIRO; University of Western Australia. Available at: http://www.csiro.au/files/files/pn0y.pdf.
Baddeley, A., Rubak, E. & Turner, R., 2016. Spatial point patterns: Methodology and applications with r, CRC Press, Taylor & Francis Group.
Bivand, R.S., Pebesma, E.J. & Gómez-Rubio, V., 2008. Applied Spatial Data Analysis with R, New York: Springer.
Diggle, P.J., 2013. Statistical Analysis of Spatial and Spatio-Temporal Point Patterns, Third Edition 3rd ed., Boca Raton: Chapman; Hall/CRC.