空間分析期末報告：R code
陳竑屹 林柏翰 陳奕儒
資料整理
library(GISTools)
library(rgdal)
library(spatstat)
library(dplyr)
library(sp)
library(spdep)
setwd('C:/Users/User/Desktop/大二下課堂/空間分析/空間分析期末報告') #改成所在的路徑
z1 = read.csv('Address_Finish-八大場所.csv',header = T)
z2 = read.csv('Address_Finish-國高中.csv',header = T)
z3 = read.csv('Address_Finish-幼稚園.csv',header = T)
z4 = read.csv('Address_Finish-電玩.csv',header = T)
z5 = read.csv('Address_Finish-網咖.csv',header = T)
maps = readOGR(dsn = ".", layer = "Taipei_Vill", encoding="utf8", verbose = F) #TWD97
spuse = readOGR(dsn = ".", layer = "toshin", encoding="utf8") #WGS84
# 八大場所
big8 = SpatialPointsDataFrame(cbind(z1$Response_X,z1$Response_Y),data=z1,proj4string = CRS(proj4string(spuse)))
big8 = spTransform(big8,maps@proj4string)
# 學校
school = SpatialPointsDataFrame(cbind(z2$Response_X,z2$Response_Y),data=z2,proj4string = CRS(proj4string(spuse)))
school = spTransform(school,maps@proj4string)
# 幼稚園
kin = SpatialPointsDataFrame(cbind(z3$Response_X,z3$Response_Y),data=z3,proj4string = CRS(proj4string(spuse)))
kin = spTransform(kin,maps@proj4string)
# 電玩
play = SpatialPointsDataFrame(cbind(z4$Response_X,z4$Response_Y),data=z4,proj4string = CRS(proj4string(spuse)))
play = spTransform(play,maps@proj4string)
# 網咖
incafe = SpatialPointsDataFrame(cbind(z5$Response_X,z5$Response_Y),data=z5,proj4string = CRS(proj4string(spuse)))
incafe = spTransform(incafe,maps@proj4string)1.台北市不可設置不良場所的地方占整體比例（法規合不合理）
school.range=gBuffer(school,width=100,byid=F) #模擬校地範圍
school.range2=gBuffer(school.range,width=200,byid=F) #國小以上學校200m Buffer
kin.range2=gBuffer(kin,width=200,byid=F) #公私立幼兒園200m Buffer
maps.tpe=gUnaryUnion(maps,id=maps$TOWN)
merge=gUnion(school.range2,kin.range2,byid=F) #不可設置區域
percent=(poly.areas(merge)/sum(poly.areas(maps.tpe))) #台北市不可設置不良場所的地方占整體比例
areapercent=poly.areas(gIntersection(maps.tpe,merge,byid=T))/poly.areas(maps.tpe) #各區比例
names(areapercent)=names(maps.tpe)
percent 1
0.2846822 areapercent 士林區 大同區 大安區 中山區 中正區 內湖區 文山區
0.1726786 0.6586509 0.6999804 0.4473994 0.6693791 0.2462723 0.3054433
北投區 松山區 信義區 南港區 萬華區
0.1468541 0.5440540 0.4543147 0.1636997 0.5284004 #範圍圖
par(mar = c(1,1,1,1))
plot(maps.tpe)
plot(merge,col='grey',add=T)
legend(315291,2788505,'限制開業範圍',cex=.8,fill='grey')
plot(big8,add=T,col='red',pch=20)
plot(incafe,add=T,col='gold',pch=20)
plot(play,add=T,col='lightpink',pch=20)
legend(315291,2782505,c('八大行業','網咖','電子遊樂場'),cex=.8,col=c('red','gold','lightpink'),pch=19,fill='white',border='white')
2.學校200m內有無不良場所的趨勢：F(d) Function
TPE = gUnaryUnion(maps)
TPE.Coord = TPE@polygons[[1]]@Polygons[[1]]@coords
x = rev(TPE.Coord[,1])
y = rev(TPE.Coord[,2])
TPE.Coord = cbind(x, y)
PTS_bnd = owin(poly = TPE.Coord)
bads = rbind(big8@coords[,1:2], incafe@coords[,1:2], play@coords[,1:2])
bads.ppp = ppp(bads[,1], bads[,2], window = PTS_bnd)
big8.ppp = ppp(big8@coords[,1], big8@coords[,2], window = PTS_bnd)
incafe.ppp = ppp(incafe@coords[,1], incafe@coords[,2], window = PTS_bnd)
play.ppp = ppp(play@coords[,1], play@coords[,2], window = PTS_bnd)
kin.ppp = ppp(kin@coords[,1], kin@coords[,2], window = PTS_bnd)
school.ppp = ppp(school@coords[,1], school@coords[,2], window = PTS_bnd)
nnd.a1 = nncross(bads.ppp, school.ppp, k = 1)
nnd.b1 = nncross(big8.ppp, school.ppp, k = 1)
nnd.c1 = nncross(incafe.ppp, school.ppp, k = 1)
nnd.p1 = nncross(play.ppp, school.ppp, k = 1)
color = rainbow(4)
plot(ecdf(nnd.b1[,1]), col = color[1],
main = "不良場所與最近學校之距離", xlab = "距離（m）", ylab = "F(d)",
xaxs = "i", yaxs = "i", cex = 0, verticals = T, xlim = c(0, 500))
plot(ecdf(nnd.c1[,1]), col = color[2], cex = 0, verticals = T , add = T)
plot(ecdf(nnd.p1[,1]), col = color[3], cex = 0, verticals = T , add = T)
plot(ecdf(nnd.a1[,1]), col = color[4], cex = 0, verticals = T , add = T)
abline(v = 300, col = "lightgrey", lty = 2)
legend("bottomright", cex = 0.7, c("八大場所", "網咖", "電玩", "所有"), fill = color, bty = 'n', y.intersp = 1)
nnd.a2 = nncross(bads.ppp, kin.ppp, k = 1)
nnd.b2 = nncross(big8.ppp, kin.ppp, k = 1)
nnd.c2 = nncross(incafe.ppp, kin.ppp, k = 1)
nnd.p2 = nncross(play.ppp, kin.ppp, k = 1)
plot(ecdf(nnd.b2[,1]), col = color[1],
main = "不良場所與最近幼稚園之距離", xlab = "距離（m）", ylab = "F(d)",
xaxs = "i", yaxs = "i", cex = 0, verticals = T, xlim = c(0, 500))
plot(ecdf(nnd.c2[,1]), col = color[2], cex = 0, verticals = T , add = T)
plot(ecdf(nnd.p2[,1]), col = color[3], cex = 0, verticals = T , add = T)
plot(ecdf(nnd.a2[,1]), col = color[4], cex = 0, verticals = T , add = T)
abline(v = 200, col = "lightgrey", lty = 2)
legend("bottomright", cex = 0.7, c("八大場所", "網咖", "電玩", "所有"), fill = color, bty = 'n', y.intersp = 1)
nnd.a3 = nncross(kin.ppp, bads.ppp, k = 1)
nnd.b3 = nncross(kin.ppp, big8.ppp, k = 1)
nnd.c3 = nncross(kin.ppp, incafe.ppp, k = 1)
nnd.p3 = nncross(kin.ppp, play.ppp, k = 1)
plot(ecdf(nnd.b3[,1]), col = color[1],
main = "幼稚園與最近不良場所之距離", xlab = "距離（m）", ylab = "F(d)",
xaxs = "i", yaxs = "i", cex = 0, verticals = T, xlim = c(0, 500))
plot(ecdf(nnd.c3[,1]), col = color[2], cex = 0, verticals = T , add = T)
plot(ecdf(nnd.p3[,1]), col = color[3], cex = 0, verticals = T , add = T)
plot(ecdf(nnd.a3[,1]), col = color[4], cex = 0, verticals = T , add = T)
abline(v = 200, col = "lightgrey", lty = 2)
legend("bottomright", cex = 0.7, c("八大場所", "網咖", "電玩", "所有"), fill = color, bty = 'n', y.intersp = 1)
nnd.a4 = nncross(school.ppp, bads.ppp, k = 1)
nnd.b4 = nncross(school.ppp, big8.ppp, k = 1)
nnd.c4 = nncross(school.ppp, incafe.ppp, k = 1)
nnd.p4 = nncross(school.ppp, play.ppp, k = 1)
color = rainbow(4)
plot(ecdf(nnd.b4[,1]), col = color[1],
main = "學校與最近不良場所之距離", xlab = "距離（m）", ylab = "F(d)",
xaxs = "i", yaxs = "i", cex = 0, verticals = T, xlim = c(0, 500))
plot(ecdf(nnd.c4[,1]), col = color[2], cex = 0, verticals = T , add = T)
plot(ecdf(nnd.p4[,1]), col = color[3], cex = 0, verticals = T , add = T)
plot(ecdf(nnd.a4[,1]), col = color[4], cex = 0, verticals = T , add = T)
abline(v = 300, col = "lightgrey", lty = 2)
legend("bottomright", cex = 0.7, c("八大場所", "網咖", "電玩", "所有"), fill = color, bty = 'n', y.intersp = 1)
3.面量圖、LISA圖
#面量圖
school.ch=gBuffer(school,width=300,byid=T)
kin.ch=gBuffer(kin,width=200,byid=T)
zall.x=c(z1$Response_X,z4$Response_X,z5$Response_X);zall.y=c(z1$Response_Y,z4$Response_Y,z5$Response_Y)
all=SpatialPointsDataFrame(cbind(zall.x,zall.y),data=data.frame(test=c(1:length(zall.x))),proj4string = CRS(proj4string(spuse)))
all=spTransform(all,maps@proj4string)
all.school=poly.counts(all,school.ch);all.kin=poly.counts(all,kin.ch)
school.qua=school[all.school>0,];kin.qua=kin[all.kin>0,]
area.count=(poly.counts(school.qua,maps.tpe)+poly.counts(kin.qua,maps.tpe))/(poly.counts(school,maps.tpe)+poly.counts(kin,maps.tpe))
shades = auto.shading(area.count, n = 5, cols = brewer.pal(5, "YlOrRd"))
choropleth(maps.tpe, area.count, shades) #附近有不良場所的學校在該區占所有學校的比例
choro.legend(315291,2782505,shades,cex=.8)
#LISA圖
Tpnb=poly2nb(maps.tpe)
Tpnb.w = nb2listw(Tpnb,zero.policy=T)
LISA = localmoran(area.count,Tpnb.w,zero.policy = T,alternative = "two.sided")
diff=area.count-mean(area.count)
z = LISA[,4]
quad = c()
quad[diff>0 & z>0] = 1 # H-H
quad[diff<0 & z>0] = 2 # L-L
quad[diff>0 & z<0] = 3 # H-L
quad[diff<0 & z<0] = 4 # L-H
quad[LISA[, 5]>0.05]=5 # 不顯著
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(maps.tpe, border="grey", col=colors[quad], main = "LISA Map")
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)