點模式分析(二)

點模式分析（一）的基礎上，繼續點模式分析。

包含ANN、G函數、K函數的R語言實現。首先需要下載rspatial包，自帶數據。

#devtools::install_github('rspatial/rspatial')

library(rspatial)

city <- sp_data('city')

crime <- sp_data('crime')

plot(city, col='light blue')

points(crime, col='red', cex=.5, pch='+')

#獲取坐標

xy<-coordinates(crime)

dim(xy)

#移除重複坐標

xy <- unique(xy)

dim(xy)

#求x,y坐標平均值

mc <- apply(xy, 2, mean)

# 標準距離

sd <- sqrt(sum((xy[,1] - mc[1])^2 +(xy[,2] - mc[2])^2) / nrow(xy))

plot(city, col='light blue')

points(crime, cex=.5)

#畫出xy坐標平均點

points(cbind(mc[1], mc[2]), pch='*',col='red', cex=5)

# 畫圓

bearing <- 1:360 * pi/180

cx <- mc[1] + sd * cos(bearing)

cy <- mc[2] + sd * sin(bearing)

circle <- cbind(cx, cy)

lines(circle, col='red', lwd=2)

 

 

#----柵格化研究區域並計算crime密度

#研究區面積

CityArea <- raster::area(city)

#研究勻內點密度

dens <- nrow(xy) / CityArea

#創建一個RasterLayer對象,extent與city相同

#此時cell沒有值

r <- raster(city)

#設置柵格解析度

res(r) <- 1000

#按1000米大小網格化city範圍

r <- rasterize(city, r)

plot(r)

quads <- as(r, 'SpatialPolygons')

#柵格值全為1，沒有計算crime點數

plot(quads, add=TRUE)

points(crime, col='red', cex=.5)

#設置fun=count，計算每個cell的crime個數

nc <- rasterize(coordinates(crime), r,fun='count', background=0)

plot(nc)

plot(city, add=TRUE)

#利用r裁剪nc

ncrimes <- mask(nc, r)

plot(ncrimes)

plot(city, add=TRUE)

 

#ANN----

#計算任意兩點之間的距離

d <- dist(xy)

dm <- as.matrix(d)

dm[1:5, 1:5]

#設置對角線數值為NA,為什麼

diag(dm) <- NA

#計算每一點到其它點的最近距離

dmin <- apply(dm, 1, min, na.rm=TRUE)

head(dmin)

ANN <- mean(dmin)

#--基於密度的計算-

#補充：A ppp objecthas the coordinates of the points and the analysis 『window』 (study region). #Toassign the points locations we need to extract the coordinates from our SpatialPoints#object. To set the window, we first need to to coerce our SpatialPolygons intoan 『owin』 #object. We need a function from the maptools package for thiscoercion.

#-基於密度計算,

library(spatstat)

library(maptools)

cityOwin <- as.owin(city)

pts <- coordinates(crime)

p <- ppp(pts[,1], pts[,2],window=cityOwin)

plot(p)

ds <- density(p)

plot(ds, main='crime density')

#人口密度

census <-sp_data("census2000.rds")

census$area <- area(census)

#persons per feet2 to persons per mile2

census$area <- census$area/27878400

census$dens <- census$POP2000 /census$area

p <- coordinates(census)

#所有多邊形合併為一個

win <- aggregate(census)

plot(census)

points(p, col='red', pch=20, cex=.25)

plot(win, add=TRUE, border='blue', lwd=3)

owin <- as.owin(win)

#Warning message:1 point was rejected

pp <- ppp(p[,1], p[,2], window=owin,marks=census$dens)

#去除這個點

sp <- SpatialPoints(p,proj4string=CRS(proj4string(win)))

library(rgeos)

i <- gIntersects(sp, win, byid=TRUE)

which(!i)

## [1] 588

plot(census)

points(sp)

#顯示該點

points(sp[!i,], col='red', cex=3, pch=20)

#把該點去除

pp <- ppp(p[i,1], p[i,2], window=owin,marks=census$dens[i])

plot(pp)

plot(city, add=TRUE)

s <- Smooth.ppp(pp)

plot(s)

plot(city, add=TRUE)

 

#---犯罪類型可視化----

par(mfrow=c(2,2), mai=c(0.25, 0.25, 0.25,0.25))

for (offense in c("Auto Theft","Drunk in Public", "DUI", "Arson")) {

 plot(city, col='grey')

 acrime <- crime[crime$CATEGORY == offense, ]

 points(acrime, col = "red")

 title(offense)

}

crime$fcat <- as.factor(crime$CATEGORY)

w <- as.owin(city)

xy <- coordinates(crime)

mpp <- ppp(xy[,1], xy[,2], window = w,marks=crime$fcat)

spp <- split(mpp)

plot(spp[1:4], main='')

plot(density(spp[1:4]), main='')

#--K函數及檢驗

spatstat.options(checksegments = FALSE)

ktheft <- Kest(spp$"AutoTheft")

ketheft <- envelope(spp$"AutoTheft", Kest)

ktheft <- Kest(spp$"Arson")

ketheft <-envelope(spp$"Arson", Kest)

par(mfrow=c(1,2))

plot(ktheft)

plot(ketheft)

KS.arson <- cdf.test(spp$Arson, ds)

KS.arson

KS.drunk <- cdf.test(spp$'Drunk inPublic', ds)

KS.drunk

參考資料：https://rspatial.org/raster/analysis/8-pointpat.html#