Final_Report

Import Data

library(spdep)
library(rgdal)
library(GISTools)
library(spatstat)
library(aspace)
library(ggplot2)
library(ggforce)
library(sf)
library(sp)

setwd("D:/OneDrive - g.ntu.edu.tw/107-2 Cartography and GIS/Final_report")

tw_poly <- readOGR(dsn = ".", layer = "TWN_no_island", encoding="utf8", verbose=FALSE)
tw_poly <- spTransform(tw_poly, CRS("+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +units=m +no_defs"))
school <- readOGR(dsn = ".", layer = "school(height)", encoding="utf8", verbose=FALSE)
school <- spTransform(school, CRS("+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +units=m +no_defs"))
library <- readOGR(dsn = ".", layer = "library(height)", encoding="utf8", verbose=FALSE)
library <- spTransform(library, CRS("+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +units=m +no_defs"))
school@data$N30M_DEM <- as.numeric(school@data$N30M_DEM) #height
library@data$RASTERVALU <- as.numeric(library@data$RASTERVALU) #height

# plot
plot(tw_poly)
points(school, col="red", pch=19)
points(library, col="green", pch=19)

1. 計算各學校與圖書館間的距離，作為圖書資源的不可近性分數

創建2D與3D的情況下，學校與圖書館之間的距離

# Create ppp
school.ppp <- ppp(school@coords[,1],school@coords[,2], tw_poly@bbox[1,],tw_poly@bbox[2,])
library.ppp <- ppp(library@coords[,1],library@coords[,2], tw_poly@bbox[1,],tw_poly@bbox[2,])
school.pp3 <- pp3(school@coords[,1],school@coords[,2], school@data[,21], tw_poly@bbox[1,],tw_poly@bbox[2,],c(0, 3893))
library.pp3 <- pp3(library@coords[,1],library@coords[,2], library@data[,11], tw_poly@bbox[1,],tw_poly@bbox[2,],c(0, 3893))

nnd_2d <- nncross(school.ppp, library.ppp, k=1);head(nnd_2d)

##         dist which
## 1  244.05596     1
## 2   33.12713     1
## 3 2273.30932     1
## 4 3600.47267     1
## 5  558.55051     3
## 6 2505.86822     3

nnd_3d <- nncross(school.pp3, library.pp3, k=1);head(nnd_3d)

##         dist which
## 1  245.65690     1
## 2   35.96396     1
## 3 2273.35243     1
## 4 3601.19250     1
## 5  558.62302     3
## 6 2506.53935     3

由於高度上的距離較難克服，我們希望利用tangent值估算它的坡度，並且以此來對3D的直線距離作加權

判斷在2D與3D的最近鄰居是否有所不同

all(nnd_2d$which == nnd_3d$which)

## [1] TRUE

既然在2D與3D的最近鄰居一樣，我們可以直接利用他們的水平距離和高度差推算坡度(Tangent)，並以坡度加權他們的3D距離

school$dist_2d <- nnd_2d$dist #水平距離
school$dif_height <- abs(school$N30M_DEM[nnd_2d$which] - school$N30M_DEM) #高度差
school$tangent <- school$dif_height / school$dist_2d
school$dist_3d <- nnd_3d$dist * (1+school$tangent)

得到的值即作為圖書資源的不可近性分數，越大代表越需要圖書資源

school$dist_3d[1:30]

##  [1]  245.65690   51.16285 2287.35269 3645.20130  583.62626 2580.55917
##  [7] 3978.43947  475.08908  867.58573 2862.03568 1612.54449  493.28912
## [13]  493.28912 5499.45986 4620.82235  898.77353 1233.69300 1961.57708
## [19] 4956.68591 5774.94088 7419.08249 6507.75457  922.25362 3961.34825
## [25] 1831.38750 3630.80176 3165.03531 2994.48890 4987.08062 4338.25653

接著我們剔除掉離圖書館範圍在4km以內的點，作為需要行動書車資源的學校

school_far <- school[school$dist_3d > 4000,]
plot(tw_poly)
points(school_far, col="red", pch=19)

利用F-function觀察所有需要行動書車資源的學校，發現這些學校有群聚的現象，因此我們希望透過將學校分群，作為可以共享行動書車資源的單位

school_far.ppp <- ppp(school@coords[,1],school@coords[,2], as.owin(tw_poly))

## Warning: data contain duplicated points

CI <- envelope(school_far.ppp, Fest, nsim=99, nrank=5)

## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.

plot(CI)

以K-medoid方法將學校進行分群

由於行動書車的服務範圍皆固定在一縣市內，我們將學校以縣市分別做k-medoid方法以距離分群

首先我們為了判斷分群的結果，以4km作為閾值，若分群後同一群內的各點離最接近中心的點距離皆在閾值內，則此分群結果才會被接受

require(factoextra)
require(cluster)

judge <- function(p, kmedoid){
  p$group <- kmedoid$clustering # 輸入分群後的結果作為分組
  dist <- rep(0, length(kmedoid$id.med)) # 分了幾群就產生幾個dist，來記錄群內距離
  for (i in 1:length(kmedoid$id.med)){ # 每一群內的點做迴圈
    point <- p[p$group == i,]
    if (length(point) == 1){
      dist[i] <- 0 # 若此群只有一個點，則記下距離為0
    }else if(length(point) == 2){
      dist[i] <- dist(method = "euclidean", matrix(point@coords, nrow=2, ncol=2)) #若有兩個點則記錄下他們之間的距離
    }else{
      school.CF <- CF(id=1, points=data.frame(point@coords[,1:2]))
      dist[i] <- max(distances(centre.xy=c(school.CF$CF.x, school.CF$CF.y), matrix(point@coords, nrow=length(point), ncol=2)))
      #若有三個以上的點，則先找出他們的central feature，並且記錄下離central feature最遠的點與其的距離
    }
  }
  if(max(dist) > 4000){ # 判斷是否最大的距離都在4km以內
    return (FALSE)
  }else{
    return (TRUE)
  }
}

接著，對於群內的學校，我們以各個群內的最靠近中心的位置，來決定在這個群內，行動書車的巡迴點應該要設在哪間學校

bookspot <- function(p){
  type <- rep('', length(p))
  for (i in 1:length(p)){ # 每一群內的點做迴圈
    point <- p[i,]
    g <- p[p$group == point$group,] # 找出跟他同群的點
    if (length(g) == 1){
      type[i] <- 'spot' # 若此群只有一個點，則該點即作為書車設置點
    }else if(length(g) == 2){
      type[i] <- ifelse(all(point$dist_3d >= g$dist_3d), 'spot', 'odinary') # 若該點在此群內的圖書資源不可近性最大，則作為書車設置點
    }else{
      school.CF <- CF(id=1, points=data.frame(g@coords[,1:2]))
      spot <- c(school.CF$CF.x, school.CF$CF.y) # 距離最近的就是書車設置點
      type[i] <- ifelse(all(point@coords == spot), 'spot', 'odinary') # 若該點書車設置點則記錄下來
    }
  }
  return(type)
}

最後，我們來看在每一群內可以服務到多少學生，作為行動書車巡迴時的考量優先程度判斷

amount <- function(p){
  num <- rep('', length(p))
  for (i in 1:length(p)){ # 每一群內的點做迴圈
    point <- p[i,]
    g <- p[p$group == point$group,] # 找出跟他同群的點
    num[i] <- as.numeric(sum(as.numeric(g$N___10)) + sum(as.numeric(g$N___11)))
  }
  return(num)
}

開始每縣市進行分群

## 新北市

school_far$N___ <- as.character(school_far$N___)
new_tp <- subset(school_far, startsWith(c(school_far$N___), '[01]'))

for (i in 1: 35){
  new_tp.kmedoid.cluster <- pam(new_tp@coords, k=i)
  tof <- judge(new_tp, new_tp.kmedoid.cluster)
  if (tof) break
}
new_tp$group <- as.character(new_tp.kmedoid.cluster$clustering)
new_tp$type <- bookspot(new_tp)
new_tp$amount <- amount(new_tp)

新北市分群結果:13群

fviz_cluster(new_tp.kmedoid.cluster,  # 分群結果
             data = new_tp@coords,    # 資料
             geom = c("point"))       # 點 (point)

new_tp.poly <- fortify(tw_poly[10,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = new_tp.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=new_tp@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=new_tp@data[new_tp@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=new_tp@data[new_tp@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 宜蘭縣
yilan <- subset(school_far, startsWith(c(school_far$N___), '[02]'))

for (i in 1:35){
  if (i < length(yilan)){
    yilan.kmedoid.cluster <- pam(yilan@coords, k=i)
    tof <- judge(yilan, yilan.kmedoid.cluster)
    if (tof) break
  }else{
    print("Can't do the grouping.")
    break
  }
}
yilan$group <- as.character(seq(1, length(yilan)), 1)
yilan$type <- bookspot(yilan)
yilan$amount <- amount(yilan)

宜蘭縣分群結果:9群

yilan.poly <- fortify(tw_poly[1,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = yilan.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=yilan@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=yilan@data[yilan@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=yilan@data[yilan@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 桃園市
taoyuan <- subset(school_far, startsWith(c(school_far$N___), '[03]'))

for (i in 1: 35){
  taoyuan.kmedoid.cluster <- pam(taoyuan@coords, k=i)
  tof <- judge(taoyuan, taoyuan.kmedoid.cluster)
  if (tof) break
}
taoyuan$group <- as.character(taoyuan.kmedoid.cluster$clustering)
taoyuan$type <- bookspot(taoyuan)
taoyuan$amount <- amount(taoyuan)

桃園市分群結果:8群

fviz_cluster(taoyuan.kmedoid.cluster,  # 分群結果
             data = taoyuan@coords,    # 資料
             geom = c("point"))       # 點 (point)

taoyuan.poly <- fortify(tw_poly[6,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = taoyuan.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=taoyuan@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=taoyuan@data[taoyuan@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=taoyuan@data[taoyuan@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 新竹縣
hsinchu <- subset(school_far, startsWith(c(school_far$N___), '[04]'))

for (i in 1: 35){
  hsinchu.kmedoid.cluster <- pam(hsinchu@coords, k=i)
  tof <- judge(hsinchu, hsinchu.kmedoid.cluster)
  if (tof) break
}
hsinchu$group <- as.character(hsinchu.kmedoid.cluster$clustering)
hsinchu$type <- bookspot(hsinchu)
hsinchu$amount <- amount(hsinchu)

新竹縣分群結果:12群

fviz_cluster(hsinchu.kmedoid.cluster,  # 分群結果
             data = hsinchu@coords,    # 資料
             geom = c("point"))       # 點 (point)

hsinchu.poly <- fortify(tw_poly[12,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = hsinchu.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=hsinchu@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=hsinchu@data[hsinchu@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=hsinchu@data[hsinchu@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 苗栗縣
miaoli <- subset(school_far, startsWith(c(school_far$N___), '[05]'))

for (i in 1: 35){
  miaoli.kmedoid.cluster <- pam(miaoli@coords, k=i)
  tof <- judge(miaoli, miaoli.kmedoid.cluster)
  if (tof) break
}
miaoli$group <- as.character(miaoli.kmedoid.cluster$clustering)
miaoli$type <- bookspot(miaoli)
miaoli$amount <- amount(miaoli)

苗栗縣分群結果:15群

fviz_cluster(miaoli.kmedoid.cluster,  # 分群結果
             data = miaoli@coords,    # 資料
             geom = c("point"))       # 點 (point)

miaoli.poly <- fortify(tw_poly[5,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = miaoli.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=miaoli@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=miaoli@data[miaoli@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=miaoli@data[miaoli@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 台中市
taichung <- subset(school_far, startsWith(c(school_far$N___), '[06]'))

for (i in 1: 35){
  taichung.kmedoid.cluster <- pam(taichung@coords, k=i)
  tof <- judge(taichung, taichung.kmedoid.cluster)
  if (tof) break
}
taichung$group <- as.character(taichung.kmedoid.cluster$clustering)
taichung$type <- bookspot(taichung)
taichung$amount <- amount(taichung)

台中市分群結果:16群

fviz_cluster(taichung.kmedoid.cluster,  # 分群結果
             data = taichung@coords,    # 資料
             geom = c("point"))       # 點 (point)

taichung.poly <- fortify(tw_poly[16,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = taichung.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=taichung@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=taichung@data[taichung@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=taichung@data[taichung@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 彰化縣
changhua <- subset(school_far, startsWith(c(school_far$N___), '[07]'))

for (i in 1: 35){
  changhua.kmedoid.cluster <- pam(changhua@coords, k=i)
  tof <- judge(changhua, changhua.kmedoid.cluster)
  if (tof) break
}
changhua$group <- as.character(changhua.kmedoid.cluster$clustering)
changhua$type <- bookspot(changhua)
changhua$amount <- amount(changhua)

彰化縣分群結果:9群

fviz_cluster(changhua.kmedoid.cluster,  # 分群結果
             data = changhua@coords,    # 資料
             geom = c("point"))       # 點 (point)

changhua.poly <- fortify(tw_poly[15,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = changhua.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=changhua@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=changhua@data[changhua@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=changhua@data[changhua@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 南投縣
nantou <- subset(school_far, startsWith(c(school_far$N___), '[08]'))

for (i in 30:35){
  nantou.kmedoid.cluster <- pam(nantou@coords, k=i)
  tof <- judge(nantou, nantou.kmedoid.cluster)
  if (tof) break
}
nantou$group <- as.character(nantou.kmedoid.cluster$clustering)
nantou$type <- bookspot(nantou)
nantou$amount <- amount(nantou)

南投縣分群結果:33群

fviz_cluster(nantou.kmedoid.cluster,  # 分群結果
             data = nantou@coords,    # 資料
             geom = c("point"))       # 點 (point)

nantou.poly <- fortify(tw_poly[3,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = nantou.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=nantou@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=nantou@data[nantou@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=nantou@data[nantou@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 雲林縣
yunlin <- subset(school_far, startsWith(c(school_far$N___), '[09]'))

for (i in 10:35){
  yunlin.kmedoid.cluster <- pam(yunlin@coords, k=i)
  tof <- judge(yunlin, yunlin.kmedoid.cluster)
  if (tof) break
}
yunlin$group <- as.character(yunlin.kmedoid.cluster$clustering)
yunlin$type <- bookspot(yunlin)
yunlin$amount <- amount(yunlin)

雲林縣分群結果:16群

fviz_cluster(yunlin.kmedoid.cluster,  # 分群結果
             data = yunlin@coords,    # 資料
             geom = c("point"))       # 點 (point)

yunlin.poly <- fortify(tw_poly[9,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = yunlin.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=yunlin@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=yunlin@data[yunlin@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=yunlin@data[yunlin@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 嘉義縣
chiayi <- subset(school_far, startsWith(c(school_far$N___), '[10]'))

for (i in 20:35){
  chiayi.kmedoid.cluster <- pam(chiayi@coords, k=i)
  tof <- judge(chiayi, chiayi.kmedoid.cluster)
  if (tof) break
}
chiayi$group <- as.character(chiayi.kmedoid.cluster$clustering)
chiayi$type <- bookspot(chiayi)
chiayi$amount <- amount(chiayi)

嘉義縣分群結果:26群

fviz_cluster(chiayi.kmedoid.cluster,  # 分群結果
             data = chiayi@coords,    # 資料
             geom = c("point"))       # 點 (point)

chiayi.poly <- fortify(tw_poly[14,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = chiayi.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=chiayi@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=chiayi@data[chiayi@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=chiayi@data[chiayi@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 台南市
tainan1 <- subset(school_far, startsWith(c(school_far$N___), '[11]'))
tainan2 <- subset(school_far, startsWith(c(school_far$N___), '[21]'))
tainan <- rbind(tainan1,tainan2)

for (i in 10:35){
  tainan.kmedoid.cluster <- pam(tainan@coords, k=i)
  tof <- judge(tainan, tainan.kmedoid.cluster)
  if (tof) break
}
tainan$group <- as.character(tainan.kmedoid.cluster$clustering)
tainan$type <- bookspot(tainan)
tainan$amount <- amount(tainan)

台南市分群結果:16群

fviz_cluster(tainan.kmedoid.cluster,  # 分群結果
             data = tainan@coords,    # 資料
             geom = c("point"))       # 點 (point)

tainan.poly <- fortify(tw_poly[19,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = tainan.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=tainan@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=tainan@data[tainan@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=tainan@data[tainan@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 高雄市
kaohsiung <- subset(school_far, startsWith(c(school_far$N___), '[12]'))

for (i in 1:35){
  kaohsiung.kmedoid.cluster <- pam(kaohsiung@coords, k=i)
  tof <- judge(kaohsiung, kaohsiung.kmedoid.cluster)
  if (tof) break
}
kaohsiung$group <- as.character(kaohsiung.kmedoid.cluster$clustering)
kaohsiung$type <- bookspot(kaohsiung)
kaohsiung$amount <- amount(kaohsiung)

高雄市分群結果:15群

fviz_cluster(kaohsiung.kmedoid.cluster,  # 分群結果
             data = kaohsiung@coords,    # 資料
             geom = c("point"))       # 點 (point)

kaohsiung.poly <- fortify(tw_poly[7,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = kaohsiung.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=kaohsiung@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=kaohsiung@data[kaohsiung@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=kaohsiung@data[kaohsiung@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 屏東縣
pintung <- subset(school_far, startsWith(c(school_far$N___), '[13]'))

for (i in 15:35){
  pintung.kmedoid.cluster <- pam(pintung@coords, k=i)
  tof <- judge(pintung, pintung.kmedoid.cluster)
  if (tof) break
}
pintung$group <- as.character(pintung.kmedoid.cluster$clustering)
pintung$type <- bookspot(pintung)
pintung$amount <- amount(pintung)

屏東縣分群結果:17群

fviz_cluster(pintung.kmedoid.cluster,  # 分群結果
             data = pintung@coords,    # 資料
             geom = c("point"))       # 點 (point)

pintung.poly <- fortify(tw_poly[4,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = pintung.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=pintung@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=pintung@data[pintung@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=pintung@data[pintung@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 台東縣
taitung <- subset(school_far, startsWith(c(school_far$N___), '[14]'))

for (i in 15:35){
  taitung.kmedoid.cluster <- pam(taitung@coords, k=i)
  tof <- judge(taitung, taitung.kmedoid.cluster)
  if (tof) break
}
taitung$group <- as.character(taitung.kmedoid.cluster$clustering)
taitung$type <- bookspot(taitung)
taitung$amount <- amount(taitung)

台東縣分群結果:19群

fviz_cluster(taitung.kmedoid.cluster,  # 分群結果
             data = taitung@coords,    # 資料
             geom = c("point"))       # 點 (point)

taitung.poly <- fortify(tw_poly[18,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = taitung.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=taitung@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=taitung@data[taitung@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=taitung@data[taitung@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

## 花蓮縣
hualien <- subset(school_far, startsWith(c(school_far$N___), '[15]'))

for (i in 15:35){
  hualien.kmedoid.cluster <- pam(hualien@coords, k=i)
  tof <- judge(hualien, hualien.kmedoid.cluster)
  if (tof) break
}
hualien$group <- as.character(hualien.kmedoid.cluster$clustering)
hualien$type <- bookspot(hualien)
hualien$amount <- amount(hualien)

花蓮縣分群結果:19群

fviz_cluster(hualien.kmedoid.cluster,  # 分群結果
             data = hualien@coords,    # 資料
             geom = c("point"))       # 點 (point)

hualien.poly <- fortify(tw_poly[2,])

## Regions defined for each Polygons

ggplot() + geom_polygon(data = hualien.poly, aes(x = long, y = lat, group = group), fill="gray") +
  geom_point(data=hualien@data, aes(x=X, y=Y, color=group, shape=type), size=3) + coord_fixed(1.0) +
  geom_circle(data=hualien@data[hualien@data$type=='spot',], aes(x0=X, y0=Y, color=group, fill=group, r=4000), alpha=0.3) +
  geom_text(data=hualien@data[hualien@data$type=='spot',], aes(x=X, y=Y, label=amount), hjust=1, vjust=1, size=3)

Final_Report

Cartus You

2019年5月31日

1. 計算各學校與圖書館間的距離，作為圖書資源的不可近性分數

判斷在2D與3D的最近鄰居是否有所不同

既然在2D與3D的最近鄰居一樣，我們可以直接利用他們的水平距離和高度差推算坡度(Tangent)，並以坡度加權他們的3D距離

得到的值即作為圖書資源的不可近性分數，越大代表越需要圖書資源

接著我們剔除掉離圖書館範圍在4km以內的點，作為需要行動書車資源的學校

利用F-function觀察所有需要行動書車資源的學校，發現這些學校有群聚的現象，因此我們希望透過將學校分群，作為可以共享行動書車資源的單位

以K-medoid方法將學校進行分群

首先我們為了判斷分群的結果，以4km作為閾值，若分群後同一群內的各點離最接近中心的點距離皆在閾值內，則此分群結果才會被接受

接著，對於群內的學校，我們以各個群內的最靠近中心的位置，來決定在這個群內，行動書車的巡迴點應該要設在哪間學校

最後，我們來看在每一群內可以服務到多少學生，作為行動書車巡迴時的考量優先程度判斷

開始每縣市進行分群

新北市分群結果:13群

宜蘭縣分群結果:9群

桃園市分群結果:8群

新竹縣分群結果:12群

苗栗縣分群結果:15群

台中市分群結果:16群

彰化縣分群結果:9群

南投縣分群結果:33群

雲林縣分群結果:16群

嘉義縣分群結果:26群

台南市分群結果:16群

高雄市分群結果:15群

屏東縣分群結果:17群

台東縣分群結果:19群

花蓮縣分群結果:19群