# Week 8 Distanced-based Methods

> outline:
> F(d) function
> Ripley's K function
>

## F(d) function

概念：
隨機點到觀察值距離小於d的個數/觀察點數
也是百分比的概念

如果觀察點本身是clustered，
則隨機點到觀察點的距離會是長距離的比較多

計算步驟：
1. 產生隨機點
2. 計算隨機點到觀察點的最小距離
3. 計算F(d)

### F(d) 的解讀

1. 若觀察點的分布是群聚的，F(d)的圖形一開始會緩慢上升，在長距離時才快速上升
2. 若觀察點的分布是分散的，F(d)的圖形一開始會快速上升

### F(d) 的顯著性檢定

跟G(d)的檢定方式一樣

R 實作 F(d) 與套件比較：

# Lab: F function 實作

library(sp)
library(rgdal)
library(spatstat)
library(maptools)
library(GISTools)

# 1. read file and convert to ppp
School <- readOGR(dsn =  "./data", layer = "Schools", encoding="utf8")
TWCounty <- readOGR(dsn = "./data", layer = "Taiwan_county", encoding="utf8")

x.coor = School$coords.x1
y.coor = School$coords.x2

index = TWCounty$COUNTY=="台南市"
TN_BND = TWCounty[index,]
TN = gBuffer(TN_BND, width = 500)
TN_BND2 = as.owin.SpatialPolygons(TN)

school.ppp = ppp(x.coor, y.coor, window = TN_BND2)

產生隨機點：

rdppp = rpoint(200, win =  TN_BND2)

計算兩點之間的最短距離

nnd = nncross(rdppp, school.ppp)

計算F(d)

F_1 = ecdf(nnd[,1])

蒙地卡羅顯著性檢定：

plot(F_1, main="F function", col="red", cex=0, verticals=T, xaxs="i", yaxs="i", xlab="distance")

for(i in 1:99){
  rp = rpoint(424, win = TN_BND2)
  nnd.s = nncross(rdppp, rp)
  F_RND = ecdf(nnd.s[,1])
  lines(F_RND, col='gray', cex=0, verticals=T, xaxs="i", yaxs="i")
}


套件比較：

plot(envelope(school.ppp, fun = Fest, nsim = 90, nrank = 10))


## Ripley's K function

G F 都只考慮 nearest neighbor，淡化了點之間的實際分布，為了能看見不同鄰近程度之間的效果，而設計了 K(d) function

計算步驟：

1. 以每個點中心畫半徑為d的圓
2. 計算每個環域包含多少除了自己(中心點)以外的點
3. 計算每個環域包含點數的平均
4. 除上點密度 (點個數除以面積) (不討論密度的不均勻性)

### K(d)圖形的解讀

1. 當環域半徑越大，K(d)值也會越大
2. 若有群聚，會出現d上升而K值平緩的情況，一直達到clustered size(群聚的最大環域)

### 邊緣校正

邊緣校正發生在點出現在邊界附近時，因為畫環域範圍時會超過指定邊界，邊界外也可能有他的鄰近點，因此需要進行邊緣校正


常見的邊緣校正方法：
1. Simulate outer boundary values
將每個需要邊緣校正的點依邊界鏡射到邊界外面，計算需要邊界校正的點時計入那些鏡射點

2. Reduce analysis area

### K(d)的顯著性檢定


在完全隨機分布(CSR)下，L(d)=0

推導K(d)的期望值

以 K(d)=0 為判斷標準，以上為clustered，以下為 dispersed

Q: 為什麼 L(d) 會出現負值？
ANS: 到最後 K(d) 會抵達上限


R 實作 K function 與套件比較

# HK: K function 實作
setwd("C:/Users/user/Desktop/SpatialAnalysis")
# 1. read file and convert to ppp
School <- readOGR(dsn =  "./data", layer = "Schools", encoding="utf8")

bnd = bbox(School)
x.coor = School$coords.x1
y.coor = School$coords.x2
x.range = bnd[1,] # units: m
y.range = bnd[2,]


school.ppp = ppp(x.coor, y.coor, xrange=x.range, yrange=y.range)

# 2. generate buffer with radius d
d = 100
school.buffer = gBuffer(School, width = d, byid = T)

# 3. count points in each buffer without self
points = poly.counts(School, school.buffer)-1

# 4. 每個buffer內的點個數加起來，再除以環域數
mean = sum(points)/length(points)

# 5. 再除以點密度 (點個數/面積)
K = mean / ((x.range[2]-x.range[1])*(y.range[2]-y.range[1]))

L = sqrt(K/pi)-d

# 重複計算2~5
d=0
xside=x.range[2]-x.range[1]
yside=y.range[2]-y.range[1]
area=xside*yside
density = 424/area

while(d<12000) {
  school.buffer = gBuffer(School, width = d, byid = T)
  points = poly.counts(School, school.buffer)-1
  mean = sum(points)/424
  K = mean / density
  if(d == 0){
    L[1] = sqrt(K/pi)-d
  }else{
    L[(d/100)+1] = sqrt(K/pi)-d
  }
  d=d+100
}

plot(L, type='l', xlab='Distance (百公尺)', ylim=c(-500,4000))

# 模擬99次
for(i in 1:99){
  d=0
  L.sim=c()
  while(d<12000){
    rp = as.SpatialPoints.ppp(rpoint(424, win = school.ppp))
    rp.buffer = gBuffer(rp, width=d, byid=T)
    points = poly.counts(rp, rp.buffer)-1
    mean = sum(points)/424
    K=mean/density
    if(d==0){
      L.sim[1]=sqrt(K/pi)-d
    }else{
      L.sim[(d/100)+1]=sqrt(K/pi)-d
    }
    d=d+100
  }
  lines(L.sim, col='gray')
}


套件實作：

plot(envelope(school.ppp, fun=Lest), .-r~r)