RIP

讀取資料

library(rgdal);library(spatstat);library(GISTools);library(aspace)
accidentA2<-read.csv("/Users/yangyuchen/Desktop/空間分析/107年度A2交通事故資料.csv",header = T,stringsAsFactors = FALSE)
accidentA1<-read.csv("/Users/yangyuchen/Desktop/空間分析/107年度A1交通事故資料.csv",header = T,stringsAsFactors = FALSE)
RIP<-read.csv("/Users/yangyuchen/Desktop/空間分析/北北基殯葬設施 - 工作表1.csv",header=T,stringsAsFactors = FALSE)
RIP<-RIP[,1:10]
RIP<-na.omit(RIP)
accidentA1<-accidentA1[(substring(accidentA1$發生地點,1,3)=="臺北市"|substring(accidentA1$發生地點,1,3)=="新北市"|substring(accidentA1$發生地點,1,3)=="基隆市"),]
accidentA1$受傷人數<-substring(accidentA1$死亡受傷人數,7,8)
accidentA1$月份<-substring(accidentA1$發生時間,5,6)
accidentA2<-accidentA2[(substring(accidentA2$發生地點,1,3)=="臺北市"|substring(accidentA2$發生地點,1,3)=="新北市"|substring(accidentA2$發生地點,1,3)=="基隆市"),]
accidentA2$受傷人數<-substring(accidentA2$死亡受傷人數,7,8)
accidentA2$死亡人數<-substring(accidentA2$死亡受傷人數,3,3)
accidentA2$月份<-substring(accidentA2$發生時間,5,6)
TW <- readOGR(dsn = "/Users/yangyuchen/Desktop/空間分析/Data2", layer = "Popn_TWN2", encoding="big5",verbose = F)
TWN<-TW[(TW$COUNTY=="臺北市"|TW$COUNTY=="基隆市"|TW$COUNTY=="新北市"),]
A1<-SpatialPointsDataFrame(cbind(accidentA1$經度,accidentA1$緯度),data = accidentA1,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
A1<-spTransform(A1,proj4string(TWN))
A2<-SpatialPointsDataFrame(cbind(accidentA2$經度,accidentA2$緯度),data = accidentA2,proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
A2<-spTransform(A2,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))

par(mar = c(0,0,1,0),family = "黑體-繁 中黑")
plot(TWN, border="darkgrey",col="whitesmoke",main="交通事故和殯葬設施分佈圖",xlim=c(281209.5,349718.8),ylim=c(2727013,2802321))
plot(RIP,col="green",pch=20,add=T,cex=0.5)
plot(A2,col="blue",pch=20,add=T,cex=0.25)
plot(A1,col="red",pch=20,add=T,cex=0.25)
legend("bottomright",legend = c("殯葬設施","受傷車禍","死亡車禍"),col = c("green","blue","red"),pch=20)

nearest(依據類別做 F function & Bivariate K-function)

linguta<-RIP[RIP$類別=="骨灰(骸)存放設施",]
gonmu<-RIP[RIP$類別=="公墓",]
Windows <- as.owin(TWN)
RIP_ppp<-ppp(RIP@coords[,1],RIP@coords[,2],Windows)
linguta_ppp<-ppp(linguta@coords[,1],linguta@coords[,2],Windows)
gonmu_ppp<-ppp(gonmu@coords[,1],gonmu@coords[,2],Windows)
A1_ppp<-ppp(A1@coords[,1],A1@coords[,2],Windows)
A2_ppp<-ppp(A2@coords[,1],A2@coords[,2],Windows)
all_accident<-data.frame(c(A2@coords[,1],A1@coords[,1]), c(A2@coords[,2],A1@coords[,2]),c(A1$發生時間,A2$發生時間))
all_accident_ppp<-ppp(all_accident$c.A2.coords...1...A1.coords...1..,all_accident$c.A2.coords...2...A1.coords...2..,Windows)

f-function
all:不同設施比較

par(mfrow=c(1,2))
nnd<-nncross(linguta_ppp,all_accident_ppp)
plot(ecdf(nnd$dist),main="linguta f-function")
abline(v=c(500,1000,1500),col="red")
nnd<-nncross(gonmu_ppp,all_accident_ppp)
plot(ecdf(nnd$dist),xlim=c(0,2500),main="gonmu f-function")
abline(v=c(500,1000,1500),col="red")

A1:不同設施比較

par(mfrow=c(1,2))
nnd<-nncross(linguta_ppp,A1_ppp)
plot(ecdf(nnd$dist),main="linguta f-function")
abline(v=c(1000,2000,4000),col="red")
nnd<-nncross(gonmu_ppp,A1_ppp)
plot(ecdf(nnd$dist),xlim=c(0,7000),main="gonmu f-function")
abline(v=c(1000,2000,4000),col="red")

蒙地卡羅

nnd<-nncross(all_accident_ppp,RIP_ppp,k=3)
nnd.r<-data.frame(row.names = seq(1:(nrow(all_accident)-100)))
for(i in 1:100){
  .<-rpoint(nrow(RIP), win = Windows)
  r<-nncross(all_accident_ppp,.,k = 3)
  nnd.r[,i]<-r$dist
}
plot(ecdf(nnd$dist),main="F-function(k=3)",xlab="distance",col="red",lwd="3")
for(i in 1:100){
  G1 = ecdf(nnd.r[,i])
  lines(G1,col='#00000044')
}

最初紅線落在模擬的區間之上，而後隨即沒入包絡之間。由此現象，可以推估是有部分殯葬設施和車禍十分群聚，但有部分殯葬設施和車禍相鄰甚遠，因此，我們需要再更進一步找出這些群聚的殯葬設施。

K-bivariate & L-bivariate

a<-data.frame(c(RIP@coords[,1],A2@coords[,1],A1@coords[,1]),c(RIP@coords[,2],A2@coords[,2],A1@coords[,2]))
a[1:327,3]<-"RIP"
a[328:nrow(a),3]<-"accident"
all<-ppp(a$c.RIP.coords...1...A2.coords...1...A1.coords...1..,a$c.RIP.coords...2...A2.coords...2...A1.coords...2..,Windows,marks = factor(a$V3))

all_L<-(Lcross(all,"RIP","accident"))
all_K<-(Kcross(all,"RIP","accident"))
plot(all_L,main="L-cross function")

plot(all_K,main="K-cross function")