109-2 空間分析 期中考二
助教 杜承軒 2021.04.19
討論一區域內供給與需求的空間關係時，常衡量可近性指標，來評估資源的分配問題。當一地可近性指標越低，代表獲得服務的中間阻隔越大，因此能偵測出資源缺乏的地區。
過去常用流動搜尋法 (FCA; floating catchment area method) 的概念計算可近性指標，如 Peng (1997) 計算每個行政區中心，固定距離內有多少供給量及需求量，並將兩者之和相除，得到可近性指標。
Luo & Wang (2003) 進一步改良評估方法，提出兩階段流動搜尋法 (2SFCA; two-step floating catchment area method) 來衡量可近性指標，兩階段步驟如下：
第一階段對於每一個供給點\(j\)，以\(d_0\)為半徑的畫出搜尋區，計算區域內的人口(\(P_k\))，並用自身的供給量(\(S_j\))，除以區域內的人口，作為能提供服務量的分數(\(R_j\))。
\(R_j=\frac{S_j}{\sum_{k∈\{d_{kj}≤d_0\}} { P_k}}\)
第二階段對於每一個需求點\(i\)，以\(d_0\)為半徑的畫出搜尋區，加總區域內每一個供給點\(j\)能提供服務量的分數(\(R_j\))，得到可近性的分數(\(A_i\))。
\(A_i=\sum_{j∈\{d_{ij}≤d_0\}}R_j\)
第一部分
從健康促進的觀點，小學附近若有較多的速食店，容易取得高熱量食物可能影響學生的飲食健康。以台北市小學與速食店的分布作為研究區對象，來理解兩者之間的空間關係。
＊初步環境建置與讀取檔案
library(sf);library(tmap);library(aspace);library(dplyr);library(units);library(gridExtra)
setwd("D:/1092SA/Mid2/Data")
TP=st_read("TP_District.shp")
fastfood=st_read("TP_Fastfood.shp")
school=st_read("TP_School.shp")
windowsFonts(TOP = windowsFont("Topedia Sans TW Beta"))
- 以學生人數作加權，分別找出公立小學與私立小學之加權平均中心點，計算兩點距離相差幾公尺（回答整數）？[10%]
school.mean=function(type){
pts=school[school$type==type,]
pts_mean=mean_centre(points=st_coordinates(pts),weighted=T,weights=pts$student)
pts_mean[,2:3]%>%as.numeric%>%st_point%>%st_sfc%>%st_sf%>%st_set_crs(3826)
}
public.mean=school.mean("public")
private.mean=school.mean("private")
ppdist=st_distance(public.mean,private.mean)sprintf("公立小學與私立小學之學生人數加權平均中心點相差%.f公尺。",ppdist)## [1] "公立小學與私立小學之學生人數加權平均中心點相差1067公尺。"
- 市政府決定找一間小學召開學生健康的會議，由各校校長出席。假設校長都由各自的小學出發，欲使直線距離總和最小化，應該要選擇哪一間小學來當作開會地點（回答小學名稱或編號）？[10%]
school.xy=school%>%st_coordinates%>%data.frame
school.center=CF(points=school.xy)
meet.id=which(school.xy$X==school.center$CF.x&school.xy$Y==school.center$CF.y)sprintf("編號：%d，名稱：%s。",school$SID[meet.id],school$name[meet.id])## [1] "編號：29，名稱：中正國小。"#或
school.median=median_centre(points=school.xy)
school.median=school.median[,2:3]%>%as.numeric%>%st_point%>%st_sfc%>%st_sf%>%st_set_crs(3826)
meet=school[which.min(st_distance(school.median,school)),]sprintf("編號：%d，名稱：%s。",meet$SID,meet$name)## [1] "編號：29，名稱：中正國小。"
- 以行政區為底圖，繪製速食店的平均中心點及1.5倍標準圓。[10%]
fast.sdd=calc_sdd(points = st_coordinates(fastfood))
fast.mean.sf=c(sddatt$CENTRE.x, sddatt$CENTRE.y)%>%st_point%>%st_sfc%>%st_sf%>%st_set_crs(3826)
fast.rad=fast.sdd$SDD.radius
fast.SDD.sf=st_buffer(fast.mean.sf,fast.rad*1.5)
qtm(TP)+tm_shape(fast.mean.sf)+tm_symbols(col="red")+tm_shape(fast.SDD.sf)+tm_borders("red")+tm_layout(frame =F)
- 以行政區為底圖，繪製每個行政區內小學的平均中心點與標準橢圓。[15%]
map=qtm(TP)
for(i in 1:nrow(TP)){
dis.school=school[TP[i,],]
calc_sde(points = st_coordinates(dis.school))
center=c(sdeatt$CENTRE.x, sdeatt$CENTRE.y)%>%st_point%>%st_sfc%>%st_sf%>%st_set_crs(3826)
SDE=data.frame(x= sdeloc$x, y= sdeloc$y)%>%st_as_sf(coords=c("x","y"))%>%
st_set_crs(3826)%>%st_combine%>%st_cast("POLYGON")%>%st_sf
map=map+tm_shape(center)+tm_symbols(col="red",size=.3)+tm_shape(SDE)+tm_borders("red")
}
map+tm_layout(frame =F)
- 利用FCA的概念，以1公里環域作為搜尋區，利用小學搜尋區內速食點座位數除以搜尋區內學生數，計算每間學校的可近性分數。篩選出可近性分數大於0.2的小學，標示為易達速食點的學校。以行政區為底圖，將這些學校標記在地圖上。[15%]
dist.SF=st_is_within_distance(school,fastfood,1000,sparse=F)
dist.SS=st_is_within_distance(school,school,1000,sparse=F)
school$FCA=(dist.SF%*%fastfood$seat)/(dist.SS%*%school$student)
qtm(TP)+tm_shape(school[school$FCA>0.2,])+tm_symbols(col="red")+tm_layout(frame =F)
SID name FCA
66 66 福星國小 1.4340278
63 63 西門國小 0.4169203
19 19 大安國小 0.3792487
45 45 忠孝國小 0.3141935
64 64 老松國小 0.2793638
48 48 日新國小 0.2310306
13 13 信義國小 0.2254098
17 17 博愛國小 0.2107809
38 38 濱江國小 0.2060653
101 101 文湖國小 0.2060653
85 85 辛亥國小 0.2023547這是錯的
school$wrong=(dist.SF%*%fastfood$seat)/school$student
qtm(TP)+tm_shape(school[school$wrong>0.2,])+tm_symbols(col="red")+tm_layout(frame =F)
***
第二部分
欲探討都市內老人安全救護的空間分布狀況，以高雄原市區與鳳山區作為研究區，老年人口為目標族群，救援單位作為提供服務的單位，找出老人安全的風險區。利用2SFCA的概念，設定2.5 公里範圍為鄰近定義，計算各村里的可近性指標。
＊初步環境建置與讀取檔案
setwd("D:/1092SA/Mid2/Data")
KH=st_read("KH_Village.shp")
Fire=st_read("KH_FireStation.shp")
- 各行政區內消防局的總救援能量，除以該行政區老年人口（萬人），可得「人均救援能量」。繪製主題地圖，以面量圖(choropleth map)表示該數值在各行政區分布（需繪製面量圖圖例）。[15%]
KH_Dist=KH%>%group_by(TOWN)%>%summarise(elder=sum(A65UP_CNT))
KH_Dist.Fire=st_contains(KH_Dist,Fire,sparse=F)
KH_Dist$resource=KH_Dist.Fire%*%Fire$resource%>%as.numeric
KH_Dist$Rpc=KH_Dist$resource/KH_Dist$elder*10^4
tm_shape(KH_Dist)+tm_polygons("Rpc",title="人均救援能量")+tm_layout(frame = F, fontfamily = "TOP")
KH_Dist%>%st_drop_geometry## # A tibble: 12 x 4
## TOWN elder resource Rpc
## * <chr> <dbl> <dbl> <dbl>
## 1 三民區 50568 22 4.35
## 2 小港區 19072 17 8.91
## 3 左營區 24355 15 6.16
## 4 前金區 5976 5 8.37
## 5 前鎮區 30570 17 5.56
## 6 苓雅區 32133 10 3.11
## 7 新興區 10648 7 6.57
## 8 楠梓區 21588 12 5.56
## 9 鼓山區 19658 12 6.10
## 10 旗津區 4421 5 11.3
## 11 鳳山區 46940 22 4.69
## 12 鹽埕區 5209 0 0
- 消防局服務量=消防局救援能量÷該消防局2.5 公里範圍內村里的老年人口總數（萬人）。
計算每個消防局的服務量分數，利用村里圖層涵蓋範圍的面積比例，來計算每個救援單位涵蓋的老年人口數，並依照編號(ID)的順序，用列表列出答案（編號、名稱、服務量；服務量請計算到小數第二位）。[15%]
KH$villArea=st_area(KH)
Fire2500=st_buffer(Fire,2500)
inter=st_intersection(KH,Fire2500)
inter$area=st_area(inter)
inter$elder=round(inter$A65UP_CNT*inter$area/inter$villArea)%>%as.numeric
Fire_POP=inter%>%group_by(ID)%>%summarise(NAME=NAME[1],resource=resource[1],elder=sum(elder))%>%st_drop_geometry
Fire_POP$Rj=round(Fire_POP$resource/Fire_POP$elder*10^4,2)
Fire_POP[,c(1,2,5)]# A tibble: 23 x 3
ID NAME Rj
<dbl> <chr> <dbl>
1 1 高桂中隊 7.03
2 2 苓雅大隊 1.37
3 3 鼓山分隊 2.21
4 4 鼎金大隊 2.84
5 5 瑞隆分隊 0.97
6 6 楠梓中隊 8.89
7 7 鳳祥大隊 7.26
8 8 鳳山中隊 2.13
9 9 旗津分隊 9.71
10 10 小港分隊 4.32
# ... with 13 more rows#或
Fire.POP=inter%>%group_by(ID)%>%summarise(elder=sum(elder))%>%st_drop_geometry
Fire=Fire%>%left_join(Fire.POP)
Fire$Rj=round(Fire$resource/Fire$elder*10^4,2)
Fire.df=Fire%>%st_drop_geometry%>%.[,c(2,1,5)]
Fire.df[order(Fire.df$ID),] ID NAME Rj
14 1 高桂中隊 7.03
13 2 苓雅大隊 1.37
20 3 鼓山分隊 2.21
19 4 鼎金大隊 2.84
18 5 瑞隆分隊 0.97
17 6 楠梓中隊 8.89
23 7 鳳祥大隊 7.26
22 8 鳳山中隊 2.13
21 9 旗津分隊 9.71
4 10 小港分隊 4.32
3 11 大林分隊 17.51
2 12 大昌中隊 1.27
1 13 十全分隊 0.77
8 14 左營小隊 2.29
7 15 右昌分隊 3.06
6 16 五甲分隊 1.21
5 17 中華中隊 1.54
12 18 前鎮中隊 2.53
11 19 前金分隊 0.93
10 20 成功分隊 1.84
9 21 左營分隊 1.86
16 22 新興中隊 0.92
15 23 新莊分隊 1.87
- 可近性分數=該里中心點 2.5 公里範圍內消防局的服務量總和。
建立村里中心的圖層，計算每個里的可近性分數，並畫出散布圖(scatter plot)呈現各里人口與可近性分數的關係（X軸：村里可近性分數｜Y軸：村里老年人口）。[15%]
KHc=st_centroid(KH)
KHc.F=st_is_within_distance(KHc,Fire,2500,F)
KH$Ai=KHc$Ai=KHc.F%*%Fire$Rj%>%as.numeric
ggplot(KH,aes(x=Ai,y=A65UP_CNT))+geom_point()+xlab("可近性分數")+ylab("老年人口")+theme_minimal()+theme(text=element_text(family="TOP"))
- 找出每個行政區中，老年人口密度最高的村里，標示這些村里中心位置，繪製泡泡地圖(bubble map)，調整圓圈的大小來表現這些村里中心的2.5 公里範圍內消防局的數量（需繪製圓圈大小的圖例）。[15%]
KH$den=KH$A65UP_CNT/KH$villArea%>%as.numeric()
TopDenID=function(x){
town=KH[KH$TOWN==x,]
town$V_ID[which.max(town$den)]
}
KH.topdenID=sapply(unique(KH$TOWN),TopDenID)
KH.topden=KH[KH$V_ID%in%KH.topdenID,]%>%st_centroid
KH.topden['2.5公里內消防局數量']=st_is_within_distance(KH.topden,Fire,2500)%>%lengths
qtm(KH)+tm_shape(KH.topden)+tm_bubbles("2.5公里內消防局數量",col="red")+tm_layout(frame = F, fontfamily = "TOP")
- 用合適的統計檢定方法，比較老年人口總數最高的五個行政區內，各里老年人口密度的平均值，是否有顯著差異？需列出假設檢定與結論。[15%]
elder_town_name=order(KH_Dist$elder,decreasing=T)[1:5]%>%KH_Dist$TOWN[.]
KHD5=KH[KH$TOWN %in% elder_town_name,]
A5=aov(den~TOWN,KHD5)%>%summary設定顯著水準為\({\alpha=0.05}\)
\(H_0:\)五個行政區之各里老年人口密度的平均值無顯著差異
\(H_a:\)五個行政區之各里老年人口密度的平均值有顯著差異檢定多母體平均值是否相同，使用AVOVA檢定。
A5 Df Sum Sq Mean Sq F value Pr(>F)
TOWN 4 0.0001063 2.658e-05 5.408 0.000316 ***
Residuals 324 0.0015927 4.916e-06
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\(F=5.408\)，\(p-value=0\)。
\(p-value<\alpha\)，拒絕虛無假設。
結論：老年人口總數最高的五個行政區內之各里老年人口密度的平均值有顯著差異。
- 利用各里人口與可近性分數的平均值，來分類該數值的高群與低群，如右表，繪製各村里可近性分數與老年人口數的雙變數面量圖(bivariate choropleth)。[15%]
oid=KH$A65UP_CNT>=mean(KH$A65UP_CNT)
aid=KH$Ai>=mean(KH$Ai)
typename=c("高需求高供給","資源缺乏","資源無虞","低需求低供給")
col=c("blue","red","green3","yellow2")
KH$type[ oid& aid]=typename[1]
KH$type[ oid&!aid]=typename[2]
KH$type[!oid& aid]=typename[3]
KH$type[!oid&!aid]=typename[4]
KH$type=KH$type%>%ordered(levels=typename)
KH.box=ggplot(KH,aes(x=Ai,y=A65UP_CNT))+geom_point(aes(col=type))+scale_color_manual("",values=col)+xlab("可近性分數")+ylab("老年人口")+theme_minimal()+theme(text=element_text(family="TOP"),legend.position = "bottom")
KH.map=ggplot(KH)+geom_sf(aes(fill=type))+scale_fill_manual("類別",values=col,labels=typename)+theme_void()+theme(text=element_text(family="TOP"))
grid.arrange(KH.box,KH.map, ncol=2)
tmap_mode("view")
tm_shape(KH)+tm_polygons("type",palette=col,title="類別")+tm_layout(frame = F, fontfamily = "TOP")