FINAL REPORT
地理二 何雨潔、地理二 劉雨函
KDE
library(splancs)
library(GISTools)
library(rgdal)
Popn.TWN <- readOGR(dsn = "D:/Class/大二下/空間分析/190304/Data2", layer = "Popn_TWN2", encoding="unicode", verbose = FALSE)
TPE = subset(Popn.TWN, Popn.TWN$COUNTY == "臺北市")
hb_107108 = read.csv(file = "D:/Class/大二下/空間分析/Presentation/Clean Data/hb_107&108.csv")
tpe_sqr_bnd = read.csv(file = "D:/Class/大二下/空間分析/190429/KDE_Lab/tpe_sqr_bnd.csv")
pts_hb_107108 = as.points(hb_107108$Response_X, hb_107108$Response_Y)
bnd = as.points(tpe_sqr_bnd$X, tpe_sqr_bnd$Y)
PTS = SpatialPoints(pts_hb_107108, TPE@proj4string)
KDE_PTS = kde.points(PTS, 1000, 1000, lims = TPE)
plot(KDE_PTS)
masker = poly.outer(KDE_PTS, TPE)
add.masking(masker, col = "white")
plot(TPE, add = T, border = "white" )
Monte Carlo Significance Test (F function)
H0: Pattern is random. Ha: Pattern is not random.
α = 0.1
obs（黑線）大致位於包絡區間右下方，因此拒絕H0，接受Ha－－[Pattern is not random.]
library(spatstat)
bnd<-bbox(TPE)
x.range<-bnd[1,]
y.range<-bnd[2,]
TPE.Windows<-owin(xrange=x.range, yrange=y.range)
hb.ppp<-ppp(pts_hb_107108[,1], pts_hb_107108[,2], TPE.Windows)
random.ppp = rpoint(466, win=TPE.Windows)
nnd = nncross(random.ppp, hb.ppp)
plot(ecdf(nnd$dist), main = "", pch = 16)
rp_1 = c()
for(i in 1:100){
nn1<-rpoint(466, win=TPE.Windows)
nnd1<-nncross(random.ppp, nn1)
#rp_1[i] = mean(nnd1)
F.random = ecdf(nnd1$dist)
lines(F.random,col=colorRampPalette(brewer.pal(8, "Set2"))(8)[1])
}
plot(ecdf(nnd$dist), add = T)
Spatial clustering of establishments
claw = read.csv(file = "D:/Class/大二下/空間分析/Presentation/Clean Data/claw.csv")
claw.ppp<-ppp(claw$Response_X, claw$Response_Y, TPE.Windows)
gas = read.csv(file = "D:/Class/大二下/空間分析/Presentation/Clean Data/gas.csv")
gas.ppp<-ppp(gas$x_coord, gas$y_coord, TPE.Windows)
hospital = read.csv(file = "D:/Class/大二下/空間分析/Presentation/Clean Data/Hospital.csv")
hospital.ppp<-ppp(hospital$Response_X, hospital$Response_Y, TPE.Windows)
hotel = read.csv(file = "D:/Class/大二下/空間分析/Presentation/Clean Data/Hotel.csv")
hotel.ppp<-ppp(hotel$Response_X, hotel$Response_Y, TPE.Windows)
food = read.csv(file = "D:/Class/大二下/空間分析/Presentation/Clean Data/FOOD.csv")
food.ppp<-ppp(food$Response_X, food$Response_Y, TPE.Windows)
retail = read.csv(file = "D:/Class/大二下/空間分析/Presentation/Clean Data/RETAIL.csv")
retail.ppp<-ppp(retail$Response_X, retail$Response_Y, TPE.Windows)
light = read.csv(file = "D:/Class/大二下/空間分析/Presentation/Clean Data/street_light.csv")
light.ppp<-ppp(light$x_coord, light$y_coord, TPE.Windows)
nnd_claw = nncross(claw.ppp, hb.ppp, k = 1)
nnd_gas = nncross( gas.ppp, hb.ppp , k = 1)
nnd_hospital = nncross( hospital.ppp, hb.ppp, k = 1)
nnd_hotel = nncross( hotel.ppp, hb.ppp, k = 1)
nnd_food = nncross( food.ppp, hb.ppp, k = 1)
nnd_retail = nncross( retail.ppp, hb.ppp, k = 1)
nnd_light = nncross(light.ppp, hb.ppp, k = 1)
par(mar = c(2,2,1,1))
plot(ecdf(nnd_hotel$dist), main = "", cex = 0, verticals = T, xaxs = "i", yaxs = "i", col = colorRampPalette(brewer.pal(8, "Set2"))(8)[1])
legend("bottomright",c("飯店","夾娃娃機臺","餐飲業","零售業", "醫院", "加油站", "路燈"),col=colorRampPalette(brewer.pal(8, "Set2"))(8)[1:7],lty=1,lwd=2, cex = 1.2, xjust = 1)
plot(ecdf(nnd_claw$dist), main = "" , add = T, cex = 0, verticals = T, xaxs = "i", yaxs = "i", col = colorRampPalette(brewer.pal(8, "Set2"))(8)[2])
plot(ecdf(nnd_food$dist), main = "", add = T, cex = 0, verticals = T, xaxs = "i", yaxs = "i", col = colorRampPalette(brewer.pal(8, "Set2"))(8)[3])
plot(ecdf(nnd_retail$dist), main = "", add = T, cex = 0, verticals = T, xaxs = "i", yaxs = "i", col = colorRampPalette(brewer.pal(8, "Set2"))(8)[4])
plot(ecdf(nnd_hospital$dist), main = "", add = T, cex = 0, verticals = T, xaxs = "i", yaxs = "i", col = colorRampPalette(brewer.pal(8, "Set2"))(8)[5])
plot(ecdf(nnd_gas$dist), main = "", add = T, cex = 0, verticals = T, xaxs = "i", yaxs = "i", col = colorRampPalette(brewer.pal(8, "Set2"))(8)[6])
plot(ecdf(nnd_light$dist), main = "", add = T, cex = 0, verticals = T, xaxs = "i", yaxs = "i", col = colorRampPalette(brewer.pal(8, "Set2"))(8)[7])
SDE
library(aspace)
library(RColorBrewer)
par(mar = c(1,1,1,1))
plot(TPE, asp=1, axes=FALSE, xlab="", ylab="", pch=16, cex=0.6, main = "")
SDE <- calc_sde(points = pts_hb_107108)
plot_sde(plotnew = F, plotpoints = T , plotcentre = T ,sde.col = colorRampPalette(brewer.pal(8, "Set2"))(8)[2] ,centre.col = colorRampPalette(brewer.pal(8, "Set2"))(8)[2], titletxt = "", points.pch = 16, sde.lwd = 3)
多元迴歸分析
###簡單線性迴歸
#竊盜
pts_hb_sp_107108 = SpatialPoints(pts_hb_107108, proj4string = TPE@proj4string)
cou = poly.counts(pts_hb_sp_107108, TPE)
TPE$HB = cou/sum(cou)
TPE$AREA = poly.areas(TPE)/1000000
#飯店
sp_hotel = SpatialPoints(as.points(hotel$Response_X, hotel$Response_Y) , proj4string = TPE@proj4string)
cou_hotel = poly.counts(sp_hotel , TPE)
TPE$hotel = cou_hotel/TPE$AREA
#夾娃娃機
sp_claw = SpatialPoints(as.points(claw$Response_X, claw$Response_Y) , proj4string = TPE@proj4string)
cou_claw = poly.counts(sp_claw , TPE)
TPE$claw = cou_claw/TPE$AREA
#餐飲業
sp_food = SpatialPoints(as.points(food$Response_X, food$Response_Y) , proj4string = TPE@proj4string)
cou_food = poly.counts(sp_food , TPE)
TPE$food = cou_food/TPE$AREA
#零售業
sp_retail = SpatialPoints(as.points(retail$Response_X, retail$Response_Y) , proj4string = TPE@proj4string)
cou_retail = poly.counts(sp_retail , TPE)
TPE$retail = cou_retail/TPE$AREA
#醫院
sp_hospital = SpatialPoints(as.points(hospital$Response_X, hospital$Response_Y) , proj4string = TPE@proj4string)
cou_hospital = poly.counts(sp_hospital , TPE)
TPE$hospital = cou_hospital/TPE$AREA
#加油站
sp_gas = SpatialPoints(as.points(gas$x_coord, gas$y_coord) , proj4string = TPE@proj4string)
cou_gas = poly.counts(sp_gas , TPE)
TPE$gas = cou_gas/TPE$AREA
#路燈
sp_light = SpatialPoints(as.points(light$x_coord, light$y_coord) , proj4string = TPE@proj4string)
cou_light = poly.counts(sp_light , TPE)
TPE$light = cou_light/TPE$AREA
try_reg=lm(HB~hotel+claw+food+retail+hospital+gas+light,data=TPE)
summary(try_reg)##
## Call:
## lm(formula = HB ~ hotel + claw + food + retail + hospital + gas +
## light, data = TPE)
##
## Residuals:
## 220 221 222 223 224 225
## -0.0023551 -0.0005476 0.0012824 -0.0030594 0.0009840 0.0012527
## 226 227 228 229 230 231
## 0.0029159 -0.0252803 0.0064986 0.0119862 -0.0193330 0.0256555
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1289535 0.0351790 3.666 0.0215 *
## hotel 0.0874335 0.0425321 2.056 0.1090
## claw 0.0218324 0.0087217 2.503 0.0665 .
## food 0.0207839 0.0334080 0.622 0.5676
## retail -0.0167828 0.0118961 -1.411 0.2311
## hospital 0.0669312 0.1273612 0.526 0.6270
## gas 0.0334724 0.0999592 0.335 0.7546
## light -0.0000736 0.0001129 -0.652 0.5501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02171 on 4 degrees of freedom
## Multiple R-squared: 0.7669, Adjusted R-squared: 0.3589
## F-statistic: 1.88 on 7 and 4 DF, p-value: 0.2828###Moran Test (Greater)
library(spdep)
tw.nb = poly2nb(TPE)
TPE$RES = residuals(try_reg)
lw = nb2listw(tw.nb, style = "B", zero.policy = T)
lm.morantest(try_reg, lw)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = HB ~ hotel + claw + food + retail + hospital +
## gas + light, data = TPE)
## weights: lw
##
## Moran I statistic standard deviate = -1.2871, p-value = 0.901
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## -0.239643000 -0.123496623 0.008142964###Moran Test (Two-sided)
lm.morantest(try_reg, lw, alternative = "two.sided")##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = HB ~ hotel + claw + food + retail + hospital +
## gas + light, data = TPE)
## weights: lw
##
## Moran I statistic standard deviate = -1.2871, p-value = 0.1981
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I Expectation Variance
## -0.239643000 -0.123496623 0.008142964