1071QG Mid1

計量地理學及實習　期中考一參考答案

助教杜承軒 2018.10.28

實作主題：建立高齡友善社會—老人安全與救護可近性評估
本次試題將評估都市老人安全的救護資源供給與需求。以高雄原市區與鳳山區作為研究區，各里的老年人口為目標族群，消防隊可作為提供救護服務的單位，利用上述的 2SFCA 方法，找出建構都市老人安全網的救護資源缺乏地區，可作為補強老人安全網的政策建議。以下將 2SFCA 方法分成四個主要的程序，依序實作以下題目。

※圖資：
＊KH_vill.shp：高雄原市區與鳳山區村里面資料（TWD97 TM2）
　 TOWN (TOWN_ID)：行政區（行政區編號）
　 VILLAGE (V_ID)：村里（村里編號）
　 A65UP_CNT：老年人口數（單位：人）
＊RescueCorps.csv：消防隊的點位資料（經緯度）
　 NAME (ENG_NAME)：消防隊的單位名稱（英文名稱）
　 LON,LAT：消防隊的經緯度

※座標參考系統 CRS 之 proj4 格式：
　[EPSG:4326] WGS84 經緯度：
　+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
　[EPSG:3826] TWD97 TM2：
　+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000
　+y_0=0 +ellps=GRS80 +units=m +no_defs

＊初步環境建置與讀取檔案
rm(list=ls())
library(GISTools);library(rgdal);library(sp);library(dplyr);library(magrittr)
setwd("D:/1071QG/Mid1")
KH=readOGR(dsn = ".", layer = "KH_vill", encoding="unicode",verbose = F)
RC.data=read.csv("RescueCorps.csv")
一、 [30%] 繪製各行政區的救護資源空間分布

[5%] 將消防隊的位點位轉換成 TWD97 TM2 的空間點資料。
#資料來源是經緯度，因此用經緯度座標來讀取，再轉換成TWD97TM2。
RC=SpatialPointsDataFrame(cbind(RC.data$LON,RC.data$LAT),RC.data,proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84"))
RC=spTransform(RC,KH@proj4string)
[5%] 將村里圖層合併成行政區圖層（透過 TOWN 或 TOWN_ID）。
KH.town=gUnaryUnion(KH,KH$TOWN)
[5%] 計算各行政區老年人口數，並將資料加進行政區圖層中。
pop=xtabs(KH$A65UP_CNT~KH$TOWN)
pop=data.frame(pop=as.vector(pop), row.names =names(pop))
KH.town=SpatialPolygonsDataFrame(KH.town,pop)
[5%] 計算行政區有多少救援單位點，並將行政區內消防隊的個數除以老年人口數（單位：萬人），定義為行政區內人均救護資源量分數。
KH.town$counts=poly.counts(RC,KH.town)
KH.town$score=KH.town$counts/KH.town$pop*10000
[5%] 呈現各行政區人均救護資源量分數繪製成面量圖，以分位數法(quantile)分成五類。

[5%] 加上地圖要素：圖名、圖例、指北針、比例尺。
par(mar=c(0,0,1,0))
quant5=auto.shading(KH.town$score,cutter=quantileCuts,n=5,col=brewer.pal(5,"GnBu"))
choropleth(KH.town,KH.town$score,quant5)
title("高雄原市區與鳳山區人均救護資源量")
choro.legend(187000,2515322,quant5,between = "~",fmt="%.2f",cex=0.8,title='人均救護資源分數')
map.scale(170000,2495983,8000, "公里",4,2) 
north.arrow(168508,2513620,800,col= '#AAAAAA') 

二、 [40%] 計算各消防隊能夠提供的服務量

[25%] 以消防隊之環域 2.5 公里為搜尋區，利用村里圖層涵蓋範圍的面積比例，來計算每個消防隊涵蓋的老年人口數。

RC2.5=gBuffer(RC,width = 2500,byid = T)
KH.RC=gIntersection(KH,RC2.5,byid=T)
KH.RC_name=strsplit(names(KH.RC), " ") 
KH.RC$KH_id=unlist(lapply(KH.RC_name,function(x) x[1]))%>%as.numeric
KH.RC$RC_id=unlist(lapply(KH.RC_name,function(x) x[2]))%>%as.numeric

KH.RC$area=as.vector(poly.areas(KH.RC))
KH.RC$KH.area=poly.areas(KH)[KH.RC$KH_id+1]
KH.RC$KH.pop=KH$A65UP_CNT[KH.RC$KH_id+1]

KH.RC$pop=round(KH.RC$area/KH.RC$KH.area*KH.RC$KH.pop)
RC$pop=xtabs(KH.RC$pop~KH.RC$RC_id)

[5%] 將消防隊區分成大隊／中隊／分隊，分別設定供給量為 10／7／5。
（若電腦無法辨識中文字元，可使用 ENG_NAME 欄位 Brigade/Squadron/Squad 來判斷）

RC$supply[substring(RC$NAME,3,3)=="大"]=10
RC$supply[substring(RC$NAME,3,3)=="中"]=7
RC$supply[substring(RC$NAME,3,3)=="分"]=5

[5%] 計算自身的供給量($S_j$)除以搜尋區內老年人口數(單位：萬人)($\sum P_k$)，為提供服務量的分數($R_j$)。

RC$Rj=RC$supply/RC$pop%>%as.numeric*10000

[5%] 以表格呈現消防隊的單位名稱與提供服務量的分數。

data.frame(NAME=RC$NAME,ENG=RC$ENG_NAME,Rj=RC$Rj)

##        NAME                ENG         Rj
## 1  十全分隊     Shi Quan Squad  0.7855583
## 2  大昌中隊  Da Chang Squadron  1.2826856
## 3  大林分隊       Da Lin Squad 17.6304654
## 4  小港分隊    Xiao Gang Squad  4.3836577
## 5  中華中隊 Zhong Hua Squadron  1.5675385
## 6  五甲分隊       Wu Jia Squad  1.2253100
## 7  右昌分隊    You Chang Squad  3.0767337
## 8  左營分隊     Zuo Ying Squad  2.3261224
## 9  左營分隊     Zuo Ying Squad  1.8774407
## 10 成功分隊   Cheng Gong Squad  1.8856539
## 11 前金分隊     Qian Jin Squad  0.9454298
## 12 前鎮中隊 Qian Zhen Squadron  2.5748547
## 13 苓雅大隊    Ling Ya Brigade  1.3923698
## 14 高桂中隊   Gao Gui Squadron  7.1001116
## 15 新莊分隊   Xin Zhuang Squad  1.8902877
## 16 新興中隊  Xin Xing Squadron  0.9322271
## 17 楠梓中隊    Nan Zi Squadron  8.9732086
## 18 瑞隆分隊     Rui Long Squad  0.9780142
## 19 鼎金大隊   Ding Jin Brigade  2.8836726
## 20 鼓山分隊      Gu Shan Squad  2.2431584
## 21 旗津分隊       Qi Jin Squad 10.0200401
## 22 鳳山中隊 Feng Shan Squadron  2.1497451
## 23 鳳祥大隊 Feng Xiang Brigade  7.4928818

三、 [35%] 計算每個里的地理可近性分數

[20%] 建立村里中心的圖層，以村里中心點找出距離 2.5 公里內的消防隊，並加總提供服務量的分數($\sum R_j$)，作為該里的地理可近性分數($A_i$)。
KH.pt=gCentroid(KH,byid = T)
KH.RC_dist=gWithinDistance(RC,KH.pt,dist = 2500,byid = T)
sumRj=function(x) sum(RC$Rj[x])
KH$Ai=apply(KH.RC_dist,1,sumRj)
[5%] 以直方圖(histogram)呈現各里地理可近性分數的分布， x 軸為地理可近性分數，y 軸為次數頻率。
par(mar=c(4.2,4.2,2,1))
hist(KH$Ai,xlab="可近性分數",ylab="次數頻率",main="可近性分數直方圖",breaks=30,col="lightblue")
[10%] 比較鳳山、前鎮、左營等三個行政區，各里地理可近性分數的平均值，是否有顯著差異？請列出虛無假設，再進行檢定得出結論。
#H0=三個區平均無差異
#Ha=三個區平均有差異
#ANOVA檢定
KH2=subset(KH@data,KH$TOWN %in% c("鳳山區","前鎮區","左營區"))
aov(KH2$Ai~KH2$TOWN)%>%summary
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## KH2$TOWN      2  246.4  123.20   13.05 5.33e-06 ***
## Residuals   171 1614.7    9.44                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#p-value<0.001
#顯著水準0.01時，拒絕H0，三個區平均有差異

四、 [20%] 可近性地圖：可近性分數 vs. 老年人口數的雙變數面量圖

[5%] 繪製各村里可近性分數與老年人口數的散布圖。
oldM=mean(KH$A65UP_CNT)
accQ1=quantile(KH$Ai,0.25)
plot(KH$A65UP_CNT~KH$Ai,pch=21,bg="#FFD70088",
    xlab="可近性分數",ylab="老年人口數",main="各村里可近性分數與老年人口數散布圖")
abline(h=oldM,v=accQ1)
[5%] 可近性分數以第一四分位數來區分低與高兩類，老年人口數以平均值區分低與高兩類。
KH$col=1
KH$col[KH$A65UP_CNT<oldM]=KH$col[KH$A65UP_CNT<oldM]+1
KH$col[KH$Ai>=accQ1]=KH$col[KH$Ai>=accQ1]+2
[10%] 請參考右圖的顏色來繪製雙變數面量圖。
紅色：可近性分數低而老年人口數高，代表此地老人安全的風險很高
黃色：可近性分數低但老年人口數低的村里
藍色：老年人口數高但可近性方便，表示救護的供給充裕
綠色：代表老年人口數不多且救護資源無虞之處。
col=c("red","yellow","blue","green")
plot(KH,col=col[KH$col])
title("可近性與老年人口數之雙變數面量圖")
legend(184000,2517500,c('資源缺乏的風險區','供給不足但需求較小','需求多但供給充裕','救護資源無虞'),fill=col,x.intersp = 0.5,text.width =15000)
map.scale(170000,2493000,8000, "公里",4,2) 
north.arrow(168508,2513620,800,col= '#AAAAAA') 

計量地理學及實習 期中考一參考答案

助教 杜承軒 2018.10.28

計量地理學及實習　期中考一參考答案

助教杜承軒 2018.10.28