國企五 B03704082 廖威凱

#環境建置與導入資料
rm(list = ls())
library(GISTools)
library(rgdal)
library(aspace)
library(spatstat)
setwd("C:/Spatial Analysis/Final Project/shapefile")

#大阪市資料(依序為:行政區、地鐵站、TOP20景點、便利商店、購物商店、旅店、販賣機)
OSA <- readOGR(dsn = ".", layer = "Osaka", encoding = "UTF-8")
OSA.Metro <- readOGR(dsn = ".", layer = "Osaka Metro", encoding = "UTF-8")
OSA.Place <- readOGR(dsn = ".", layer = "Osaka Attraction", encoding = "UTF-8")
OSA.Con <- readOGR(dsn = ".", layer = "Osaka Convenience Store", encoding = "UTF-8")
OSA.Buy <- readOGR(dsn = ".", layer = "Osaka Shopping", encoding = "UTF-8")
OSA.Stay <- readOGR(dsn = ".", layer = "Osaka Accomodation", encoding = "UTF-8")
OSA.Vend <- readOGR(dsn = ".", layer = "Osaka Vending Machine", encoding = "UTF-8")

#投影轉換(JGD -> JGD2011 / UTM zone 53N)
OSA <- spTransform(OSA ,CRS("+proj=utm +zone=53 +ellps=GRS80 +units=m +no_defs"))
OSA.Metro <- spTransform(OSA.Metro, OSA@proj4string)
OSA.Place <- spTransform(OSA.Place, OSA@proj4string)
OSA.Con <- spTransform(OSA.Con, OSA@proj4string)
OSA.Buy <- spTransform(OSA.Buy, OSA@proj4string)
OSA.Stay <- spTransform(OSA.Stay, OSA@proj4string)
OSA.Vend <- spTransform(OSA.Vend, OSA@proj4string)

#臺北市資料
TPE <- readOGR(dsn = ".", layer = "Taipei", encoding = "UTF-8")
TPE.Metro <- readOGR(dsn = ".", layer = "Taipei Metro", encoding = "UTF-8")
TPE.Place <- readOGR(dsn = ".", layer = "Taipei Attraction", encoding = "UTF-8")
TPE.Con <- readOGR(dsn = ".", layer = "Taipei Convenience Store", encoding = "UTF-8")
TPE.Buy <- readOGR(dsn = ".", layer = "Taipei Shopping", encoding = "UTF-8")
TPE.Stay <- readOGR(dsn = ".", layer = "Taipei Accomodation", encoding = "UTF-8")
  1. 用景點的Median Center決定住宿處
#界定研究範圍
Windows <- as.owin(OSA)

##景點的Median Center
OSA.Place.Median <- median_centre(points = OSA.Place@coords)
##   id median.x median.y
## 1  1 546161.9  3836060
#轉成Spatial Point
OSA.Place.Median.sp <- SpatialPoints(coords = cbind(OSA.Place.Median$median.x, OSA.Place.Median$median.y), OSA@proj4string)
#轉成ppp格式
OSA.Place.Median.ppp <- ppp(OSA.Place.Median$median.x, OSA.Place.Median$median.y, Windows)

#記錄地鐵站座標
OSA.Metro.df <- data.frame(OSA.Metro@coords)
#轉成ppp格式
OSA.Metro.ppp <- ppp(OSA.Metro.df$coords.x1, OSA.Metro.df$coords.x2, Windows)

#找出距離景點中心500m內的地鐵站
for(i in 1:length(OSA.Metro)){
  nnc <- nncross(OSA.Place.Median.ppp, OSA.Metro.ppp, k = i)
  if(nnc[,1] < 500){
    number <- nnc[,2]
    print(OSA.Metro$name_en[number])
  }
}
## [1] Nippombashi
## 98 Levels: Abeno Abiko Asashiobashi Awaza Bentencho ... Zuiko 4-chome
## [1] Namba
## 98 Levels: Abeno Abiko Asashiobashi Awaza Bentencho ... Zuiko 4-chome
## [1] Namba
## 98 Levels: Abeno Abiko Asashiobashi Awaza Bentencho ... Zuiko 4-chome

可知難波(Namba)、日本橋(Nippombashi)這兩站附近地區是適合的住宿地點

  1. 景點周遭環境的便利性
#轉成ppp格式
OSA.Place.df <- data.frame(OSA.Place@coords)
OSA.Place.ppp <- ppp(OSA.Place.df$coords.x1, OSA.Place.df$coords.x2, Windows)

OSA.Con.df <- data.frame(OSA.Con@coords)
OSA.Con.ppp <- ppp(OSA.Con.df$coords.x1, OSA.Con.df$coords.x2, Windows)

OSA.Buy.df <- data.frame(OSA.Buy@coords)
OSA.Buy.ppp <- ppp(OSA.Buy.df$coords.x1, OSA.Buy.df$coords.x2, Windows)

OSA.Vend.df <- data.frame(OSA.Vend@coords)
OSA.Vend.ppp <- ppp(OSA.Vend.df$coords.x1, OSA.Vend.df$coords.x2, Windows)

#景點與地鐵站(紅)、便利商店(黑)、購物商店(藍)、販賣機(綠)的F Function圖
plot(c(0, 1500), c(0, 1), type = 'n', xaxs = "i", yaxs = "i", xlab = "距離", ylab = "F")

for(i in 1:length(OSA.Place)){
  nnc <- nncross(OSA.Place.ppp, OSA.Metro.ppp, k = 1)
  F.func <- ecdf(nnc$dist)
  lines(F.func, col = "red", cex = 0, lwd = 3, verticals = T)
}


for(i in 1:length(OSA.Place)){
  nnc <- nncross(OSA.Place.ppp, OSA.Con.ppp, k = 1)
  F.func <- ecdf(nnc$dist)
  lines(F.func, col = "black", cex = 0, lwd = 3, verticals = T)
}

for(i in 1:length(OSA.Place)){
  nnc <- nncross(OSA.Place.ppp, OSA.Buy.ppp, k = 1)
  F.func <- ecdf(nnc$dist)
  lines(F.func, col = "blue", cex = 0, lwd = 3, verticals = T)
}

for(i in 1:length(OSA.Place)){
  nnc <- nncross(OSA.Place.ppp, OSA.Vend.ppp, k = 1)
  F.func <- ecdf(nnc$dist)
  lines(F.func, col = "green", cex = 0, lwd = 3, verticals = T)
}

#與台北比較
Windows2 <- as.owin(TPE)

TPE.Metro.df <- data.frame(TPE.Metro@coords)
TPE.Metro.ppp <- ppp(TPE.Metro.df$coords.x1, TPE.Metro.df$coords.x2, Windows2)

TPE.Place.df <- data.frame(TPE.Place@coords)
TPE.Place.ppp <- ppp(TPE.Place.df$coords.x1, TPE.Place.df$coords.x2, Windows2)

TPE.Con.df <- data.frame(TPE.Con@coords)
TPE.Con.ppp <- ppp(TPE.Con.df$coords.x1, TPE.Con.df$coords.x2, Windows2)

TPE.Buy.df <- data.frame(TPE.Buy@coords)
TPE.Buy.ppp <- ppp(TPE.Buy.df$coords.x1, TPE.Buy.df$coords.x2, Windows2)
## Warning: data contain duplicated points
plot(c(0, 1500), c(0, 1), type = 'n', xaxs = "i", yaxs = "i", xlab = "距離", ylab = "F")

for(i in 1:length(TPE.Place)){
  nnc <- nncross(TPE.Place.ppp, TPE.Metro.ppp, k = 1)
  F.func <- ecdf(nnc$dist)
  lines(F.func, col = "red", cex = 0, lwd = 3, verticals = T)
}

for(i in 1:length(TPE.Place)){
  nnc <- nncross(TPE.Place.ppp, TPE.Con.ppp, k = 1)
  F.func <- ecdf(nnc$dist)
  lines(F.func, col = "black", cex = 0, lwd = 3, verticals = T)
}

for(i in 1:length(TPE.Place)){
  nnc <- nncross(TPE.Place.ppp, TPE.Buy.ppp, k = 1)
  F.func <- ecdf(nnc$dist)
  lines(F.func, col = "blue", cex = 0, lwd = 3, verticals = T)
}