空間分析 期末報告
楊鎧綾 甘佳昀 盧淨婕
library(rgdal);library(aspace);library(GISTools);library(spatstat);library(spdep)
library(ggplot2);library(splancs);library(raster);library(dplyr);library(ggtern)
library(sp);library(classInt);library(ggExtra);library(RColorBrewer)#畫圖畫顏色的套件#讀資料
Popn.TWN <- readOGR(dsn = "E:/COURSE/SPATIAL_analysis/report2/taiwan_li", layer = "VILLAGE_MOI_121_1070312", encoding="unicode",verbose = F)
Popn.TWN2 <- readOGR(dsn = "E:/COURSE/SPATIAL_analysis/week3/Data2", layer = "Popn_TWN2", encoding="unicode",verbose = F)
rnb<-read.csv("E:/COURSE/SPATIAL_analysis/report2/RNB.csv")
hotel<-read.csv("E:/COURSE/SPATIAL_analysis/report2/hotel.csv")
travel<-read.csv("E:/COURSE/SPATIAL_analysis/report2/travel.csv")
#資料轉換成shp
newtai<-Popn.TWN[Popn.TWN$COUNTYNAME=="新北市",]
hotel.shp <- SpatialPointsDataFrame(cbind(hotel$Response_X,hotel$Response_Y), hotel ,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"))
rnb.shp <- SpatialPointsDataFrame(cbind(rnb$Response_X,rnb$Response_Y), rnb ,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"))
travel.shp <- SpatialPointsDataFrame(cbind(travel$Px,travel$Py), travel ,proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84"))
travel.shp<-spTransform(travel.shp,Popn.TWN@proj4string)用面量圖看分布
#面量圖前處理
#計算出每個村里的平均最低價格
hol.inter<-gIntersection(newtai,hotel.shp,byid = T)
rnb.inter<-gIntersection(newtai,rnb.shp,byid = T)
tmp_1 <- strsplit(rownames(data.frame(hol.inter)), " ")
li_ID<-unlist(lapply(tmp_1, function(x) x[1]))
li_ID<-as.numeric(li_ID)
hol_ID<-unlist(lapply(tmp_1, function(x) x[2]))
hol_ID<-as.numeric(hol_ID)
hol.attr=data.frame(li_ID, hol_ID )#旅館分段id處理完畢
tmp_2 <- strsplit(rownames(data.frame(rnb.inter)), " ")
li_ID.r<-unlist(lapply(tmp_2, function(x) x[1]))
li_ID.r<-as.numeric(li_ID.r)
rnb_ID<-unlist(lapply(tmp_2, function(x) x[2]))
rnb_ID<-as.numeric(rnb_ID)
rnb.attr=data.frame(li_ID.r, rnb_ID )#民宿分段id處理完畢
for (i in 1:length(li_ID)) {
hol.attr$h_price[i]<-hotel$button_price[hol_ID[i]] #旅館價格
}
for (i in 1:length(li_ID.r)) {
rnb.attr$r_price[i]<-rnb$button_price[rnb_ID[i]] #民宿價格
}
#樞紐分析
rp.tabs<-data.frame(xtabs(rnb.attr$r_price~rnb.attr$li_ID.r))#民宿房間數完成
hp.tabs<-data.frame(xtabs(hol.attr$h_price~hol.attr$li_ID))#旅館房間數完成
colnames(rp.tabs)=c("li_id","rnb_price")
colnames(hp.tabs)=c("li_id","hol_price")
#把價格結果放回newtai裡
newtai$li_id<-rownames(newtai@data)
newtai@data<- left_join(newtai@data, hp.tabs)
newtai@data<- left_join(newtai@data, rp.tabs)
#把na用0取代
newtai$rnb_price[is.na(newtai$rnb_price)]=0
newtai$hol_price[is.na(newtai$hol_price)]=0#各里平均最低房價面量圖
lm.palette=colorRampPalette(c( "white","orange","red"), space = "rgb")
spplot(newtai, zcol="hol_price", col.regions=lm.palette(20), main="各里旅館平均最低價格面量圖")
spplot(newtai, zcol="hol_price", col.regions=lm.palette(20), main="各里旅館平均最低價格面量圖(局部放大)",xlim=c( 285113.3,306783.2),ylim = c(2760155,2777151))
spplot(newtai, zcol="rnb_price", col.regions=lm.palette(20), main="各里民宿平均最低價格面量圖")
面LISA
TWN_nbh<-poly2nb(newtai)# Defining Spatial Neighbors
TWN_nbh_w<- nb2listw(TWN_nbh, zero.policy=T)# 2. Building Weighting Matrix
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
hol_price<-newtai$hol_price
LISA.h <- localmoran(hol_price, TWN_nbh_w, zero.policy=T,alternative = 'two.sided')
diff = hol_price - mean(hol_price) # diff看自己和平均比起來算是H還是L
z = LISA.h[,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.h[, 5]>0.05]=5 # 不顯著，設定雙尾所以用0.05比較就可以
plot(newtai, border="grey", col=colors[quad], main = "各里旅館平均最低價格LISA圖",xlim=c(272202.2,345783.1),ylim = c(2729675,2807715))
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
plot(newtai, border="grey", col=colors[quad], main = "各里旅館平均最低價格LISA圖(局部放大)",xlim=c( 285113.3,306783.2),ylim = c(2760155,2777151))
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
hol_price<-newtai$rnb_price
LISA.h <- localmoran(hol_price, TWN_nbh_w, zero.policy=T,alternative = 'two.sided')
diff = hol_price - mean(hol_price) # diff看自己和平均比起來算是H還是L
z = LISA.h[,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.h[, 5]>0.05]=5 # 不顯著，設定雙尾所以用0.05比較就可以
plot(newtai, border="grey", col=colors[quad], main = "各里民宿平均最低價格LISA",xlim=c(272202.2,345783.1),ylim = c(2729675,2807715))
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
點LISA
#空間自相關
price = hotel$button_price
coords = hotel.shp@coords
#定義鄰近
TW.nb = dnearneigh(coords, d1=0, d2=10000)
#建立鄰近目錄
TW.nb.w = nb2listw(TW.nb,zero.policy=T) #預設style="W"(列標準化)
#lisa
LISA = localmoran(price, TW.nb.w,zero.policy = T,alternative = "two.sided")
#區分顏色
diff = price - mean(price) # diff看自己和平均比起來算是H還是L
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 # 不顯著，設定雙尾所以用0.05比較就可以
#畫圖
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(newtai, border="grey", main = "LISA Map of hotel")
points(rnb.shp,col=colors[quad],pch=16)
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
#空間自相關
price = rnb$button_price
coords = rnb.shp@coords
#定義鄰近
TW.nb = dnearneigh(coords, d1=0, d2=10000)
#建立鄰近目錄
TW.nb.w = nb2listw(TW.nb,zero.policy=T) #預設style="W"(列標準化)
#lisa
LISA = localmoran(price, TW.nb.w,zero.policy = T,alternative = "two.sided")
#區分顏色
diff = price - mean(price) # diff看自己和平均比起來算是H還是L
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 # 不顯著，設定雙尾所以用0.05比較就可以
#畫圖
colors=c("red", "blue", "lightpink", "skyblue2", rgb(.95, .95, .95))
plot(newtai, border="grey", main = "LISA Map of RNB")
points(rnb.shp,col=colors[quad],pch=16)
legend("bottomright",legend=c("HH","LL","HL","LH","NS"),fill=colors,bty="n",cex=0.7,y.intersp=1,x.intersp=1)
f-function(旅館是否群聚在觀光景點旁)
#ppp
Windows.h= as.owin(newtai)
#hotel
x.coor.h<-hotel.shp@coords[,1]
y.coor.h<-hotel.shp@coords[,2]
hotel.ppp<-ppp(x.coor.h,y.coor.h,Windows.h)
#travel
x.coor.t<-travel.shp@coords[,1]
y.coor.t<-travel.shp@coords[,2]
travel.ppp<-ppp(x.coor.t,y.coor.t,Windows.h)
#rnb
x.coor.r<-rnb.shp@coords[,1]
y.coor.r<-rnb.shp@coords[,2]
rnb.ppp<-ppp(x.coor.r,y.coor.r,Windows.h)
#f-function
##旅館
nnd.ht=nncross(hotel.ppp, travel.ppp)
hotel$nei_di<-nnd.ht$dist
F.ht = ecdf(nnd.ht$dist)
plot(F.ht,cex=0,verticals=T,xaxs="i",yaxs="i",col='red',main="新北市旅館 V.S. 觀光景點")
for(i in 1:100){
random.h<-rpoint(844, win=Windows.h)
nn1.ht<-nncross(random.h,travel.ppp)
G.ran.ht = ecdf(nn1.ht$dist)
lines(G.ran.ht, main="F function",xlim=c(0,5000),cex=0,verticals=T,xaxs="i",yaxs="i")
}
lines(F.ht,cex=0,verticals=T,xaxs="i",yaxs="i",col='red')
##民宿
nnd.rt=nncross(rnb.ppp, travel.ppp)
rnb$nei_di<-nnd.rt$dist
F.rt = ecdf(nnd.rt$dist)
plot(F.rt,cex=0,verticals=T,xaxs="i",yaxs="i",col='red',main="新北市民宿 V.S. 觀光景點")
for(i in 1:100){
random.r<-rpoint(nrow(rnb.shp@data), win=Windows.h)
nn1.rt<-nncross(random.r,travel.ppp)
G.ran.rt = ecdf(nn1.rt$dist)
lines(G.ran.rt, main="F function",xlim=c(0,5000),cex=0,verticals=T,xaxs="i",yaxs="i")
}
lines(F.rt,cex=0,verticals=T,lwd=5,xaxs="i",yaxs="i",col='red')
相關性與回歸(與觀光景點的距離和房價的相關性)
cor.test(hotel$nei_di,hotel$button_price)##
## Pearson's product-moment correlation
##
## data: hotel$nei_di and hotel$button_price
## t = 1.9016, df = 250, p-value = 0.05837
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.004228508 0.239446446
## sample estimates:
## cor
## 0.1194067cor.test(rnb$nei_di,rnb$button_price)##
## Pearson's product-moment correlation
##
## data: rnb$nei_di and rnb$button_price
## t = 0.88983, df = 247, p-value = 0.3744
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.06826792 0.17958247
## sample estimates:
## cor
## 0.05652815results<-lm(hotel$button_price~hotel$nei_di)
plot(hotel$button_price~hotel$nei_di,pch=16,cex=1,col="navy",main="旅館最近觀光景點距離與最低房價回歸圖",xlab="near_distance",ylab="price")
abline(results,col="red")
results<-lm(rnb$button_price~rnb$nei_di)
plot(rnb$button_price~rnb$nei_di,pch=16,cex=1,col="navy",main="民宿最近觀光景點距離與最低房價回歸圖",xlab="near_distance",ylab="price")
abline(results,col="red")