import JSON files
# Load the package required to read JSON files.
#library("rjson")
# Give the input file name to the function.
#result <- fromJSON(file = "mrt.geojson")
#new=unlist(result)
Creating mrt csv file by using ‘result’
#new1=c()
#for (i in 1:3280){
#new1[i]=new[i]
#}
#exit=c()
#et=seq(7,3280,9)
#for (e in et){
#exit[e]=new1[e]
#}
#exit=na.omit(exit)
#x=c()
#xx=seq(12,3280,9)
#for (k in xx){
#x[k]=new1[k]}
#x=na.omit(x)
#y=c()
#yy=seq(13,3280,9)
#for (l in yy){
#y[l]=new1[l]
#}
#y=na.omit(y)
#mrt=data.frame(x,y)
#write.table(mrt,file="/Users/eyo/Desktop/mrt.csv",sep=",",row.names=F)
library(spdep)
library(rgdal)
library(GISTools)
library(splancs)
library(raster)
library(ggtern)
library(dplyr)
library(aspace)
library(sp)
library(ggplot2)
library(spatstat)
setwd("/Users/eyo/Desktop")
toshin <- readOGR(dsn = "SA_Mid2", layer = "toshin", encoding="big5") # For WGS84 coordinate
tp <- readOGR(dsn = "Lab1", layer = "Taipei_Vill", encoding="big5") # Taipei Village shp
TWPOP <- readOGR(dsn = "Data2", layer = "Popn_TWN2", encoding="big5") # Taiwan shp
pop <- readOGR(dsn = "pop_points", layer = "pop_points", encoding="big5") #population points shapefile
gym <- readOGR(dsn = "gym_shapefile", layer = "gym", encoding="big5") # Gym Shapefile
bus <- readOGR(dsn = "bus", layer = "busstop97", encoding="big5") # Bus shapefile
mrt = read.table("mrt.csv" , sep = "," , header = T) # MRT shapefile
ubike = read.table("ubike.csv" , sep = "," , header = T) #Ubike shapefile
newmrt <- SpatialPoints(mrt, proj4string = toshin@proj4string) #Converting mrt coordinates into spatial data
newbike <- SpatialPoints(ubike, proj4string = toshin@proj4string) #Converting ubike coordinates into spatial data
newmrt=spTransform(newmrt,pop@proj4string) # From WGS84 to TWD97
newubike=spTransform(newbike,pop@proj4string)
### Creating spatialpointsdataframe
data.mrt <- SpatialPointsDataFrame(newmrt, mrt,
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 "))
data.ubike <- SpatialPointsDataFrame(newbike, ubike,
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 "))
Creating shapefile to share with teammates
#writeOGR(data.mrt, dsn = "/Users/eyo/Desktop/mrt", layer = "mrt.shp",driver = "ESRI Shapefile")
#writeOGR(data.ubike, dsn = "/Users/eyo/Desktop/ubike", layer = "mrt.shp",driver = "ESRI Shapefile")
Creating a shapefile that includes all transportation systems (MRT,bus,Ubike)
pbus=data.frame(bus@coords)
pubike=data.frame(newubike@coords)
pmrt=data.frame(newmrt@coords)
x=c(pbus$coords.x1,pubike$lng,pmrt$LON)
y=c(pbus$coords.x2,pubike$lat,pmrt$LAT)
xy=data.frame(x,y)
newxy <- SpatialPoints(xy, proj4string = pop@proj4string)
data.xy <- SpatialPointsDataFrame(newxy, xy,
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 "))
#writeOGR(data.xy, dsn = "/Users/eyo/Desktop/xy", layer = "xy.shp",driver = "ESRI Shapefile")
Plotting Population choropleth map
Popn<-tp@data$CENSUS
Popn<-as.numeric(paste(Popn))
par(mar = c(0,5,3,3))
shades = auto.shading(Popn, n = 5, cols = brewer.pal(5, "Greens"))
choropleth(tp, Popn, shades)
choro.legend(317198.9-2500,2789179.6-10000,shades)
north.arrow( 292688,2768239
,500,col= 'black')
title('population choropleth map of Taipei')
Plotting KDE map by using gym and population points data.
#ggtern
TPE_lim=c(tp@bbox[1,], tp@bbox[2,])
tpp <- gUnaryUnion(tp,id=tp$TOWN)
pts1=as.points(coordinates(pop)[,1],coordinates(pop)[,2])
pts2=as.points(coordinates(gym)[,1],coordinates(gym)[,2])
kde.pts1 = kde2d(pts1[,1], pts1[,2], 2000, 100,TPE_lim)
kde.pts2 = kde2d(pts2[,1], pts2[,2], 2000, 100,TPE_lim)
kde.diff = kde.pts1
kde.diff$z = kde.pts1$z-kde.pts2$z
image(kde.diff, asp=1,col=brewer.pal(9,'RdYlBu'))
masker=poly.outer(as.points(tp@bbox[1,],tp@bbox[2,]),tp)
## Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
## unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4
## strings
add.masking(masker, col='white')
plot(tpp,add=T)
Creating points boxes with actual Taipei City shape
pop_x.coor<-pop@coords[,1]
pop_y.coor<-pop@coords[,2]
gym_x.coor<-gym@coords[,1]
gym_y.coor<-gym@coords[,2]
t_x.coor<-data.xy@coords[,1]
t_y.coor<-data.xy@coords[,2]
TN_BND=gBuffer(tp,width=500)
TN_BND.Coord<- TN_BND@polygons[[1]]@Polygons[[1]]@coords
x1<-rev(TN_BND.Coord[,1])
y1<-rev(TN_BND.Coord[,2])
xx<-cbind(x1, y1)
# building WINDOW file
PTS_bnd <- owin(poly=xx)
pop.ppp<-ppp(pop_x.coor,pop_y.coor, window = PTS_bnd)
gym.ppp<-ppp(gym_x.coor,gym_y.coor, window = PTS_bnd)
## Warning: data contain duplicated points
t.ppp<-ppp(t_x.coor,t_y.coor, window = PTS_bnd)
## Warning: 1825 points were rejected as lying outside the specified window
## Warning: data contain duplicated points
Bivariate F function between population and transportation (order from 1 to 5)
nnd=nncross(pop.ppp,t.ppp)
f=ecdf(nnd[,1])
plot(f,cex=0,verticals=T,xaxs="i",yaxs="i",col="red",main="The five most adjacent public transport stations from the population",xlab="distance(m)",ylab="F",lwd=8,xlim=c(0,700))
p=c(8,6,4,2,1)
color=c("red","orange","yellow","green","blue")
for (i in 1:5){
nnd=nncross(pop.ppp,t.ppp,k=i)
ff=ecdf(nnd[,1])
lines(ff,cex=0,verticals=T,xaxs="i",yaxs="i",col=color[i],lwd=p[i])
}
Bivariate F function between gym and transportation (order from 1 to 5)
nnd=nncross(gym.ppp,t.ppp)
f=ecdf(nnd[,1])
plot(f,cex=0,verticals=T,xaxs="i",yaxs="i",col="red",main="The five most adjacent public transport stations from the gym",xlab="distance(m)",ylab="F",lwd=8)
p=c(8,6,4,2,1)
color=c("red","orange","yellow","green","blue")
for (i in 1:5){
nnd=nncross(gym.ppp,t.ppp,k=i)
ff=ecdf(nnd[,1])
lines(ff,cex=0,verticals=T,xaxs="i",yaxs="i",col=color[i],lwd=p[i])
}
Analysis through spatial autocorrelation
gym.count=poly.counts(gym,tp)
for(i in 1:length(tp)){
tp$gym.desity[i]=gym.count[i]/poly.areas(tp)[i]
}
coords=coordinates(tp)
tp.nb =dnearneigh(coords,d1=0,d2=5000)
#TW_li.nb = poly2nb(TW.li)
tp.nb.in = include.self(tp.nb)
tp.nb.w.in = nb2listw (tp.nb.in)
Gi = localG(tp$gym.desity,tp.nb.w.in)
# Moran Spatial Correlogram
cor=sp.correlogram(tp.nb, tp$gym.desity, order=5, method="I", style="W",zero.policy=T)
print(cor); plot(cor,main="Moran Spatial Correlogram",xlab="Density classes",ylab="Moran's I statistic")
## Spatial correlogram for tp$gym.desity
## method: Moran's I
## estimate expectation variance standard deviate
## 1 (456) 6.7395e-02 -2.1978e-03 2.4030e-05 14.1969
## 2 (456) -8.0022e-02 -2.1978e-03 1.2728e-05 -21.8144
## 3 (456) -4.2839e-02 -2.1978e-03 4.6837e-05 -5.9385
## 4 (370) -2.5066e-02 -2.7100e-03 7.2768e-04 -0.8288
## 5 (75) 3.1601e-02 -1.3514e-02 5.8263e-03 0.5910
## Pr(I) two sided
## 1 (456) < 2.2e-16 ***
## 2 (456) < 2.2e-16 ***
## 3 (456) 2.877e-09 ***
## 4 (370) 0.4072
## 5 (75) 0.5545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
LG =as.vector(Gi)
quad = c()
quad[1.96000<LG] = 1 # cluster
quad[1.64500<LG & LG<=1.960000] = 2 # cluster
quad[-1.65<LG & LG<=1.645000] = 3 # non cluster
quad[-1.96<LG & LG<=-1.645000] = 4 # non cluster
quad[LG<=-1.960000] = 5 # non cluster
lm.palette=colorRampPalette(c("blue","green","yellow"), space = "rgb")
colors=lm.palette(5)
plot(tp, border="grey", col= colors[quad], main = "Standardized Gi * valuesTaipei of Gym Village Gym Density")
legend("bottomright", c("1.960001 - 3.000000","1.644999 - 1.960000","-1.644999 - 1.645000","-1.959999 - -1.645000","-3.038610 - -1.960000"), fill=colors, bty ="n", cex=0.7,y.intersp=1,x.intersp=1)
