setwd("C:/Users/Wu/Desktop/107-2 course/spatial analysis/project")
library(rgdal)
library(GISTools)
library(ggtern)
library(splancs)
library(raster)
library(sp)
accidentA1<-read.csv("./Data/107年度A1交通事故資料.csv",header = T,stringsAsFactors = FALSE)
accidentA2<-read.csv("./Data/107年度A2交通事故資料.csv",header = T,stringsAsFactors = FALSE)
RIP<-read.csv("./Data/北北基殯葬設施.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)
TW <- readOGR(dsn = "./Data/Popn_TWN2.shp", layer = "Popn_TWN2", encoding="utf8", verbose = F)
TWN<-TW[(TW$COUNTY=="臺北市"|TW$COUNTY=="基隆市"|TW$COUNTY=="新北市"),]
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))
#RIP<-na.omit(RIP)
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))
純kde
#殯葬設施
Kde.RIP = kde.points(SpatialPoints(cbind(RIP$coords.x1, RIP$coords.x2)), h=2000,n = 1000, lims = TWN)
Kde.R1=raster(Kde.RIP)
plot(Kde.R1, col = brewer.pal(9,"YlGn"),xlim=c(269460.4,362639.3),ylim=c(2721067,2808525)) #KDE圖
masker=poly.outer(Kde.R1, TWN) #建立遮罩
## Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
## unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4
## strings
add.masking(masker, col="white") #覆蓋遮罩
plot(TWN, add=T)
#交通事故
Kde.acc = kde.points(SpatialPoints(cbind(c(A1$coords.x1, A2$coords.x1), c(A1$coords.x2, A2$coords.x2))), h=2000,lims = TWN,n = 1000)
Kde.R2=raster(Kde.acc)
plot(Kde.R2, col = brewer.pal(9,"YlGn"),xlim=c(269460.4,362639.3),ylim=c(2721067,2808525)) #KDE圖
masker=poly.outer(Kde.R2, TWN) #建立遮罩
## Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
## unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4
## strings
add.masking(masker, col="white") #覆蓋遮罩
plot(TWN, add=T)
weighted kde
RIP.buf <- gBuffer(RIP, width = 1000, byid = T, id = rownames(RIP@data))
A1.count <- poly.counts(A1, RIP.buf) #big=which(A1.count>3);RIP@data$名稱[big]
A2.count <- poly.counts(A2, RIP.buf)
RIP$A1count=A1.count;RIP$A2count=A2.count
weights = RIP$A1count+RIP$A2count
#以事故數量加權殯葬設施
Kde.lim=c(TWN@bbox[1,], TWN@bbox[2,])
Kde.w = kde2d.weighted(RIP$coords.x1, RIP$coords.x2, 2000, 1000, lims = Kde.lim, weights) #kde
image(Kde.w, asp=1, main = "KDE", col = brewer.pal(9,"YlGn"),xlim=c(269460.4,362639.3),ylim=c(2721067,2808525))
masker = poly.outer(as.points(TWN@bbox[1,], TWN@bbox[2,]), TWN, extend = 300)
## Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
## unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4
## strings
add.masking(masker, col="white")
plot(TWN, add=T)
dual kde
### 事故和設施 ###
kde1 = kde2d(RIP$coords.x1, RIP$coords.x2, 2000, 1000, Kde.lim) #確認範圍一致
kde2 = kde2d(c(A1$coords.x1, A2$coords.x1), c(A1$coords.x2, A2$coords.x2), 2000, 1000, Kde.lim)
# 設施減事故
diff = kde2$z - kde1$z
kde.diff = list(x = kde1$x, y = kde1$y, z = diff)
image(kde.diff, asp=1, main = "Dual KDE diff", col = brewer.pal(9,"RdYlBu"),xlim=c(269460.4,362639.3),ylim=c(2721067,2808525))
masker = poly.outer(as.points(TWN@bbox[1,], TWN@bbox[2,]), TWN, extend = 300)
## Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
## unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4
## strings
add.masking(masker, col="white")
plot(TWN, add=T)
# 設施加事故(不知道可不可以)
add = kde1$z + kde2$z
kde.add = list(x = kde1$x, y = kde1$y, z = add)
image(kde.add, asp=1, main = "Dual KDE add", col = brewer.pal(9,"RdYlBu"),xlim=c(269460.4,362639.3),ylim=c(2721067,2808525))
masker = poly.outer(as.points(TWN@bbox[1,], TWN@bbox[2,]), TWN, extend = 300)
## Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
## unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4
## strings
add.masking(masker, col="white")
plot(TWN, add=T)
