# Week 9 Density-based Methods

> outline:
> Kernel density estimation (KDE)
> Dual KDE
>
## KDE

計算過程：

1. 在研究區建立均勻網格
2. 設定搜尋半徑(bandwidth, h)
3. 選擇核密度函數(Kernel function)
* Uniform distribution
* Normal distribution
* Negative exponential
* Quartic
* Triangular

依特性選擇適合的核密度函數：

### 常態分佈
* peaks & declines rapidly
* No defined radius; continues across entire grid

### Uniform
* represented by cylinder; all points **in the same radius**

### Quartic(spherical) distribution
* Gradual curve; density highest over point; falls to limit of radius
* 例如：住宅竊盜(也可是uniform)，附近住宅風險高，越遠平滑遞減

### Triangular (conical) distribution
* peaks above the point; falls off in a linear manner to edges of radius
* 例如：商業搶劫(也可是負指數型態)，通常都在商業區域，讓附近區域暴露在相同的風險之下

### Negative exponential distribution
* curve tha falls off rapidly from the peak to a specified radius ---> tightly focused
* 例如：家庭暴力案件的發生點，通常家暴的影響範圍很小，家庭之間不太會互相影響

**不同的方法會產生不同的結果**

### 選擇合適的搜尋半徑

Fixed vs. Adaptive intervals

依照選擇的核密度函數去做變化

若選擇Normal kernel function，則以他的標準差為搜尋半徑

其他方式則用研究區域大小為準

* Fixed intervals:
1. Mean Integrated Squared Error (MISE)


2. Algorithms

3. Range of influence (K-order NNA)

* Adaptive intervals

Conclusion:
* 針對不同的資料型態選擇合適的方法
* 運用各種統計分析方法亦需瞭解其各種方法的基本假設、分析模式的侷限性以及分析結果的正確解讀。
* 透過各種資料型態與不同統計方法的分析，在各種方法間之相互比較，若都有共同的分析結果，較能合理的反映出該地區確實具有地理群聚現象。

## 實習：Dual KDE 實作
分別使用 **splancs** 與 **GISTools** 實作

splancs:

rm(list = ls())
library(rgdal)
library(splancs)
setwd("C:/Users/user/Desktop/SpatialAnalysis")

# spatial points
data1 = read.table("./data/point1.csv", header = T, sep = ",")
data2 = read.table("./data/point2.csv", header = T, sep = ",")
pts1 = as.points(data1[,1],data1[,2])
pts2 = as.points(data2[,1],data2[,2])

# boundary file
data_bnd = read.table("./data/tpe_sqr_bnd.csv", header = T, sep = ",")
Pts_bnd = as.points(data_bnd[,2],data_bnd[,3])

tpekde1 = kernel2d(pts1, Pts_bnd, 2000,200,200)

tpekde2 = kernel2d(pts2, Pts_bnd, 2000,200,200)

diff = (tpekde1$z - tpekde2$z)

tpekde_diff = list(x=tpekde1$x, y=tpekde1$y, z=diff, h0=tpekde1$h0, kernel=tpekde1$kernel)
polymap(Pts_bnd);image(tpekde_diff, add=T)


GISTools:

# shapefile
TaiwanCounty = readOGR(dsn="./data", layer = "Taiwan_county")
index = TaiwanCounty$COUNTY=="台北市"
TaipeiTown = TaiwanCounty[index,]

proj = CRS("+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +units=m +no_defs")

pts_gis1 = SpatialPoints(pts1, proj4string = proj)
pts_gis2 = SpatialPoints(pts2, proj4string = proj)

pts1.kernel = kde.points(pts_gis1, 2000,200, lims = TaipeiTown)
pts2.kernel = kde.points(pts_gis2, 2000,200, lims = TaipeiTown)

pts1KDE.raster = raster(pts1.kernel)
pts2KDE.raster = raster(pts2.kernel)

KDE.Diff = pts1KDE.raster - pts2KDE.raster

plot(KDE.Diff, axes=FALSE)
masker.gis = poly.outer(KDE.Diff, TaipeiTown)
add.masking(masker.gis, color = "aliceblue");plot(TaipeiTown, add=T)