Final report
jie-en
2019年5月31日
檔案處理
library(rgdal)
library(GISTools)
library(spdep)
accidentA2<-read.csv("C:/Users/ufmpo/Downloads/歷史交通事故A2.csv",header = T,stringsAsFactors = FALSE)
accidentA1<-read.csv("C:/Users/ufmpo/Downloads/歷史交通事故A1.csv",header = T,stringsAsFactors = FALSE)
accidentAAG = read.csv("C:/Users/ufmpo/Desktop/鬼月.csv",header = T,stringsAsFactors = FALSE)
accidentNIT = read.csv("C:/Users/ufmpo/Desktop/accident.csv",header = T,stringsAsFactors = FALSE)
RIP<-read.csv("C:/Users/ufmpo/Downloads/殯葬設施_北北基.csv",header=T,stringsAsFactors = FALSE)
RIP<-RIP[,1:10]
accident_TWN_A1<-accidentA1[(substring(accidentA1$發生地點,1,3)=="臺北市"|substring(accidentA1$發生地點,1,3)=="新北市"|substring(accidentA1$發生地點,1,3)=="基隆市"),]
accident_TWN_A1$受傷人數<-substring(accident_TWN_A1$死亡受傷人數,7,8)
accident_TWN_A1$月份<-substring(accident_TWN_A1$發生時間,5,6)
accident_TWN_A2<-accidentA2[(substring(accidentA2$發生地點,1,3)=="臺北市"|substring(accidentA2$發生地點,1,3)=="新北市"|substring(accidentA2$發生地點,1,3)=="基隆市"),]
accident_TWN_A2$受傷人數<-substring(accident_TWN_A2$死亡受傷人數,7,8)
accident_TWN_A2$死亡人數<-substring(accident_TWN_A2$死亡受傷人數,3,3)
accident_TWN_A2$月份<-substring(accident_TWN_A2$發生時間,5,6)
accident_TWN_AAG<-accidentAAG[(substring(accidentAAG$發生地點,1,3)=="臺北市"|substring(accidentAAG$發生地點,1,3)=="新北市"|substring(accidentAAG$發生地點,1,3)=="基隆市"),]
accident_TWN_AAG$受傷人數<-substring(accident_TWN_AAG$死亡受傷人數,7,8)
accident_TWN_AAG$死亡人數<-substring(accident_TWN_AAG$死亡受傷人數,3,3)
accident_TWN_AAG$月份<-substring(accident_TWN_AAG$發生時間,5,6)
accident_TWN_NIT<-accidentNIT[(substring(accidentAAG$發生地點,1,3)=="臺北市"|substring(accidentAAG$發生地點,1,3)=="新北市"|substring(accidentAAG$發生地點,1,3)=="基隆市"),]
accident_TWN_NIT$受傷人數<-substring(accident_TWN_NIT$死亡受傷人數,7,8)
accident_TWN_NIT$死亡人數<-substring(accident_TWN_NIT$死亡受傷人數,3,3)
accident_TWN_NIT$月份<-substring(accident_TWN_NIT$發生時間,5,6)
accident_TWN_NIT$時間 = substring(accident_TWN_NIT$發生時間,12,13)
accident_TWN_NIT = subset(accident_TWN_NIT, accident_TWN_NIT$時間 %in% c("21","22","23","00","01","02","03"))
.<-xtabs(~accident_TWN_A2$月份)
#barplot(.)
.<-xtabs(~accident_TWN_A1$月份)
#barplot(.)
TW <- readOGR(dsn = "C:/Users/ufmpo/Desktop/空間分析/2.Geospatial visualization/Data2/Popn_TWN2.shp", layer = "Popn_TWN2", encoding="utf8",verbose = F)
TWN<- subset(TW, TW$COUNTY %in% c("臺北市","新北市","基隆市"))
A1<-SpatialPointsDataFrame(cbind(accident_TWN_A1$經度,accident_TWN_A1$緯度),data = accident_TWN_A1,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
A1<-spTransform(A1,proj4string(TWN))
A2<-SpatialPointsDataFrame(cbind(accident_TWN_A2$經度,accident_TWN_A2$緯度),data = accident_TWN_A2,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
A2<-spTransform(A2,proj4string(TWN))
AAG<-SpatialPointsDataFrame(cbind(accident_TWN_AAG$經度,accident_TWN_AAG$緯度),data = accident_TWN_AAG,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
AAG<-spTransform(AAG,proj4string(TWN))
ANIT<-SpatialPointsDataFrame(cbind(accident_TWN_NIT$經度,accident_TWN_NIT$緯度),data = accident_TWN_NIT,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
ANIT<-spTransform(ANIT,proj4string(TWN))
RIP<-SpatialPointsDataFrame(cbind(RIP$經度,RIP$緯度),data = RIP,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
RIP<-spTransform(RIP,proj4string(TWN))
#plot(TWN, border="darkgrey")
#plot(RIP,col="green",pch=20,add=T)
#plot(A2,col="blue",pch=20,add=T,cex=0.35)
#plot(A1,col="red",pch=20,add=T,cex=0.35)Local Indicator of Spatial Association(LISA)
bf = gBuffer(RIP,byid = T ,width = 1000)
bf$A1 = poly.counts(A1,bf)
bf$A2 = poly.counts(A2,bf)
bf$AAG = poly.counts(AAG,bf)
bf$ANIT = poly.counts(ANIT,bf)
for (i in 1:nrow(bf)) bf$id[i] = as.character(TWN$TOWN_ID[bf$縣市[i] == TWN$COUNTY & bf$鄉鎮市區[i] == TWN$TOWN] )
A11 = data.frame(xtabs(bf$A1~bf$id),stringsAsFactors = F)
A22 = data.frame(xtabs(bf$A2~bf$id),stringsAsFactors = F)
AAG = data.frame(xtabs(bf$AAG~bf$id),stringAsFactors = F)
ANIT = data.frame(xtabs(bf$ANIT~bf$id),stringsAsFactors = F)
A11$bf.id = as.character(A11$bf.id)
A22$bf.id = as.character(A22$bf.id)
AAG$bf.id = as.character(AAG$bf.id)
ANIT$bf.id = as.character(ANIT$bf.id)
for (i in 1:nrow(TWN)){
if (any(TWN$TOWN_ID[i] == A11$bf.id)==TRUE){
TWN$A1[i] = A11$Freq[TWN$TOWN_ID[i] == A11$bf.id]
}else{
TWN$A1[i] = 0
}
if (any(TWN$TOWN_ID[i] == A22$bf.id)==TRUE){
TWN$A2[i] = A22$Freq[TWN$TOWN_ID[i] == A22$bf.id]
}else{
TWN$A2[i] = 0
}
if (any(TWN$TOWN_ID[i] == AAG$bf.id)==TRUE){
TWN$AAG[i] = AAG$Freq[TWN$TOWN_ID[i] == AAG$bf.id]
}else{
TWN$AAG[i] = 0
}
if (any(TWN$TOWN_ID[i] == ANIT$bf.id)==TRUE){
TWN$ANIT[i] = ANIT$Freq[TWN$TOWN_ID[i] == ANIT$bf.id]
}else{
TWN$ANIT[i] = 0
}
}
TWN$A3 = TWN$A1 + TWN$A2
bf$A3 = bf$A1 + bf$A2########################
#LISA(unit: per buffer)#
########################
coords = coordinates(bf)
bf.nb = knn2nb(knearneigh(coords, k=2))
bf.nb.w = nb2listw(bf.nb,zero.policy=T)
LISA1 = localmoran(bf$A1, bf.nb.w,zero.policy = T,alternative = "two.sided")
LISA2 = localmoran(bf$A2, bf.nb.w,zero.policy = T,alternative = "two.sided")
LISA6 = localmoran(bf$AAG, bf.nb.w,zero.policy = T,alternative = "two.sided")
LISA7 = localmoran(bf$ANIT, bf.nb.w,zero.policy = T,alternative = "two.sided")
LISA8 = localmoran(bf$A3, bf.nb.w,zero.policy = T,alternative = "two.sided")
#LISA of A1
diff1 = bf$A1 - mean(bf$A1) # diff看自己和平均比起來算是H還是L
z1 = LISA1[,4]
quad1 = c()
quad1[diff1>0 & z1>0] = 1 # H-H
quad1[diff1<0 & z1>0] = 2 # L-L
quad1[diff1>0 & z1<0] = 3 # H-L
quad1[diff1<0 & z1<0] = 4 # L-H
quad1[LISA1[, 5]>0.05]=5 # 不顯著，設定雙尾所以用0.05比較就可以
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(TWN, main = "LISA Map",xlim = c(271884.7,359423.0),ylim = c(2727079,2805226),border = "grey")
plot(bf, border="grey", col=colors[quad1],add = T)
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
#LISA of A2
diff2 = bf$A2 - mean(bf$A2) # diff看自己和平均比起來算是H還是L
z2 = LISA2[,4]
quad2 = c()
quad2[diff2>0 & z2>0] = 1 # H-H
quad2[diff2<0 & z2>0] = 2 # L-L
quad2[diff2>0 & z2<0] = 3 # H-L
quad2[diff2<0 & z2<0] = 4 # L-H
quad2[LISA2[, 5]>0.05]=5 # 不顯著，設定雙尾所以用0.05比較就可以
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(TWN, main = "LISA Map",xlim = c(271884.7,359423.0),ylim = c(2727079,2805226),border = "grey")
plot(bf, border="grey", col=colors[quad2], main = "LISA Map",add = T)
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
#LISA of AAG
diff6 = bf$A6 - mean(bf$A6) # diff看自己和平均比起來算是H還是L
z6 = LISA6[,4]
quad6 = c()
quad6[diff6>0 & z6>0] = 1 # H-H
quad6[diff6<0 & z6>0] = 2 # L-L
quad6[diff6>0 & z6<0] = 3 # H-L
quad6[diff6<0 & z6<0] = 4 # L-H
quad6[LISA6[, 5]>0.05]=5 # 不顯著，設定雙尾所以用0.05比較就可以
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(TWN, main = "LISA Map",xlim = c(271884.7,359423.0),ylim = c(2727079,2805226),border = "grey")
plot(bf, border="grey", col=colors[quad6], main = "LISA Map",add = T)
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
#LISA of ANIT
diff7 = bf$A7 - mean(bf$A7) # diff看自己和平均比起來算是H還是L
z7 = LISA7[,4]
quad7 = c()
quad7[diff7>0 & z7>0] = 1 # H-H
quad7[diff7<0 & z7>0] = 2 # L-L
quad7[diff7>0 & z7<0] = 3 # H-L
quad7[diff7<0 & z7<0] = 4 # L-H
quad7[LISA7[, 5]>0.05]=5 # 不顯著，設定雙尾所以用0.05比較就可以
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(TWN, main = "LISA Map",xlim = c(271884.7,359423.0),ylim = c(2727079,2805226),border = "grey")
plot(bf, border="grey", col=colors[quad7], main = "LISA Map",add = T)
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
#LISA of both
diff8 = bf$A3 - mean(bf$A3) # diff看自己和平均比起來算是H還是L
z8 = LISA8[,4]
quad8 = c()
quad8[diff8>0 & z8>0] = 1 # H-H
quad8[diff8<0 & z8>0] = 2 # L-L
quad8[diff8>0 & z8<0] = 3 # H-L
quad8[diff8<0 & z8<0] = 4 # L-H
quad8[LISA8[, 5]>0.1]=5 # 不顯著，設定雙尾所以用0.05比較就可以
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(TWN, main = "LISA Map",xlim = c(271884.7,359423.0),ylim = c(2727079,2805226),border = "grey")
plot(bf, border="grey", col=colors[quad8], main = "LISA Map of both (unit: district)",xlim = c(271884.7,359423.0),ylim = c(2727079,2805226), add = T)
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
##########################
#LISA(unit: per district)#
##########################
TWN.nb = poly2nb(TWN)
TWN.nb.w = nb2listw(TWN.nb,zero.policy=T)
LISA3 = localmoran(TWN$A1, TWN.nb.w,zero.policy = T,alternative = "two.sided")
LISA4 = localmoran(TWN$A2, TWN.nb.w,zero.policy = T,alternative = "two.sided")
LISA5 = localmoran(TWN$A3, TWN.nb.w,zero.policy = T,alternative = "two.sided")
#LISA of A1
diff3 = TWN$A1 - mean(TWN$A1) # diff看自己和平均比起來算是H還是L
z3 = LISA3[,4]
quad3 = c()
quad3[diff3>0 & z3>0] = 1 # H-H
quad3[diff3<0 & z3>0] = 2 # L-L
quad3[diff3>0 & z3<0] = 3 # H-L
quad3[diff3<0 & z3<0] = 4 # L-H
quad3[LISA3[, 5]>0.1]=5 # 不顯著，設定雙尾所以用0.05比較就可以
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(TWN, border="grey", col=colors[quad3], main = "LISA Map of A1 (unit: district)",xlim = c(271884.7,359423.0),ylim = c(2727079,2805226))
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
#LISA of A2
diff4 = TWN$A2 - mean(TWN$A2) # diff看自己和平均比起來算是H還是L
z4 = LISA4[,4]
quad4 = c()
quad4[diff4>0 & z4>0] = 1 # H-H
quad4[diff4<0 & z4>0] = 2 # L-L
quad4[diff4>0 & z4<0] = 3 # H-L
quad4[diff4<0 & z4<0] = 4 # L-H
quad4[LISA4[, 5]>0.1]=5 # 不顯著，設定雙尾所以用0.05比較就可以
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(TWN, border="grey", col=colors[quad4], main = "LISA Map of A2 (unit: district)",xlim = c(271884.7,359423.0),ylim = c(2727079,2805226))
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
#LISA of combination
diff5 = TWN$A3 - mean(TWN$A3) # diff看自己和平均比起來算是H還是L
z5 = LISA5[,4]
quad5 = c()
quad5[diff5>0 & z5>0] = 1 # H-H
quad5[diff5<0 & z5>0] = 2 # L-L
quad5[diff5>0 & z5<0] = 3 # H-L
quad5[diff5<0 & z5<0] = 4 # L-H
quad5[LISA5[, 5]>0.1]=5 # 不顯著，設定雙尾所以用0.05比較就可以
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(TWN, border="grey", col=colors[quad5], main = "LISA Map of both (unit: district)",xlim = c(271884.7,359423.0),ylim = c(2727079,2805226))
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
global moran’s I
M=moran.test(bf$A3,listw=bf.nb.w, zero.policy=T)
mc=moran.mc(bf$A3,listw=bf.nb.w,nsim=999,zero.policy=T)
hist(mc$res)
abline(v=M$estimate[1], col="red")
moran.plot (bf$A3, bf.nb.w, zero.policy=T)
cor=sp.correlogram(bf.nb, bf$A3, order=5, method="I", style="W",zero.policy=T)
print(cor); plot(cor)## Spatial correlogram for bf$A3
## method: Moran's I
## estimate expectation variance standard deviate Pr(I) two sided
## 1 (327) 0.8621161 -0.0030675 0.0023973 17.6703 < 2.2e-16
## 2 (198) 0.4050538 -0.0050761 0.0041673 6.3533 2.108e-10
## 3 (99) 0.1950512 -0.0102041 0.0079590 2.3007 0.02141
## 4 (69) 0.1270259 -0.0147059 0.0114233 1.3261 0.18481
## 5 (45) 0.0659295 -0.0227273 0.0160179 0.7005 0.48361
##
## 1 (327) ***
## 2 (198) ***
## 3 (99) *
## 4 (69)
## 5 (45)
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1