第六章 找出點事件的群聚熱區

讀取資料與初步整理

library(sf);library(GISTools);library(dbscan)
setwd("D:/SA_Book/Lab_Data_Ch6")
dengue=st_read(dsn="Dengue_Case.shp")
town=st_read(dsn="KAOH_town.shp")
#轉換成sp格式
dengue.sp=as(dengue,Class="Spatial")
town.sp=as(town,Class="Spatial")
par(mar=c(0,0,0,0))

問題6-1：核密度估計

請呈現高雄市各地區下半年度的登革熱疫情熱區趨勢分布。其中，不同的參數設定會對結果產生何種影響？

# kde.points函數：計算核密度（設定搜尋半徑為2000公尺）
kde=kde.points(dengue.sp,2000,lims=town.sp)
# 繪圖
par(mar=c(0,0,0,0))
plot(st_geometry(town))
image(kde,col=heat.colors(n=10,alpha=0.7,rev=T),add=T)

# 比較不同搜尋半徑的結果
par(mar=c(0,0,1,0),mfrow=c(1,3))
for(width in c(500,1000,2000)){
  kde=kde.points(dengue.sp,width,lims=town.sp)
  plot(st_geometry(town),xlim=c(162,198)*10^3,ylim=c(249,254)*10^4,main=width)
  image(kde,col=heat.colors(n=10,alpha=0.7,rev=T),add=T)
}

問題6-2：應用DBSCAN進行分群

初探登革熱疫情的空間分布之後，如果想要將所有確診病例做分群，評估哪些病例事件屬於同一個疫情熱區，藉此瞭解傳染病可能的傳播鏈，可使用何種分析方法？例如：知道哪些確診病例的分布位置接近，推估相同感染源的可能性。其中，不同的參數設定會對結果產生何種影響？

# dbscan函數：應用DBSCAN進行分群（設定搜尋半徑為500公尺、分群最少點個數為3）
dbscan=dbscan(st_coordinates(dengue),eps=500,minPts=3)
# 繪圖
par(mar=c(0,0,0,0))
hullplot(st_coordinates(dengue),dbscan,asp=1,main ="",axes =F)
plot(st_geometry(town),add=T)

# 比較不同搜尋半徑的結果
par(mar=c(0,0,1,0),mfrow=c(1,3))
for(eps in c(300,500,1000)){
  dbscan=dbscan(st_coordinates(dengue),eps=eps,minPts=3)
  hullplot(st_coordinates(dengue), dbscan,asp=1,main =eps,axes =F,xlim=c(165,195)*10^3,ylim=c(249,253)*10^4)
  plot(st_geometry(town),add=T)
}

問題6-3：應用OPTICS來區分不同密度的群聚

登革熱病例在不同都市環境下，可能出現不同密集程度的群聚熱區。區分此特徵才能更細緻、更貼近真實情況地來描述登革熱的地區流行情況。請使用合適的分析方法，找出高雄市登革熱在空間中的群聚病例。

# optics函數：應用OPTICS來區分不同密度的群聚（先設定搜尋半徑為1000公尺、分群最少點個數為3）
opt=optics(st_coordinates(dengue),eps=1000,minPts=3)
# 畫出可達距離圖(reachability plot)
plot(opt)
abline(h=400,lty=2,col='red')

par(mar=c(0,0,0,0))
#設定搜尋半徑為400公尺，並使用DBSCAN分群
#P.S.也可以透過extractXi函數，利用陡度(steepness)來分群。
optDBSCAN=extractDBSCAN(opt, eps_cl = 400)
hullplot(st_coordinates(dengue), optDBSCAN,asp=1, main = "",axes =F)
plot(st_geometry(town),add=T)

# 依照400公尺分群的可達距離圖
plot(opt$reachdist[opt$order],col=optDBSCAN$cluster[opt$order]+1,ylab='reachability dist.',pch=20)
abline(h=400,lty=2,col='red')