空間分析　 作業五

地理四 B03208011 杜承軒

＊初步環境建置與讀取檔案

library(GISTools)
library(rgdal)
library(splancs)
library(raster)
library(dbscan)
setwd("C:/Users/user/Desktop/SA")
#讀取資料與格式轉換
tpe=readOGR(dsn = "Data", layer = "TPE_TOWN", encoding="utf8",verbose = F)
#速食店
FF=readOGR(dsn = "Homework", layer = "Tpe_Fastfood", encoding="utf8",verbose = F)
MIC=subset(FF,FF@data$STORE=="MIC")
KFC=subset(FF,FF@data$STORE=="KFC")
FF.pt=as.points(FF@coords[,1],FF@coords[,2]) 
MIC.pt=as.points(MIC@coords[,1],MIC@coords[,2]) 
KFC.pt=as.points(KFC@coords[,1],KFC@coords[,2]) 
#邊界資料
bondary <- read.table("Data/tpe_sqr_bnd.csv", header=TRUE, sep=",")
BDP <- as.points(bondary[,2], bondary[,3])
#畫圖
plot(tpe, col="#DDDDDD")
pointmap(MIC.pt, pch=21, bg="yellow", add=T)
pointmap(KFC.pt, pch=21, bg="red", add=T)
legend(313000,2787000,c("麥當勞","肯德基"),pch=20,col=c("yellow","red"),box.col="white")

問題一、
參考 Reading_Dual.KDE.pdf 這篇論文關於market dominance的定義，用dual KDE分析台北市MIC 或 KFC 市場主導的空間分布。

算出MIC-KFC的dual KDE
＊(1)找出合適的寬度

Mse2d <- mse2d(FF.pt,BDP, nsmse=30, range=6000)
plot(Mse2d$h[1:30],Mse2d$mse[1:30], type="l", 
    xlab="bandwidth", ylab="MSE",main="Find Optimal Bandwidth")
Mse2d <- mse2d(MIC.pt,BDP, nsmse=30, range=6000)
lines(Mse2d$h[1:30],Mse2d$mse[1:30], type="l",col="red",lty=3)
Mse2d <- mse2d(KFC.pt,BDP, nsmse=30, range=6000)
lines(Mse2d$h[1:30],Mse2d$mse[1:30], type="l",col="blue",lty=3)
legend(4000,20,c("所有速食店","MIC","KFC"),col = c(1,2,4),pch=20)
abline(v=1000,col="gold",lty=2)

＊(2)做Dual KDE

kde.mic =kde.points(MIC,1000,n=2000, lims=tpe) 
kde.kfc =kde.points(KFC,1000,n=2000, lims=tpe) 
kde.mic_r=raster(kde.mic)  
kde.kfc_r=raster(kde.kfc)
kde.diff=kde.mic_r-kde.kfc_r  
plot(kde.diff, axes=FALSE,col=rev(heat.colors(15)),bty="n")  
masker=poly.outer(kde.diff,tpe)
add.masking(masker, col="white")
plot(tpe,add=TRUE) 

紅色表示麥當勞市場主導，橘色為兩者實力相當，偏黃色表示肯德基市場主導。

---

問題二、
利用DBSCAN找出 MIC 與 KFC 的空間群聚。並討論不同參數設定 (eps, minPts)，對於群聚結果的敏感性。

＊寫成函數

xeps=function(pt,u,v0){
 par(mfrow=c(2,2));par(mar=c(0,0,0,0))
 for (i in u){
   res=dbscan(pt, eps = i, minPts = v0)
   hullplot(pt, res, asp=1,add=T,xlim=c(bbox(tpe)[1,]), ylim=c(bbox(tpe)[2,]),main ="",axes =F)
   plot(tpe,col=rgb(0,0,0,0), add=T)
   pointmap(pt, col = res$cluster +1 , pch = res$cluster +1 , add=T)
   legend(310000,2787000,paste("搜尋半徑=",i,"\n臨界點數=",v0),box.col=rgb(0,0,0,0))
 }
 par(mfrow=c(1,1));par(mar=c(5,4,4,2.5))
}

xminPts=function(pt,u0,v){
 par(mfrow=c(2,2));par(mar=c(0,0,0,0))
 for (i in v){
   res=dbscan(pt, eps = u0, minPts = i)
   hullplot(pt, res, asp=1,add=T,xlim=c(bbox(tpe)[1,]), ylim=c(bbox(tpe)[2,]),main ="",axes =F)
   plot(tpe,col=rgb(0,0,0,0), add=T)
   pointmap(pt, col = res$cluster +1 , pch = res$cluster +1 , add=T)
   legend(310000,2787000,paste("搜尋半徑=",u0,"\n臨界點數=",i),box.col=rgb(0,0,0,0))
 }
 par(mfrow=c(1,1));par(mar=c(5,4,4,2.5))
}

＊麥當勞

kNNdistplot(MIC.pt, k = 3);mtext("麥當勞\n")
abline(h=c(1000,1500,2000,3000), col = "red", lty=2)

xeps(MIC.pt,c(1000,1500,2000,3000),3)

xminPts(MIC.pt,1500,c(2,3,5,7))

＊肯德基

kNNdistplot(KFC.pt, k = 3);mtext("肯德基\n")
abline(h=c(1000,2000,3000,5000), col = "red", lty=2)

xeps(KFC.pt,c(1000,2000,3000,5000),3)

xminPts(KFC.pt,2000,c(2,3,5,7))

＊速食店

kNNdistplot(FF.pt, k = 3);mtext("速食店\n")
abline(h=c(1000,1500,2000,3000), col = "red", lty=2)

xeps(FF.pt,c(1000,1500,2000,3000),3)

xminPts(FF.pt,1500,c(2,3,5,7))

＊固定臨界點數，逐漸擴大搜尋半徑，會使群聚的群數有變少的趨勢，然而每一群的範圍也有擴大的走向。
＊固定搜尋半徑，逐漸增加臨界點數，會讓群聚的群數有變少，越不易偵測出群聚的現象。