空間分析 作業六
地理四 B03208011 杜承軒
＊初步環境建置與讀取檔案
library(sp);library(rgdal);library (splancs);library(GISTools);library(ggplot2) setwd("C:/Users/user/Desktop/SA") source("ST_functions.R")
問題一、建立自訂函數
Q1-1.建立距離與時間矩陣Dist_Matrix=function(X,Y){ n=length(X);M=matrix(0,n,n) for(i in seq(1,n-1)){ for(j in seq(i+1,n)){ M[i,j]=sqrt((X[i]-X[j])^2+(Y[i]-Y[j])^2) }} return(M) } Time_Matrix=function(t){ n=length(t);M=matrix(0,n,n) for(i in seq(1,n-1)){ for(j in seq(i+1,n)){ M[i,j]=abs(t[i]-t[j]) }} return(M) }Q1-2.建立時空關連圖
sptplot=function(X,Y,t){ n=length(t); dismat=Dist_Matrix(X,Y); timmat=Time_Matrix(t) plot(0,0,"n",xlim=c(0,max(timmat)),ylim=c(0,max(dismat)), xlab="Time-Distance",ylab="Space-Distance",main="Space-Time Association") for(i in seq(1,n-1)){ for(j in seq(i+1,n)){ points(timmat[i,j],dismat[i,j],col="red",pch=20) }} abline(h=sum(dismat)*2/(n-1)^2,col="blue",lwd=2) #以平均距離來劃分 abline(v=sum(timmat)*2/(n-1)^2,col="blue",lwd=2) #以平均時間間隔來劃分 }Q1-3.建立 Jacquez‘s k-NN曲線，含隨機區間
Jacquez.plot=function(X,Y,t){ res.T=data.frame();mon.T=data.frame() for(k in 1:10){ out=Jacquez.test(X,Y,t,k,99) res.T=rbind(res.T,data.frame(k=factor(k),t=out$Jacquez.T)) mon.T=rbind(mon.T,data.frame(k=factor(k),t=out$Freq)) } ggplot()+geom_boxplot(data=mon.T,aes(x=k, y=t,col=k),fill="#AAAAAA",size=.8)+ geom_point(data=res.T,aes(x=k, y=t),col="red",size=3)+ geom_line(data=res.T,aes(x=c(1:10), y=t),col="red")+ labs(title="Jacquez's k-Nearest Neighbor Plot",x="k-nearest neighbor",y="Jacquez T")+ guides(col=FALSE)+ theme_light()+theme(plot.title = element_text(size=16,hjust = 0.5)) }
問題二、利用台北市某疾病的疫情資料進行時空風險評估
#匯入資料 pts=readOGR(dsn = "ST_Data", layer = "point_tpe", encoding="big5",verbose = F) pdata=pts@data pd.x=pdata$X pd.y=pdata$Y pd.t=pdata$WEEKQ2-1.使用Q1的函數進行時空分析 2-1-1. 距離與時間矩陣
pd.dis=Dist_Matrix(pd.x,pd.y) pd.tim=Time_Matrix(pd.t) pd.tim[1:10,1:10]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] ## [1,] 0 13 6 6 8 5 3 2 12 0 ## [2,] 0 0 19 19 5 18 10 15 1 13 ## [3,] 0 0 0 0 14 1 9 4 18 6 ## [4,] 0 0 0 0 14 1 9 4 18 6 ## [5,] 0 0 0 0 0 13 5 10 4 8 ## [6,] 0 0 0 0 0 0 8 3 17 5 ## [7,] 0 0 0 0 0 0 0 5 9 3 ## [8,] 0 0 0 0 0 0 0 0 14 2 ## [9,] 0 0 0 0 0 0 0 0 0 12 ## [10,] 0 0 0 0 0 0 0 0 0 0pd.dis[1:7,1:7]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] ## [1,] 0 4062.697 9329.303 9329.303 2914.570 9613.9633 9489.2173 ## [2,] 0 0.000 5339.513 5339.513 3584.772 5610.8525 5441.1224 ## [3,] 0 0.000 0.000 0.000 8635.163 307.8246 758.1669 ## [4,] 0 0.000 0.000 0.000 8635.163 307.8246 758.1669 ## [5,] 0 0.000 0.000 0.000 0.000 8870.2181 8559.5073 ## [6,] 0 0.000 0.000 0.000 0.000 0.0000 646.4093 ## [7,] 0 0.000 0.000 0.000 0.000 0.0000 0.00002-1-2. 時空關連圖
sptplot(pd.x,pd.y,pd.t)
藍線為平均距離與時間間隔。從關聯圖觀察時空資料分散，看似沒有時空相關的趨勢。
2-1-3. Jacquez’s k-NN曲線
Jacquez.plot(pd.x,pd.y,pd.t)
紅點為計算出來的觀測值，箱型圖為模擬99次的結果。除了k=1時不明顯，其餘紅點都落在盒狀圖中，以Jacquez’s方法來說並無顯著時空相關。
Q2-2.其他時空分析
2-2-1. Knox Test
KT=KnoxM.test(pd.x,pd.y,pd.t,2000,10,99) paste("Knox.T 統計值為：",KT$Knox.T,"； 模擬後的p-value為：",KT$Simulated.p.value)## [1] "Knox.T 統計值為： 65 ； 模擬後的p-value為： 0.86"hist(KT$Freq, freq=T,main="Monte Carlo of Knox.T",xlab="Knox.T") abline(v=KT$Knox.T, col="red")
不顯著，在Knox Test中方法來說並無顯著時空群聚。
2-2-2. Mantel Test
MT=Mantel.test(pd.x,pd.y,pd.t,0.1,0.1,99) paste("Mantel.T 統計值為：",MT$Mantel.T,"； 模擬後的p-value為：",MT$Simulated.p.value)## [1] "Mantel.T 統計值為： 102.060167431139 ； 模擬後的p-value為： 0.01"hist(MT$Freq, freq=T,main="Monte Carlo of Mantel.T",xlab="Mantel.T") abline(v=MT$Mantel.T, col="red")
顯著，在Mantel Test中方法來說，有時空群聚的現象。
2-2-3. Jacquez Test (k=1時)
JT=Jacquez.test(pd.x,pd.y,pd.t,1,99) paste("Jacquez.T 統計值為：",JT$Jacquez.T,"； 模擬後的p-value為：",JT$Simulated.p.value)## [1] "Jacquez.T 統計值為： 1.5 ； 模擬後的p-value為： 0.09"hist(JT$Freq, freq=T,main="Monte Carlo of Jacquez.T",xlab="Jacquez.T") abline(v=JT$Jacquez.T, col="red")
不顯著，在Jacquez Test中方法來說並無顯著時空群聚。
2-2-4. Space-time K function
#邊界資料 bondary=read.table("Data/tpe_sqr_bnd.csv", header=TRUE, sep=",") BDP=as.points(bondary[,2], bondary[,3]) pd.pt=as.points(pd.x,pd.y) #畫出D圖與D0圖 Kst=stkhat(pd.pt, pd.t, BDP, c(0,50),seq(0,2000,50),seq(0,20,1)) Ks=matrix(Kst$ks) Kt=matrix(Kst$kt) KsKt=Ks%*% t(Kt) Dst=Kst$kst - KsKt persp(Kst$s,Kst$t,Dst,theta=-30, phi =15, expand = .6, xlab="spatial distance", ylab="temporal distance", zlab="D", ticktype ="detailed" ,main="D")
D0st=Dst/KsKt persp(Kst$s,Kst$t,D0st,theta=-30,phi = 15, expand = .6, xlab="spatial distance", ylab="temporal distance", zlab="D0", ticktype ="detailed",main="D0" )
#畫出標準化後的殘差 se=stsecal(pd.pt, pd.t, BDP, c(0,50),seq(0,2000,50),seq(0,20,1)) res=Dst/ se plot(KsKt,res, xlab="K(s)K(t)", ylab=" standardized residuals R(s,t)",main="standardized residual") abline(h=c(-2,2), lty=2,col="red")
觀察D0的圖，能看出在距離在500公尺以內，時間間隔5單位以下，時空群聚的現象越明顯。隨時空距離上升，群聚就不明顯了。配合殘差的圖，可以看出只有在兩者乘積很小時才有顯著的現象。
