1
使用 R 軟體 繪製統計地圖的介紹與應用(1)
蔡靜雯 副統計分析師 地區性的統計資料結合地圖方式呈現,可以更容易觀察到不同區域間的變化, 和是否有地埋位置的影響。本期的 eNews 內容,分成 2 大部分,第 1 部分介紹 地圖的相關資料以及 shp 地圖檔案格式。第 2 部分則是使用台灣地圖檔案和戶 政各地的統計資料,使用 R 軟體實際操作,繪製台灣統計地圖。第1部分 地圖資料介紹
1-1、坐標系統 目前台灣常見的地圖座標系統有以下3種:• TWD67 - 採用國際上Geodetic Reference System 1967 橢球體, 作為臺灣地區 參考橢球體,以南投埔里虎子山為大地基準。
• TWD97 - 採用國際上Geodetic Reference System 1980 橢球體, 作為臺灣地區 參考橢球體,以八個衛星追蹤站為大地基準,1997年完成以 GPS 重 新計算坐標基準。 • WGS84 - 美國國防部於西元1984採用新的地球標準物理模型,以地球的質量 中心為中心點,加上了分布在全世界各地的1500個地理座標參考點, 為GPS而制定的世界座標系統,目前為國際通用的經緯度坐標。 TWD67 是與該地區最密合的「區域性」大地基準,只適用於台灣地區的座標系統, 與 WGS84 的坐標基準相差約1公里,TWD97 的坐標基準則與 WGS84 的坐標基 準差異不大,相差約幾公分至數十公分。目前線上的地圖系統幾乎都使用WGS84的 座標系統(Google Map, UrMap, Google Earth)。
2 1-2、地圖資料格式
ESRI Shapefile (shp),或簡稱shapefile,是美國環境系統研究所公司(ESRI)開發的 空間資料開放格式。一個 Shapefile 檔案通常包含多個檔案,同組地圖資料其主檔 名應該相同,如下圖範例所示,相同主檔名代表為同一組資料, gadm36_TWN_0、 gadm36_TWN_1、gadm36_TWN_2 分別為三組地圖資料。 Shapefile 檔案中有三個必須檔案,其副檔名和儲存的資料內容,分別為 - .shp:圖形格式,用以儲存資料的圖形, - .shx:圖形索引格式,記錄圖形的位置, - .dbf:屬性資料格式,儲存每個圖形的屬性資料。 如果要讀取 shapefile 檔,以上三個檔案都必須同時具備,否則會無法讀取 成功。 有時可能還會有其它非必須的檔案: - .prj:投影格式,用於儲存地理座標系統與投影資訊., - .sbn:空間索引檔,可以加快開啟shapefile的速度。.
3
第2部分 繪製統計地圖
2-1、台灣地圖資料下載
以下兩個地方的網址都可以下載台灣的地圖資料
1. GADM maps and data https://gadm.org/download_country_v3.html
GADM是全球行政區域資料庫,採用座標系統 WGS84,提供了所有國家及 其行政區域細分的地圖和空間資料,可免費用於學術和其他非商業用途。未經 GADM事先許可,不得用於商業用途。 下載步驟:點選網址後,在Country欄位選取Taiwan或直接輸入Taiwan,點選 Shapefile 即可下載台灣地圖檔案。 下載的壓縮檔,共有16個檔案,除了 license.txt 外有3種檔案名稱,分別代表3組 不同區域劃分的地圖資料: gadm36_TWN_0 - 無區分縣市,只有台灣形狀, gadm36_TWN_1 - 區分六都區域, gadm36_TWN_2 - 區分22個縣市區域。 從壓縮檔內容,也同時確認繪製地圖必須的3個檔案 .shp、.shx 和 .dbf 檔都有 包含在其中。
4 2. 鄉鎮市區界線(TWD97經緯度) https://data.gov.tw/dataset/7441 為我國各鄉(鎮、市、區)行政區域界線圖資,採用座標系統 TWD97 由內 政部國土測繪中心提供,在政府資料開放平臺提供免費下載,可直接點網址下 載檔案,或依照以下路徑:政府資料開放平臺->公共資訊->內政部國土測繪中 心->鄉鎮市區界線(TWD97經緯度),找到該檔案的所在位置。除此之外,也有其 它台灣相關的地圖資料可下載。 下載步驟:點選網址頁面的 SHP 下載。
5
下載的壓縮檔,共有7個檔案,其中也確認有包含繪製地圖必須的3個檔案 .shp、.shx 和 .dbf 檔。
6 2-2、讀入地圖資料檔案 以下範例使用 GADM 區分22縣市的地圖檔案 gadm36_TWN_1。 [程式碼] - 讀入 shp 檔案 install.packages("rgdal") #安裝 讀取shp檔案的套件 library(rgdal) taiwan_shap2=readOGR("C:/Users/CWTsai/Desktop/plot -繪製地圖/gadm 36_TWN_shp/gadm36_TWN_2.shp") #讀入shp檔案 - output - 顯示此檔案有22筆資料,13個欄位 [程式碼] - 使用 head 函數,查看前6筆資料的內容(位於地圖檔案的data屬性裡) head(taiwan_shap2@data) - output - 欄位名稱NL_NAME_1、VARNAME_2、NL_NAME_2、TYPE_2的資料呈現亂 碼,需進一步處理亂碼問題。
7
[程式碼] - 使用 iconv 函數,對亂碼欄位執行轉換
taiwan_shap2@data$NL_NAME_1 <- iconv(taiwan_shap2@data$NL_NAME_ 1, from = "UTF-8", to="UTF-8")
taiwan_shap2@data$VARNAME_2 <- iconv(taiwan_shap2@data$VARNAME_ 2, from = "UTF-8", to="UTF-8")
taiwan_shap2@data$NL_NAME_2 <- iconv(taiwan_shap2@data$NL_NAME_ 2, from = "UTF-8", to="UTF-8")
taiwan_shap2@data$TYPE_2 <- iconv(taiwan_shap2@data$TYPE_2, from = "UTF-8", to="UTF-8")
head(taiwan_shap2@data)
- output - 主要的中文縣市名稱欄位 NL_NAME_1 已成功轉換成繁體中文,雖然 VARNAME_2 欄位還是有亂碼情況,但因後續不會用到此欄位,可忽略 亂碼問題。
8 2-3、繪製地圖
[程式碼] - 使用ggplot() 函數繪製台灣地圖,其中,x 放經度(long),y 放緯度(lat), group 為22個縣市 install.packages("ggplot2") #安裝 讀取shp檔案的套件 library(ggplot2) ggplot(taiwan_shap2, aes(x=long,y=lat,group=group))+geom_path() - output - 區分22個縣市的台灣地圖 2-3、統計資料檔案 使用內政統計指標縣市排名資料https://www.moi.gov.tw/stat/node.aspx?sn=6718, 分別下載100~106年的縣市指標檔案,採用戶政檔案裡的出生登記嬰兒性比例資料 (性比例係指每百個女嬰的男嬰數量),再匯整在一起存成CSV檔。
9
原始檔案內容:
10
[程式碼] - 讀入匯整完成的嬰兒性比例統計資料檔案 index.csv
index.data=read.csv("C:/Users/CWTsai/Desktop/plot -繪製地圖/政府開 放資料/縣市指標歷年資料/index.csv", head=T)
View(index.data)
- output - 縣市欄位名稱city,bsex_ratio100 ~ bsex_ratio106 欄位為100~106年嬰兒 性別比例 2-4、地圖資料和統計資料整理合併 地圖資料(taiwan_shap2)的縣市名稱和嬰兒性比例統計資料( index.data )的縣市名 稱要一致,才能將資料合併在一起。 [程式碼] - 地圖資料(taiwan_shap2)的縣市名稱 taiwan_shap2$NL_NAME_2 - output - 地圖資料的縣市名稱欄位 NL_NAME_2 屬性為文字字串,名稱頭尾有雙 引號 [程式碼] - 嬰兒性比例統計資料( index.data )的縣市名稱 index.data$city - output - 嬰兒性比例資料的縣市名稱欄位city屬性為因子變數(factor),Levels顯示 有幾種分類變數
11 雖然在合併檔案時,factor屬性欄位和字串屬性欄位合併沒有問題,但建議還是將 欄位屬性轉成一致再合併。 [程式碼] - 地圖資料的縣市名稱和嬰兒性比例資料的縣市名稱整理成一致 #另外建立一個 zh_NAME 欄位,存放修正後的地圖資料縣市名稱 taiwan_shap2$zh_NAME=taiwan_shap2$NL_NAME_2 #將地圖和嬰兒性比例資料的5個不一致縣市名稱,重新命名與嬰兒性比例的縣市名稱 一樣 taiwan_shap2$zh_NAME[taiwan_shap2$zh_NAME=="馬祖列島"]="連江縣" taiwan_shap2$zh_NAME[taiwan_shap2$zh_NAME=="台中"]="臺中市" taiwan_shap2$zh_NAME[taiwan_shap2$zh_NAME=="台南"]="臺南市" taiwan_shap2$zh_NAME[taiwan_shap2$zh_NAME=="台北市"]="臺北市" taiwan_shap2$zh_NAME[taiwan_shap2$zh_NAME=="台東縣"]="臺東縣" #將嬰兒性比例資料的縣市名稱欄位city屬性轉成文字字串 index.data$city=as.character(index.data$city) #排序後print修正後的地圖和嬰兒性比例資料,確認縣市名稱是否一致 [程式碼] - print修正後的地圖資料其縣市名稱欄位 zh_NAME,排序sort是為了方便 比對 sort(taiwan_shap2$zh_NAME) - output - [程式碼] - print修正後的嬰兒性比例資料其縣市名稱欄位city sort(index.data$city) - output -
12
[程式碼] - 合併地圖資料和嬰兒性比例資料
taiwan_shap2@data=merge(taiwan_shap2@data, index.data, by.x="zh_ NAME", by.y="city", sort=F)
head(taiwan_shap2@data) - output - 100~106年的嬰兒性比例資料與地圖資料成功合併在同一個檔案 2-5、繪製合併後的統計地圖資料 2-5-1、方法1 - 使用qtm() 函數 install.packages("tmap") #安裝 qtm 函數案的套件 library(tmap)
qtm(taiwan_shap2,fill="bsex_ratio100", text="zh_NAME", text.size =0.7, fill.title="birth sex ratio", fill.palette="Reds")
13
- output -
[程式碼] - 100年和106年的嬰兒性比例統計地圖
qtm(taiwan_shap2,fill="bsex_ratio100", fill.title="birth sex rat io(100)")
qtm(taiwan_shap2,fill="bsex_ratio106", fill.title="birth sex rat io(106)")
- output - qtm()函數會根據嬰兒性比例的數值分佈,自動分類比例組別和給定組別 顏色,100年和106年的圖例,雖然性別比例都是以每 5% 分組,但給定 的組別顏色不一樣,同時放一起比較時,會因圖例不一致而不易比較。
14
2-5-2、方法2 – plot() 函數 調色盤介紹
Palettes {grDevices} 套件裡,有rainbow、heat.colors、terrain.colors、topo.colors、 cm.colors 5種類型的調色盤。
[程式碼] - 呈現5種類型的調色盤顏色 #Display a Color Image
display.col=function(color,title){
image(1:10,1,as.matrix(1:10),col=color,xlab=title, ylab="",xaxt="n",yaxt="n",bty="n") } par(mfcol=c(5,1), mai=c(0.1,0.1,0.5,0)) display.col(rainbow(10),"rainbow"); display.col(heat.colors(10),"heat.colors"); display.col(terrain.colors(10),"terrain.colors"); display.col(topo.colors(10),"topo.colors"); display.col(cm.colors(10),"cm.colors" ) - output - 由上至下分別為rainbow、heat.colors、terrain.colors、topo.colors、cm.colors 5種調色盤的顏色
15 使用plot() 函數,繪製統計地圖步驟,依序程式碼1~程式碼3, [程式碼 1] - 建立繪圖色盤,100 年~106 年的性別比例介於 85%~135% 之間,每 5 個百分比為一組,共分成 10 組,指定繪圖色盤有 10 種顏色。 my.color=rev(heat.colors(10)) [程式碼 2] - 建立一致性的性別比例資料分類,100 年~106 年的性別比例介於 85%~135% 之間,每 5 個百分比為一組,並從建立的繪圖色盤給每 一種比例組別一個固定的顏色 ratio.color100=my.color[cut(taiwan_shap2@data$bsex_ratio100, bre aks=seq(85,135,5),labels=F)] [程式碼 3] - 繪製性別比地圖 #設定畫布邊界 par(mai=c(0.1,0.1,0.5,0)) #繪製統計地圖 plot(taiwan_shap2, col=ratio.color100) #在22個縣市區域標註縣市名稱 text(coordinates(taiwan_shap2)[,1], coordinates(taiwan_shap2)[, 2], taiwan_shap2$zh_NAME, cex=0.6, col="gray5")
#圖例設定
legend("bottomleft",
title="birth sex ratio(100)",
legend=c("85~90", "90~95", "95~100", "100~105", "105~110", "110~115", "115~120", "120~125", "125~130", "130~135"), fill=my.color,
16
- output -
[程式碼] - 呈現 100~106 年 縣市性別比例地圖 attach(taiwan_shap2@data)
#將繪製統計地圖步驟(程式碼1~程式碼3)寫成function,年度為參數 ratio.plot=function(year){
bsex_ratio=paste("bsex_ratio", year, sep="") #建立繪圖色盤
ratio.color=my.color[cut(get(bsex_ratio), breaks=seq(85,135,5),l abels=F)]
#繪製統計地圖
plot(taiwan_shap2, col=ratio.color, main=paste(year, "年 出生性別比 ", sep="") )
#圖例設定
legend("bottomleft",
title=paste("birth sex ratio", sep="") ,
17 "110~115", "115~120", "120~125", "125~130", "130~135"), fill=my.color, cex=0.8) } #設定畫布邊界 par(mai=c(0.1,0.1,0.5,0), mfrow=c(2,4)) #將繪製統計地圖步驟的function ratio.plot(year),使用迴圈,畫出100年~1 06年的統計地圖 for (year in 100:106){ ratio.plot(year) } - output - 每個年度的性別比例分組和指定的分組顏色都一致,就可直接從顏色觀 察 22 個縣市,從 100 年~106 年的出生嬰兒性別比例的變化。
18 2-5-3、方法3 – ggplot() 函數 先將 .shp 檔案中的地區屬性資料和經緯度資料分別抓出,再合併在一起 [程式碼] - 抓出 .shp 檔案中的地區屬性資料(包含地區名稱和先前已合併的統計 資料) taiwan_shap2.data=taiwan_shap2@data head(taiwan_shap2.data) output
-19 [程式碼] - 在taiwan_shap2.data(地區屬性資料)中增加一個 id 欄位(用於合併經緯 度資料) taiwan_shap2.data$id=as.character(0:21) head(taiwan_shap2.data) output -[程式碼] - 使用 fortify() 函數將 .shp 檔案中的經緯度抓出來 taiwan_shap2.latitude_ongitude=fortify(taiwan_shap2) head(taiwan_shap2.latitude_ongitude)
20 [程式碼] - 合併 taiwan_shap2.data(地區屬性+統計資料) 和 taiwan_shap2.latitude_ongitude (經緯度資料) taiwan_mapdata=merge(taiwan_shap2.data, taiwan_shap2.latitude_on gitude, by="id")
ggplot(taiwan_mapdata, aes(x=long, y=lat, group=group, fill=bsex_ ratio100 )) + geom_polygon(colour="black")
-21 2-6、使用鄉鎮市區界線(TWD97經緯度)地圖資料,依照前面執行步驟,繪製台灣368 個鄉鎮區域地圖 [程式碼] – 讀入縣市鄉鎮地圖檔案 town_shap=readOGR(dsn="C:/Users/CWTsai/Desktop/plot -繪製地圖/鄉鎮 地圖資料/TOWN_MOI_1071226.shp") - output - 縣市鄉鎮地圖檔案中,有368筆資料(368個鄉鎮),7個欄位 [程式碼] # 中文亂碼轉換 iconv {base} town_shap@data$COUNTYNAME=iconv(town_shap@data$COUNTYNAME , from = "UTF-8", to="UTF-8")
town_shap@data$TOWNNAME=iconv(town_shap@data$TOWNNAME , from = "U TF-8", to="UTF-8")
head(town_shap@data) # 中文正常顯示 output
-22 [程式碼] #繪製縣市鄉鎮地圖-方法1 qtm() 函數 library(tmap) qtm(town_shap) # 未給定鄉鎮顏色 output
-23 [程式碼] - #以鄉鎮市區界線(TWD97經緯度)地圖檔案中,現有的22個縣市 COUNTYCODE欄位編 碼,做為fill的指定縣市顏色 qtm(town_shap, fill="COUNTYCODE") output
-24
[程式碼] -
#以鄉鎮市區界線(TWD97經緯度)地圖檔案中,現有的365個鄉鎮 TOWNCODE 欄位編 碼,做為fill的指定鄉鎮顏色
qtm(town_shap, fill="TOWNCODE") #fill最多只能指定30種顏色,分類超過 30個組別會自動合併分成30組
-25
[程式碼]
#繪製縣市鄉鎮地圖-方法2 plot()函數 plot(town_shap) #未給定鄉鎮區域顏色 output
-plot(town_shap, col=rainbow(365)) #使用rainbow指令給定鄉鎮區域顏色 output
-26 參考資料 坐標系統相關 • Taiwan datums https://wiki.osgeo.org/wiki/Taiwan_datums#TWD97_vs._WGS84 • 上河文化 – 大地坐標系統漫談 http://www.sunriver.com.tw/grid_tm2.htm • 中華民國內政部國士測繪中心 – 平面控制 https://www.nlsc.gov.tw/Home/MakePage/42?level=42 • 台灣常用坐標系統簡介 https://blog.xuite.net/lwkntu/blog/364536963-%E5%8F%B0%E7%81%A3%E5%B8%B8%E7%94%A8%E5%9D%90%E6%A8%99%E7 %B3%BB%E7%B5%B1%E7%B0%A1%E4%BB%8B • 坐標系統 Coordinate systems http://140.121.160.124/GEO/%E5%BA%A7%E6%A8%99%E7%B3%BB%E7%B5%B1 .pdf shp 檔案 • 維基百科 Shapefile https://zh.wikipedia.org/wiki/Shapefile • Ultimanual 的 ArcGIS https://ultimanualarcgis.wordpress.com/2016/06/03/2-fu-dang-ming-de-yi-yi/ 繪製地圖相關 • qtm function | R Documentation https://www.rdocumentation.org/packages/tmap/versions/2.2/topics/qtm • 深入淺出繪製統計地圖 4 – R https://ariheart2011.wordpress.com/2017/05/10/%E6%B7%B1%E5%85%A5%E6% B7%BA%E5%87%BA%E7%B9%AA%E8%A3%BD%E7%B5%B1%E8%A8%88%E5%9C
27
%B0%E5%9C%96-4-r/
• 主題式地圖(Thematic map) - 政府開放資料為例