空間分析　 期中考一

助教 杜承軒 2020.04.03

一、概念題 (12點)

現有東京都行政區面資料，以及東京都內國鐵車站點位資料。將兩筆資料讀取至R中，設定前者名稱為Tokyo，資料格式為SpatialPolygonsDataFrame，參考投影座標為EPSG:4326，後者名稱為JR，資料格式為SpatialPointsDataFrame，參考投影座標為EPSG:3095，又已知東京都的範圍約介於139°E~140°E、35.5°~36°N之間。
請分別回答以下各小題的程式，東京面資料、車站點位或底圖（各題紅色行數的資料）是否能成功疊合？若答案為否，請更正程式碼使其能夠成功執行疊圖（更正正確才給分）。

  1. [2點]
proj4string(JR)=proj4string(Tokyo)
plot(Tokyo)
points(JR)
## 否，第1行改成：JR=spTransform(JR,Tokyo@proj4string)

正確寫出spTransform來更改CRS：2點
寫出spTransform但不完整(如第2題)：1點

  1. [2點]
spTransform(Tokyo,JR@proj4string)
plot(Tokyo)
points(JR)
## 否，第1行改成：Tokyo=spTransform(Tokyo,JR@proj4string)

正確寫出spTransform來更改CRS：2點

  1. [2點]
TK.f=fortify(spTransform(Tokyo,CRS('+init=epsg:3095')))
ggplot()+geom_polygon(data=TK.f,aes(x=long,y=lat,group=group))+
          geom_point(aes(x=JR@coords[,1],y=JR@coords[,2]))
## 是

是：2點／其他：0點

  1. [3點]
map=openmap(c(139,36),c(140,35.5),12,'esri-topo')
plot(map)
points(spTransform(JR,CRS('+init=epsg:4326')))
## 否
## 第1行改成：map=openmap(c(36,139),c(35.5,140),12,'esri-topo')→經緯度相反
## 第3行改成：points(spTransform(JR,osm())
## 或，第3行改成：points(spTransform(JR,CRS('+init=epsg:3857'))

有更改經緯度：+1點
有更改後面的座標：+2點

  1. [3點]
new.Tokyo=spTransform(Tokyo,CRS('+init=epsg:4326'))
new.JR=spTransform(JR,CRS('+init=epsg:4326'))
min=c(new.Tokyo@bbox[2,1],new.Tokyo@bbox[1,1])
max=c(new.Tokyo@bbox[2,2],new.Tokyo@bbox[1,2])
map=openmap(min,max,13,'osm')
plot(map)
plot(new.Tokyo,add=T)
plot(new.JR,add=T)
## 否
## 第5~6行之間加入：map=openproj(map,CRS(+init=epsg:4326))
## 或，第7行改成：plot(spTransform(new.Tokyo,osm()),add=T)
## 　　第8行改成：plot(spTransform(new.JR,osm()),add=T)

更改前五行內容：+0點
更改後面的座標：+3點

二、實作題 (38點)

＊初步環境建置與讀取檔案

library(ggplot2);library(GISTools);library(rgdal);library(dplyr)
windowsFonts(JH = windowsFont("微軟正黑體"))
setwd("D:/1082SA/Mid1/Data/")
SW=readOGR("Taiwan_SW.shp",encoding = "unicode")
PH=readOGR("pharmacy.shp")
sample=read.csv("sample.csv")
  1. ggplot繪製面量圖

定義鄉鎮藥局服務量（間／萬人）＝該鄉鎮藥局數量（間）÷該鄉鎮總人口（萬人）。
目標繪製出臺南市各鄉鎮市區的服務量面量圖，以及臺南市藥局分布點位。

(1) [5點]計算出臺南市各鄉鎮市區服務量，並以平均─標準差分類法分成五類。

轉換成一致的座標系統後，計算、整理各鄉鎮的人口、面積、藥局數、服務量

SW=spTransform(SW,PH@proj4string)
SW$pop=SW$A0A14_CNT+SW$A15A64_CNT+SW$A65UP_CNT
SW$area=poly.areas(SW)/10^6
SW$pharmacy=poly.counts(PH,SW)
SW$service=SW$pharmacy/SW$pop*10^4

計算出鄉鎮服務量(臺南市即可)：2點

以「平均─標準差法」分組

TN=SW[SW$COUNTY=='臺南市',]
TN.shade=auto.shading(TN$service,4,sdCuts,cols = brewer.pal(5,"RdBu"))
TN$class=cut(TN$service,c(-Inf,TN.shade$breaks,Inf),labels = 1:5)
#或是
m=mean(TN$service);sd=sd(TN$service)
TN$class=cut(TN$service,c(-Inf,m-1.5*sd,m-0.5*sd,m+0.5*sd,m+1.5*sd,Inf),labels = 1:5)
xtabs(~TN$class)
## TN$class
##  1  2  3  4  5 
##  2 11 12  9  3

以「平均─標準差法」正確地分組：3點

(2) [9點]透過ggplot2套件，繪製臺南市服務量面量圖，以及臺南市藥局分布點位。

篩選出台南市的藥局點位

PH.TN=PH[colSums(gWithin(PH,TN,byid=T))==1,]
#或是
#PH.TN=gIntersection(PH,TN)

篩選出台南市藥局點位：2點

一般方法繪圖

par(mar=c(0,1,1,1),family ="JH")
choropleth(TN,TN$service,TN.shade, main = "臺南市藥局服務量與藥局分布地圖")
map.scale(152946,2537068,20000, "公里",2,10) #比例尺
north.arrow(152416,2583206,1000,lab='N',col='black') #指北針
choro.legend(205459,2550674,TN.shade,under="<",over=">",between="~",fmt = "%.2f",cex=0.8,title="鄉鎮藥局服務量") #圖例
points(PH.TN,pch=20)


一般方法繪圖：最高4點
缺少面資料或點資料：各扣2點
面量圖有誤：扣1點
缺少圖名、圖例、比例尺、指北針或有誤：各扣1點(最多2點)
(分組方式有誤：扣在前面資料處理的部分)
(沒有先篩選出台南市的資料：扣在前面資料處理的部分)

ggplot繪圖

#以經緯度網格來繪圖需先轉換CRS
TN=spTransform(TN,CRS("+init=epsg:4326"))
PH.TN=spTransform(PH.TN,CRS("+init=epsg:4326"))
#進行fortify，再將資料對回後繪圖
TN.f=fortify(TN,region="TOWN")
TN.f=left_join(TN.f,TN@data,by=c("id"="TOWN"))
ggplot()+geom_polygon(data=TN.f,aes(x=long,y=lat,group=group,fill=class),col='black')+coord_equal()+
  scale_fill_manual("鄉鎮藥局服務量",values=brewer.pal(5,"RdBu"),labels=c("<-1.5sd","-1.5~-0.5sd","-0.5~0.5sd","0.5~1.5sd",">1.5sd"))+
  geom_point(aes(x=PH.TN@coords[,1],y=PH.TN@coords[,2]),size=0.8)+
  theme_minimal()+labs(title ="臺南市藥局服務量與藥局分布地圖",x="經度",y="緯度")+
  theme(legend.position = "right",plot.title = element_text(size=16,hjust = 0.5),text=element_text(family="JH"))


ggplot繪圖：最高7點
缺少面資料或點資料：各扣2點
面量圖有誤(fortify等部分有錯誤)：扣2點
缺少圖名、圖例、經緯度網格或有誤：各扣1點(最多2點)
(分組方式有誤：扣在前面資料處理的部分)
(沒有先篩選出台南市的資料：扣在前面資料處理的部分)

  1. [12點]計算服務受限人數

假設藥局的服務範圍為1 km，且服務範圍不受行政區邊界限制，並預設在鄉鎮內的人口是均勻分布的。請利用鄉鎮市區服務範圍的面積比例，計算出「嘉義縣、嘉義市無法被藥局所服務到的人數」分別為多少人？

「藥局一公里環域」與「嘉義縣市」面資料，進行截切，計算各鄉鎮服務範圍面積

CY=SW[SW$COUNTY%in%c("嘉義縣","嘉義市"),]  #將嘉義縣市的鄉鎮取出
PH1K=gBuffer(PH,width=1000) #藥局一公里環域，融合環域結果即可，因此byid=F
CY1K=gIntersection(CY,PH1K,byid=T) #計算各鄉鎮x各藥局服務範圍
CY.id=unlist(lapply(strsplit(names(CY1K)," "),function(x) x[1])) #取出鄉鎮id
CY1K.area=poly.areas(CY1K)/10^6 #計算各鄉鎮x各藥局服務範圍面積

篩選出嘉義縣市：1點
藥局一公里環域(注意CRS)：2點
透過gIntersection，正確計算各鄉鎮服務範圍面積，與ID等等：4點

服務範圍面積對回CY資料

#將服務範圍面積對回CY資料
CY$area1K=0 
for(i in 1:length(CY.id))  CY[CY.id[i],'area1K']=CY1K.area[i]
#或是用left_join將CY資料對回截切資料
CY.data=data.frame(CY.id,CY1K.area)
CY$CY.id=rownames(CY@data)
CY.data=left_join(CY.data,CY@data)
#注意：這個步驟會有一個鄉鎮資料(大埔鄉)消失，只能算影響人口，直接算未被影響人口會有誤(參考下面的方法)

將Intersect資料，正確的與原始鄉鎮面積、人口對齊：3點

計算出各鄉鎮不在服務範圍的人口，再加總至縣市

CY$nosevrice=round((1-CY$area1K/CY$area)*CY$pop)
CY$COUNTY=as.character(CY$COUNTY)
data.frame(xtabs(nosevrice~COUNTY,CY))
##   COUNTY   Freq
## 1 嘉義市 119964
## 2 嘉義縣 440919
#或是
CY.data$popsevrice=round((CY.data$area1K/CY.data$area)*CY.data$pop)
total=xtabs(pop~COUNTY,SW@data)
popsevrice=xtabs(popsevrice~COUNTY,CY.data)
nosevrice=data.frame(total-popsevrice)
nosevrice[nosevrice$COUNTY%in%c("嘉義市","嘉義縣"),]
##   COUNTY   Freq
## 2 嘉義市 119964
## 3 嘉義縣 440919

計算出各鄉鎮不在服務範圍的人口：1點
加總至縣市：1點
※答案可能因為四捨五入等方式的選擇，有些微差距是正常的。
【整大題的方法、順序可能很多，以這個作答的比例為原則，詳細點數可以自己拿捏。】

  1. [12點]評估藥局可近性差異

隨機抽樣藥局來調查空間中供需分配關係，結果如sample.csv。
定義藥局可近性分數＝藥品供給量÷（服務範圍內人口數÷服務範圍內藥局總數量），其中，藥品供給量及服務範圍內人口數請直接使用csv內的supply、pop欄位。服務範圍內藥局總數量，則請自行計算該藥局服務範圍 (1 km環域) 的藥局總數量。
請從抽樣的藥局中，再篩選出臺南市與雲林縣的藥局(csv內的COUNTY欄位)，利用統計檢定方法，評估臺南市與雲林縣藥局之平均可近性分數是否有顯著差異。

將抽樣的藥局資料對回pharmacy.shp，並篩選出來。計算抽樣藥局服務範圍內的藥局總數，再計算可近性分數。

sample$id=as.character(sample$id)
PH@data=left_join(PH@data,sample)
PH.sample=PH[!is.na(PH$pop),]
PH.sample$count=rowSums(gWithinDistance(PH,PH.sample,dist=1000,byid=T))
PH.sample$access=PH.sample$supply/(PH.sample$pop/PH.sample$count)

從pharmacy.shp資料中選出sample的點位：2點
計算藥局服務範圍內的藥局總數：2點
計算藥局可近性：1點

  1. \(設定顯著水準為{\alpha=0.05}\)
    \(H_0:臺南市與雲林縣藥局之平均可近性分數相同\)
    \(H_a:臺南市與雲林縣藥局之平均可近性分數不同\)
    \(若令臺南市與雲林縣藥局之平均可近性分數分別為\mu_1及\mu_2\)
    \(則H_0:\mu_1=\mu_2，H_a:\mu_1\neq\mu_2\)
    統計假設正確：1點

  2. \(因假設兩地區藥局可近性分數之分布皆呈常態分配，且兩者變異數相等 → 使用t檢定(pooled)即可。\)

two.t=t.test(access~COUNTY,data=PH.sample,var.equal=T)
two.t
## 
##  Two Sample t-test
## 
## data:  access by COUNTY
## t = -3.4093, df = 77, p-value = 0.001039
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.0929080 -0.8121283
## sample estimates:
## mean in group 雲林縣 mean in group 臺南市 
##             7.673326             9.625844
sprintf("t分數=%.4f",two.t$statistic)
## [1] "t分數=-3.4093"

正確統計方法：2點(需用pooled t-test，手算也可)
使用unpooled t-test：1點(忘記寫var.equal=T)

  1. \(計算p值(雙尾)。\)
sprintf("p值=%.4f",two.t$p.value)
## [1] "p值=0.0010"

正確p-value：2點
手算未考慮雙尾或多乘以2：1點

  1. \(p值=0.001<0.05=\alpha，因此拒絕虛無假設H_0。\)
    拒絕虛無假設：1點

  2. \(結論：臺南市與雲林縣藥局之平均可近性分數有顯著地不同。\)
    顯著不同：1點
    臺南市高於雲林縣：0點