library(ggmap)
## Loading required package: ggplot2
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggmap':
## 
##     wind
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(highcharter)
## Highcharts (www.highcharts.com) is a Highsoft software product which is
## not free for commercial and Governmental use
library(treemap)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.3-4, (SVN revision 766)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: D:/Program Files/R/R-3.5.1/library/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: D:/Program Files/R/R-3.5.1/library/rgdal/proj
##  Linking to sp version: 1.3-1
tw<-readOGR("C:/Users/dr476825/Desktop/file/106_2_spatial_analysis/twcounty","TOWN_MOI_1070330",encoding="utf8")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\dr476825\Desktop\file\106_2_spatial_analysis\twcounty", layer: "TOWN_MOI_1070330"
## with 368 features
## It has 7 fields
dp<-tw[tw$COUNTYNAME=="臺北市"|tw$COUNTYNAME=="新北市",]
mrtcount<-read.csv("C:/Users/dr476825/Desktop/intern_proj/mrtcount.csv",encoding="utf8")
dp@data<-left_join(dp@data,mrtcount)
## Joining, by = "TOWNNAME"
dpf<-fortify(dp,region="TOWNNAME")
dpf<-merge(dpf,dp@data,by.x="id",by.y="TOWNNAME")
tpemrtcount<-ggplot(dpf,aes(long,lat,group=group,fill=COUNT,color="white"))+
  geom_polygon()+coord_equal()+
  scale_fill_continuous(high="red",low="white")+
    theme(plot.title = element_text(size=14, face="bold"),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
  theme(axis.title.y=element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())+
  ggtitle("雙北市各行政區捷運站數量圖")+
  guides(color=FALSE)
tpemrtcount

首先以一些基本的視覺畫作為overview

首先這個面量圖是以雙北市行政區作為單位

計算並呈現各行政區內擁有的捷運站數量

mrt<-read.csv("C:/Users/dr476825/Desktop/taipeimertro/mrt_pop_holi.csv")
gg<-ggplot(mrt,aes(TOWN,fill=routec))+
geom_bar(position="stack")+
coord_flip() +
scale_fill_manual(name="路線別",values=c("#f8b61c","#c48c31","#888888","#008659","#0070bd","#e3002c"))+
scale_y_continuous(breaks=seq(0,10,2))+
theme(panel.background = element_blank())+
labs(x = "雙北市行政區",y="捷運站數量")+
labs(title = "雙北市各行政區捷運站數量(並以路線區別)")+
  theme(plot.title = element_text(size=14, face="bold"))
ggplotly(gg)

將下來以長條圖表示

並將路線區別出來

map <- get_map(location =c(121.5155386,25.0457962), zoom = 12)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=25.045796,121.515539&zoom=12&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
mrtsam<-read.csv("C:/Users/dr476825/Desktop/taipeimertro/mrtsample.csv")
head(mrtsam)
##       date time      in.      out count
## 1 2018/3/1    0 松山機場 松山機場     0
## 2 2018/3/1    0 松山機場 中山國中     0
## 3 2018/3/1    0 松山機場 南京復興     0
## 4 2018/3/1    0 松山機場 忠孝復興     0
## 5 2018/3/1    0 松山機場     大安     0
## 6 2018/3/1    0 松山機場 科技大樓     0

這是一開始的原始資料式例

如第二筆資料所代表的意思是

2018/3/1，0時從松山機場進站，中山國中出站的人次數量為0

所以每一小時會有108*107/2+108(同站進同站出)=5886筆資料

就3月有31天而言總共會產生約438萬筆資料

接下來進入資料處理的部分

有了這種形式的原始資料

使我感興趣的是每一站在每一個小時之內總共有多少人潮經過？（出站和入站的人次的總和）

還有人潮的流動方向與情形？（出站和入站人次的差異）

而這個問題將上方的資料每一小時每一站的資料合併就可以被解答（該小時內所有從中山國中站出發的搭乘紀錄都會是被判斷為自

中山國中入站的人次，亦然地，所有抵達中山國中站的搭乘紀錄都會是被判斷為在中山國中出站的人次）

主要的迴圈進行結果會產生台北捷運每個站點每小時平均的進出站人次總和（具有各站人潮量的代表性）

還有各站出站人次減去入站人次之於進出站人次總和的比例（具有各站人潮移動方向的代表性）

因此在下方我就分別針對資料的這兩個面相進行初步的視覺化

colors <- c("#f8b61c", "#c48c31", "#888888","#008659","#0070bd","#e3002c")
routechinese<-aggregate(. ~routec , data=mrt, sum)[,1]
routesum<-aggregate(. ~routec , data=mrt, sum)[,27]
piee<-data.frame(routechinese,routesum)
p <- plot_ly(piee, labels =routechinese, values =routesum, type = 'pie',
        textposition = 'inside',
        textinfo = 'label+percent',
        insidetextfont = list(color = '#FFFFFF'),
        hoverinfo = 'text',
        text = ~paste('人次:', routesum),
        marker = list(colors = colors,
                      line = list(color = '#FFFFFF', width = 1)),
        showlegend = FALSE) %>%
  layout(title = '台北捷運各路線假日一日平均運量',
         xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
         yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE))
p

以pie chart 表現假日的台北捷運各路線一日平均運量

mrt$value<-as.numeric(as.character(mrt$sum))
tm<-treemap(mrt,
index=c("route","X"),
vSize = "value",
vColor = "palettee",
type="color",
title="假日台北捷運各站人潮量（單位：人次）")

hctreemap(tm)
## Warning: 'hctreemap' is deprecated.
## Use 'hctreemap2' instead.
## See help("Deprecated")

以treemap表現假日的台北捷運各站一日平均運量

以路線分別排序（並將兩個不同路線的交會站獨立表示）

並以不同的面積顯現出各站的人潮數量大小之差異(較pie chart更能比較同路線站間之差異)

hours <- rep(0:23,times=2)
inorout<-c(rep("out",24),rep("in",24))
out <- as.numeric(mrt[which(mrt$X=="台北車站"),c(3:26)])*((1+as.numeric(mrt[which(mrt$X=="台北車站"),c(30:53)]))/2)
inm <- as.numeric(mrt[which(mrt$X=="台北車站"),c(3:26)])*((1-as.numeric(mrt[which(mrt$X=="台北車站"),c(30:53)]))/2)
population<-c(out,inm)
tpems <- data.frame(hours,population,inorout)
gg <- ggplot(tpems, aes(x=hours, y=population))
gg <- gg + geom_area(aes(colour=, fill=inorout))+
  labs(colour = "進出站情形")+
  labs(x = "時間(時)")+
  labs(y = "人次數量")+
  labs(title = "台北車站2018年3月假日每日出入站平均人次")+
  scale_fill_manual(values=c("#5599FF", "#FF3333")) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA))
ggplotly(gg)

在這邊先以台北車站作為範例

列出其在3月假日逐時的人次進站與出站的分別情形(也可看出總人潮量的變化趨勢)

attach(mrt)
mappy2<-ggmap(map, darken = c(0.55, "black")) + 
  geom_point(aes(x = lat, y = lon,
                 size =V8,
                 color=V8f,
                 alpha=0.8), 
             data = mrt)+
  scale_color_gradient2(breaks=c(-0.4,0,0.4),
                        label=c("-0.4(乘客移動方向以入站為主)","0","0.4(乘客移動方向以出站為主)"),
                        name="出站與入站人口佔該站人潮數量之差(出-入)",
                        low="blue",high="red",mid="white")+
  scale_size(name="人潮數量(人次)",
             range = c(0, 20),
             breaks=c(1000,3000,5000))+
  scale_alpha_continuous(guide=FALSE)+
  theme(plot.title = element_text(size=14, face="bold"),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())+
  ggtitle("台北捷運假日7時各站人潮數量與出入站方向")
mappy2

接下來以假日7時人潮量與移動方向做為範例

以視覺化表現出各站在7時的人潮數量與移動方向

圓圈則代表台北捷運108個站點

這張圖主要有兩個面向去解釋

一是圓圈大小代表該時該站的人潮量（進站與出站人次的總和）

而另外這邊以out-in（出站人次－入站人次佔該站該小時人潮量的比例）做為方向的指標

圓圈的顏色越偏藍紫代表這個指標越低

說明這個時入該站的移動方向是入站大於出站

另外顏色越偏紅代表這個指標越高

說明這個時間該站的移動方向是出站大於入站

另外白色的話代表此段時間入站跟出站佔總出入站的人次比例幾乎沒有差異

res<-c()
tmp<-c()
td<-c(1:18)
tmpd<-c()
for (j in c(7:24)){
tmpd[j]<-get(paste("V",j,"f",sep=""))[75]
}
tmpd<-tmpd[c(7:24)]
treg<-lm(tmpd~c(1:18))
samtpms<-data.frame(td,tmpd)
reg1<-ggplot(samtpms,aes(td, tmpd))+
  geom_point()+
  geom_abline(linetype="dashed",intercept = treg$coefficients[1], slope = treg$coefficients[2])+
  labs(x = "時間(時),(1代表早上7時)")+
  labs(y = "進出站趨勢(正代表該時出站人口大於入站人口)")+
  labs(title = "海山捷運站2018年3月逐時之進出站趨勢圖")+
  geom_text(aes(x =3 , 
                y = treg$coefficients[1]-0.01,
                label = paste("slope:",
                              round(treg$coefficients[2],3),
                              "p-value:<0.01")))
ggplotly(reg1)

完成了視覺化的呈現

最後一個步驟是以上面已經擁有的各站每小時流動資料進行分析

如果一地主要是住宅區的話

直覺是人群會傾向會是在早上、上午離開，並前往商業、休閒中心

因此早上、上午的移動方向是會偏向入站（代表out-in流動比例應該是負號）

晚上的移動方向是會偏向出站（代表out-in應該是正號）

這邊就可以藉由各站隨時間out-in流動比例的變化進行簡單迴歸分析去判斷是否各站的流動方向是不是隨時間顯著變化且有特定的變化模式

於是統計上的做法是將每一站7~24時的out-in流動比例做為應變數，虛擬時間變數做為自變數，進行簡單回歸

如果跑出來的結果自變數的係數顯著異於零

則代表該地out-in的比例隨時間就統計上來說是有顯著變化的

又如果該地的係數是正

代表從早到晚這個out-in的比例有逐漸上升（或者說是由負轉正）的趨勢（最佳配適線也會明顯地隨時間上升）

則在這邊可以初步的推估此地的型態應該是以住宅區為主（上午坐捷運離開前往商業密集區，晚上再做捷運回家）

如上圖（為海山捷運站的資料）

則海山捷運站附近的土地利用方式就會被我歸類成「住宅區」

tmpd<-c()
for (j in c(7:24)){
tmpd[j]<-get(paste("V",j,"f",sep=""))[32]
}
tmpd<-tmpd[c(7:24)]
treg<-lm(tmpd~c(1:18))
samtpms<-data.frame(td,tmpd)
reg2<-ggplot(samtpms,aes(td, tmpd))+
  geom_point()+
  geom_abline(linetype="dashed",intercept = treg$coefficients[1], slope = treg$coefficients[2])+
    labs(x = "時間(時),(1代表早上7時)")+
  labs(y = "進出站趨勢(正代表該時出站人口大於入站人口)")+
  labs(title = "台北車站捷運站2018年3月逐時之進出站趨勢圖")+
  geom_text(aes(x =3 ,
                y = treg$coefficients[1]-0.01,
                label = paste("slope:",
                              round(treg$coefficients[2],3),
                              "p-value:<0.01")))
ggplotly(reg2)

另外反之如果是out-in的比例隨時間有逐漸下降（或者說是由正轉負）的趨勢

這代表此處附近的土地利用應該是以商業中心或是休閒遊憩區為主

原因是直觀上來說商業中心在早上和上午的進出站方向會傾向於出站

而下午、晚上會傾向於入站（民眾享受完娛樂機能準備搭捷運回家）

如上圖（為台北車站捷運站之資料）

tmpd<-c()
for (j in c(7:24)){
tmpd[j]<-get(paste("V",j,"f",sep=""))[46]
}
tmpd<-tmpd[c(7:24)]
treg<-lm(tmpd~c(1:18))
summary(treg)
## 
## Call:
## lm(formula = tmpd ~ c(1:18))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14412 -0.09965  0.00772  0.06945  0.23907 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.024240   0.056196   0.431    0.672
## c(1:18)     0.002930   0.005192   0.564    0.580
## 
## Residual standard error: 0.1143 on 16 degrees of freedom
## Multiple R-squared:  0.01951,    Adjusted R-squared:  -0.04177 
## F-statistic: 0.3184 on 1 and 16 DF,  p-value: 0.5804
samtpms<-data.frame(td,tmpd)
reg2<-ggplot(samtpms,aes(td, tmpd))+ 
  geom_point()+
  geom_abline(linetype="dashed",intercept = treg$coefficients[1], slope = treg$coefficients[2])+
    labs(x = "時間(時),(1代表早上7時)")+
  labs(y = "進出站趨勢(正代表該時出站人口大於入站人口)")+
  labs(title = "西湖捷運站2018年3月逐時之進出站趨勢圖")+
  theme(legend.position="none")+
  geom_text(aes(color="red",
                x =3 ,
                y = treg$coefficients[1]-0.01,
                label = paste("slope:",
                              round(treg$coefficients[2],3),
                              "p-value:",
                              round(summary(treg)$coefficients[2,4],3))))
ggplotly(reg2)

如果一開始的自變數（時間虛擬變數）係數就沒有顯著異於零

如上圖（西湖捷運站之資料）

代表該地我們沒有辦法用這個方法判斷該站附近主要的土地利用方式

該地可能屬於交通樞紐（大眾通常會在該地轉乘其他交痛工具，因此進出捷運站並無隨時間產生明顯規律）

或甚至該地周遭之夜間之娛樂機能較為豐富（因為捷運的營業時間受限，我們觀察不到1~5時人潮的移動方向）

在這邊可能的情況很多

而避免單純主觀的認定是否趨勢線夠斜而具有向上、向下的趨勢

在這邊藉由迴歸分析中斜率的p統計值判斷

在線性迴歸中斜率的p統計值代表根據模型認為此樣本代表的真實資料的斜率不為0之「犯錯機率」

特別將樣本與真實分開的原因是因為統計都是在用我們觀察到的樣本資料去反推所謂的真實

換句話說如果p統計值越低(通常以0.05作為一判斷門檻)

代表我們越有理由相信真實資料的斜率是異於0的(代表隨著時間會有向上或向下趨勢)

而在這邊我採用的門檻一樣是以0.05為準

而西湖捷運站的斜率p值為0.58

代表這筆資料並不能支撐我們上面的假說

在這個指標中西湖捷運站並沒有單純線性的向上或向下趨勢

for (i in c(1:nrow(mrt)))
  {for (j in c(7:24))
    {tmp[j]<-get(paste("V",j,"f",sep=""))[i]}
  reg<-lm(tmp[c(7:24)]~c(1:18))
  if (summary(reg)$coefficients[2,4]>0.05)
    {res[i]<-"交通樞紐/無法判斷"} 
  else 
    {if (summary(reg)$coefficients[2,1]<0)
      {res[i]<-"商業中心/觀光休閒區"} else{res[i]<-"住宅區"}
    }
}
mrt$res<-res
cbind(as.character(X),res)
tpemetromap2<-ggmap(map, darken = c(0.7, "black")) + 
  geom_point(aes(x = lat, y = lon, color=routec,shape=res,guide=FALSE),size=2.5, data = mrt)+
  scale_colour_manual(values = c("#f8b61c", "#c48c31", "#FFFFFF","#008659","#0070bd","#e3002c"))+
  theme(plot.title = element_text(size=14, face="bold"),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())+
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())+ggtitle("台北捷運各站周圍土地利用之預測結果")
ggplotly(tpemetromap2)

而以迴圈判斷完108個站點周邊可能的土地利用形式後一樣將結果視覺化成上圖

不同顏色代表不同的捷運路線（一樣地把不同路線的交會站獨立出來）

而不同的形狀代表被分類到的土地利用方式（分別有三種：住宅區、商業中心／觀光休閒區、無法判斷）

而同樣的判斷的標準都是：

斜率並不具有顯著性（資料沒有單純線性的向上或向下趨勢）→→　「交通運輸／無法判斷」

斜率具有顯著性　→→　斜率之係數若為正，即判斷為「商業／休閒中心」，反之斜率之係數若為負，則判斷為「住宅區」