Summary Functions and Kernel Density Estimation
WEN, Tzai-Hung (NTU Geography)
Loading Packages and Reading Sample Data
setwd("D:/_Lab_Data/R/R_Labs")
library(sp);library(rgdal)
library (splancs) #Objects: SpatialPoints
data <- read.table("pts_data/tpedata.csv", header=TRUE, sep=",")
Pts <- as.points(data[,2], data[,3])
data_bnd <- read.table("pts_data/tpe_sqr_bnd2.csv", header=TRUE, sep=",")
Pts_bnd2 <- as.points(data_bnd[,2], data_bnd[,3])
### Sample data for Bivariate K functions
data1 <- read.table("pts_data/point1.csv", header=TRUE, sep=",")
Pts1 <- as.points(data1[,1], data1[,2])
data2 <- read.table("pts_data/point2.csv", header=TRUE, sep=",")
Pts2 <- as.points(data2[,1], data2[,2])
### Plotting sample points and the boundary
polymap(Pts_bnd2, col="gray")
pointmap(Pts, xlab="X", ylab="Y", pch=21, bg="yellow", add=T)
G function is defined as the cumulative frequency distribution of the nearest-neighbor distances
- The shape of G-function tells us the way the events are spaced in a point pattern
- Clustered：G increases rapidly at short distance
- Uniform：G increases slowly up to distance where most events spaced, then increases rapidly
tpe.s<-seq(1,5000,100)
tpe.G<-Ghat(Pts, tpe.s)
plot(tpe.s,tpe.G, type="l", xlab="distance", ylab="Estimated G", main="G function")
# Generating the envelopes
nsim=29
U.Ghat <- -99999
L.Ghat <- 99999
for(i in 1:nsim) {
tpe.CSR<-csr(Pts_bnd2, npts(Pts))
tpe.G.CSR<-Ghat(tpe.CSR, tpe.s)
U.Ghat<-pmax(tpe.G.CSR, U.Ghat)
L.Ghat<-pmin(tpe.G.CSR, L.Ghat)
}
lines(tpe.s, U.Ghat, col="red", lty=3)
lines(tpe.s, L.Ghat, col="blue", lty=3)
F function is similar to the G function, but instead of the events a sample of point locations is selected randomly from anywhere in the study area
- Clustered：F(r) rises slowly at first, but more rapidly at longer distances
- Uniform：F(r) rises rapidly at first, then slowly at longer distances
- Examine significance by simulating “envelopes”
### F Function
tpe.CSR<-csr(Pts_bnd2, npts(Pts))
tpe.F<-Fhat(Pts,tpe.CSR,tpe.s)
plot(tpe.s,tpe.F, type="l", xlab="distance", ylab="Estimated F", main="F function")
# Generating the envelopes
nsim=29
U.Fhat <- -99999
L.Fhat <- 99999
for(i in 1:nsim) {
tpe.CSR2<-csr(Pts_bnd2, npts(Pts))
tpe.F.CSR<-Fhat(tpe.CSR2, tpe.CSR, tpe.s)
U.Fhat<-pmax(tpe.F.CSR, U.Fhat)
L.Fhat<-pmin(tpe.F.CSR, L.Fhat)
}
lines(tpe.s, U.Fhat, col="red", lty=3)
lines(tpe.s, L.Fhat, col="blue", lty=3)
Ripley’s K function is a statistical method for point pattern analysis.
- Summary of local dependence of spatial process: the second order process
- Expresses the number of expected events within given distance of randomly chosen event
- K(t) describes characteristics of the point processes at many distance scales
### Ripley's K-Function
tpe.s<-seq(1,10000,500)
tpe.K <-khat(Pts,Pts_bnd2,tpe.s); tpe.K## [1] 276149 2492473 8422971 15801829 22874884 34986132 49794730
## [8] 65949187 78874014 100111057 121415131 141850788 163056393 187865833
## [15] 210429907 229850882 250841512 273074367 294308303 320173954tpe.env<-Kenv.csr(npts(Pts), Pts_bnd2, nsim=29, tpe.s)## Doing simulation 1
## Doing simulation 2
## Doing simulation 3
## Doing simulation 4
## Doing simulation 5
## Doing simulation 6
## Doing simulation 7
## Doing simulation 8
## Doing simulation 9
## Doing simulation 10
## Doing simulation 11
## Doing simulation 12
## Doing simulation 13
## Doing simulation 14
## Doing simulation 15
## Doing simulation 16
## Doing simulation 17
## Doing simulation 18
## Doing simulation 19
## Doing simulation 20
## Doing simulation 21
## Doing simulation 22
## Doing simulation 23
## Doing simulation 24
## Doing simulation 25
## Doing simulation 26
## Doing simulation 27
## Doing simulation 28
## Doing simulation 29plot(tpe.s,tpe.K, type="l", xlab="distance", ylab="K(d)", main="K function")
lines(tpe.s, tpe.env$upper, col="red", lty=3)
lines(tpe.s, tpe.env$lower, col="blue", lty=3)
# L(d) Fuction
plot(tpe.s, sqrt(tpe.K/pi)-tpe.s, type="l",ylim=c(-2000,3000), xlab="distance", ylab="L(d)", main="L function")
lines(tpe.s, sqrt(tpe.env$upper/pi)-tpe.s, col="red", lty=3)
lines(tpe.s, sqrt(tpe.env$lower/pi)-tpe.s, col="blue", lty=3)
Bivariate K function (or Cross K function)
The function assesses the spatial relationship between two event points under the null hypothesis of spatial independence and uses the same method as a univariate analysis to test statistical significance. The bivariate K function is defined as the expected number of points for event A within a given distance to arbitrary points of event B divided by the overall density of event A points.
# Plotting points and the boundary
par(mfrow=c(1,2))
polymap(Pts_bnd2, col="light gray")
pointmap(Pts1, xlab="X", ylab="Y", pch=21, bg="yellow", add=T)
polymap(Pts_bnd2, col="light gray")
pointmap(Pts2, xlab="X", ylab="Y", pch=21, bg="red", add=T)
s <- seq(0, 5000 ,50)
tpe.k12 <-k12hat(Pts1, Pts2, Pts_bnd2, s)
plot(s, tpe.k12, type="l",xlab="distance", ylab="K(d)", main="Bivariate K function")
env12<-Kenv.tor(Pts1, Pts2, Pts_bnd2, nsim=49, s)## Doing shift 1 / 49
## Doing shift 2 / 49
## Doing shift 3 / 49
## Doing shift 4 / 49
## Doing shift 5 / 49
## Doing shift 6 / 49
## Doing shift 7 / 49
## Doing shift 8 / 49
## Doing shift 9 / 49
## Doing shift 10 / 49
## Doing shift 11 / 49
## Doing shift 12 / 49
## Doing shift 13 / 49
## Doing shift 14 / 49
## Doing shift 15 / 49
## Doing shift 16 / 49
## Doing shift 17 / 49
## Doing shift 18 / 49
## Doing shift 19 / 49
## Doing shift 20 / 49
## Doing shift 21 / 49
## Doing shift 22 / 49
## Doing shift 23 / 49
## Doing shift 24 / 49
## Doing shift 25 / 49
## Doing shift 26 / 49
## Doing shift 27 / 49
## Doing shift 28 / 49
## Doing shift 29 / 49
## Doing shift 30 / 49
## Doing shift 31 / 49
## Doing shift 32 / 49
## Doing shift 33 / 49
## Doing shift 34 / 49
## Doing shift 35 / 49
## Doing shift 36 / 49
## Doing shift 37 / 49
## Doing shift 38 / 49
## Doing shift 39 / 49
## Doing shift 40 / 49
## Doing shift 41 / 49
## Doing shift 42 / 49
## Doing shift 43 / 49
## Doing shift 44 / 49
## Doing shift 45 / 49
## Doing shift 46 / 49
## Doing shift 47 / 49
## Doing shift 48 / 49
## Doing shift 49 / 49lines(s, env12$upper, col="red", lty=3)
lines(s, env12$lower, col="blue",lty=3)
plot(s, sqrt(tpe.k12/pi)-s, type="l", ylim=c(-4000,4000),xlab="distance", ylab="D(d)", main="D function")
lines(s, sqrt(env12$upper/pi)-s, col="red", lty=3)
lines(s, sqrt(env12$lower/pi)-s, col="blue",lty=3)
Inhomogeneous Poisson process: Estimation of the Intensity (non-parametric)
### Plotting Kernel Density Map
tpe.kernel<-kernel2d(Pts,Pts_bnd2, 1500, 50,50)## Xrange is 295265 316374
## Yrange is 2761740 2789393
## Doing quartic kernel#polymap(Pts_bnd);
polymap(Pts_bnd2)
image(tpe.kernel, add=T)
# Panel the Maps
par(mfrow=c(1,2))
polymap(Pts_bnd2, col="light gray")
pointmap(Pts, xlab="X", ylab="Y", pch=21, bg="yellow", add=T)
polymap(Pts_bnd2)
image(tpe.kernel, add=T)
dev.off()## RStudioGD
## 2# Finding the optimal bandwidth for kernel smoothing
Mse2d <- mse2d(Pts, Pts_bnd2, nsmse=30, range=6000)
plot(Mse2d$h[1:30],Mse2d$mse[1:30], type="l", xlab="bandwidth", ylab="MSE")
Mse2d$h[which.min(Mse2d$mse)]## [1] 1200