空間分析 期末報告

讀取資料

台灣鄉鎮區底圖

rm(list = ls())
library(aspace);library(rgdal);library(spdep);library(GISTools);library(aspace);library(spatstat);library(choroplethr)
TWPOP <- readOGR(dsn="D:/school/Spatial Analysis/data/Popn_TWN2/Popn_TWN2.shp", layer = "Popn_TWN2", encoding="utf8")
TWPOP.LongLat <- spTransform(TWPOP, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
proj=TWPOP.LongLat@proj4string
TWPOP.d <- TWPOP@data

i郵箱

ibox <- read.csv("D:/school/Spatial Analysis/data/i郵箱位置.csv")
XCOOR.i = ibox[,8]
YCOOR.i = ibox[,9]
ibox.p <- SpatialPointsDataFrame(coords=cbind(x=XCOOR.i,y=YCOOR.i), data=ibox, proj4string=proj)

便利超商-選出7-11與全家

store <- read.csv("D:/school/Spatial Analysis/data/簡單版超商位置.csv")
XCOOR.s = store[,2]
YCOOR.s = store[,3]
store.p <- SpatialPointsDataFrame(coords=cbind(x=XCOOR.s,y=YCOOR.s), data=store, proj4string=CRS("+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +units=m +no_defs"))
store.pt<-spTransform(store.p,TWPOP@proj4string)
con.store <-store.p[store.p$store=="seven"|store.p$store=="family",]

郵務局資料(包含全台分局以及投遞總件數、失敗件數)

post.off <- read.csv("D:/school/Spatial Analysis/data/Post_All_new1 - 複製.csv")
BranchMailCount<-read.csv("D:/school/Spatial Analysis/data/BranchMailCount.csv")

total_post<-merge(post.off,BranchMailCount,by="branch")
total_post<-total_post[,-(3:8)]
total_post<-total_post[,-5]
total_post<-total_post[,-(11:22)]#這邊是整理資料，把不相關的資料（電話等等）拿掉
total_post<-total_post[,-13]
total_post<-total_post[-468,] #這個郵局的經緯度異常小，google地圖上也沒有郵局了，所以我把他拔掉

total_post$服務時間<-(total_post$平常日提供服務+total_post$平常日延時提供服務+total_post$週六提供服務+total_post$週日提供服務+total_post$提供郵務業務延時服務郵局.週一至週五+total_post$提供郵務及儲匯業務延時服務郵局.週一至週五)

XCOOR.p = total_post[,5]
YCOOR.p = total_post[,6]
total_post.p <- SpatialPointsDataFrame(coords=cbind(x=XCOOR.p,y=YCOOR.p), data=total_post, proj4string=proj)

分析

投遞失敗率面量圖

#計算各鄉鎮市區投遞失敗率
for (i in 1:length(TWPOP.LongLat)){
  s=total_post.p$all..[total_post.p$鄉鎮市區 %in% TWPOP.LongLat$TOWN[i]]
  TWPOP.LongLat$ALL[i]=sum(s)
  sS=total_post.p$NotAtHome..[total_post.p$鄉鎮市區 %in% TWPOP.LongLat$TOWN[i]]
  TWPOP.LongLat$NotAtHome[i]=sum(sS)
  TWPOP.LongLat$ratio[i]=TWPOP.LongLat$NotAtHome[i]/TWPOP.LongLat$ALL[i]
}
TWPOP.LongLat$ratio[is.na(TWPOP.LongLat$ratio)] <- 0 
#繪製面失敗率量圖
lm.palette <- colorRampPalette(c("white","orange", "red"), space = "rgb")
spplot(TWPOP.LongLat, zcol="ratio", col.regions=lm.palette(20), main="投遞失敗率面量圖")

投遞失敗率-cluster map與i郵箱點位圖

TW_nb<-poly2nb(TWPOP.LongLat) 
Ratio=TWPOP.LongLat$ratio
TW_nb_in<-include.self(TW_nb)
TW_nb_in_w<- nb2listw(TW_nb_in, zero.policy=T)
LG<-localG(Ratio, TW_nb_in_w)
LG1<-0
for (i in 1:length(TWPOP.LongLat)){LG1[i]<-LG[i]}
TWPOP.LongLat$LG<-LG1
lm.palette <- colorRampPalette(c("white","orange", "red"), space = "rgb")
LGV<-TWPOP.LongLat$LG
chk<-LGV-1.65
quadrant <- vector(mode="numeric",length(TWPOP.LongLat))
quadrant[chk>0] <- 1 
quadrant[chk<0] <- 2 
colors <- c("red", "lightgray")
plot(TWPOP.LongLat, border="darkgrey", col=colors[quadrant], main = "Cluster Map: 投遞失敗率")
legend("bottomright",legend=c("Cluster","Non-cluster"), fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
#將i郵箱點上
plot(ibox.p,add=T,pch=16,cex=0.8,col='blue')

便利商店服務人數面量圖

store.count=poly.counts(store.p, TWPOP)
for (i in 1:length(TWPOP.LongLat)){
  TWPOP.LongLat$store.density[i]=(store.count[i]/poly.areas(TWPOP.LongLat)[i])*10^6
}

for (i in 1:length(TWPOP.LongLat)){
  TWPOP.LongLat$store.service[i]=((TWPOP.LongLat$A0A14_CNT[i]+TWPOP.LongLat$A15A64_CNT[i]+TWPOP.LongLat$A65UP_CNT[i])/store.count[i])
}

for (i in 1:length(TWPOP.LongLat)){
  if (TWPOP.LongLat$store.service[i]==Inf){
    TWPOP.LongLat$store.service[i]=0
  }
}

lm.palette <- colorRampPalette(c("white", "green"), space = "rgb") # 便利商店服務人數圖
spplot(TWPOP.LongLat, zcol="store.service", col.regions=lm.palette(20), main="便利商店服務人數圖")

便利商店服務人數LISA與i郵箱點位圖

Ser=TWPOP.LongLat$store.service
TW_nb_w<- nb2listw(TW_nb, zero.policy=T)
LISA.S <- localmoran(Ser, TW_nb_w, zero.policy=T,alternative = 'two.sided')
TWPOP.LongLat$z.li <- LISA.S[,4]
TWPOP.LongLat$pvalue <- LISA.S[,5]
chk<-Ser - mean(Ser)
zi<- LISA.S[,4]
quadrant <-c()
quadrant[chk>0 & zi>0] <- 1 # H-H
quadrant[chk<0 & zi>0] <- 2 # L-L
quadrant[chk>0 & zi<0] <- 3 # H-L
quadrant[chk<0 & zi<0] <- 4 # L-H
quadrant[LISA.S[,5]>0.05]<-5
signif <- 0.05
quadrant[LISA.S[, 5]> signif] <- 5
colors <- c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(TWPOP.LongLat, border="grey", col=colors[quadrant], main = "LISA Cluster Map: 便利商店服務人數")
legend("bottomright",legend=c("High-High","Low-Low","High-Low","Low-High","NS"), fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
#將i郵箱點上
plot(ibox.p,add=T,pch=16,cex=0.8,col='green')