---
GA: UA-159972578-2
---
###### tags: `計量經濟` `地理資訊繪圖` `Econometrics` `Visualization` `tmap` `sf` `Map` `Geometry`
# 使用Tmap繪製地理資料地圖
## 地理資料的構成
1. 點：每一點(如房價、地標)
2. 線：線段(如道路、河流)
3. 面：區域(如村里、工業區等)
4. 網格：多使用於地形、自然資源，由每一網格所構成的一張圖(類似解析度)
而地理資料若加上當地區域或與四大資料有關的外接資料時，此時所能呈現在地圖上的效果將非常強大，而今日主角tmap便是讓功能放大的好套件。
## TGOS(地址轉經緯度服務平台)
+ 申請帳號後，可以申請[批次地址轉經緯度服務](https://www.tgos.tw/tgos/Web/Service/Apply/TGOS_Apply_Service.aspx?keyword=%E6%89%B9%E6%AC%A1%E5%9C%B0%E5%9D%80%E6%AF%94%E5%B0%8D%E6%9C%8D%E5%8B%99)
+ 批准後會得到一個API Key
+ 再到[批次門牌地址比對服務](https://www.tgos.tw/tgos/Web/Address/TGOS_Address.aspx)貼上API Key
+ 依照範本將地址資料整理成csv檔並上傳
+ 即可進行批次比對
# Tmap套件介紹
## 範例：高雄輕軌捷運沿線房價研究
```{r}
# 台灣地理經緯度資料
load("geotw108.rdata") # twCounty：縣市資料 # twTown：鄉鎮市區資料 # twVill：里資料
county = st_simplify(twCounty,dTolerance=100) # 將地圖邊界變得較為平滑
# 高雄捷運輕軌線資料
L = st_read("kaoh/light_rail/LRT_1090102.shp", crs=3826) # 是"線"的要用3826
L = subset(L,LRTSYS=="高雄捷運") # 選取高雄捷運的路線
st_crs(L) # 檢查坐標系轉換是否成功
L = L %>% mutate_if(is.factor,as.character)
# 房產時價登錄資料
load("kaoh/house_price.rdata") # 在資料框中加入geometry的資訊
A = st_as_sf(D, coords=c("E","N"), crs=4326) %>% # "E"東經 "N"北緯
st_transform(crs=3826) # 當資料本身是經緯度時，不能直接轉3826，要先轉4326
# 計算房子到輕軌路線的距離
dx = st_distance(A, L) # 房子到輕軌的距離
dim(dx) # 49898筆資料，2條輕軌線
A$dx = apply(dx,1,min) # 房子到其中一條輕軌線的距離
range(A$dx) %>% round # 距離介於0-87104
# 篩選資料
Houses = A %>% filter(dx < 2000) %>% sample_n(1000) # 由於房屋的資料相當龐大，故進行抽樣
# 找出有輕軌穿過的行政區
Town = twTown[st_intersects(st_union(L), twTown)[[1]], ] # 回傳match到的行政區的index # st_union 將兩條輕軌線合併成一條
# 畫圖
subset(twVill, TOWN_ID %in% Town$TOWN_ID) %>%
tm_shape()+
tm_polygons("P_DEN", palette = "Greens", title="人口密度")+
tm_shape(Houses) +
tm_dots(col='type',size=0.7, alpha=0.8 ,n=8, title="房屋類型", palette = "Set3") +
tm_shape(L) +
tm_lines(col='LRTID',lwd=10,palette=c("green","orange"), title.col="輕軌路線") +
tm_shape(Town) + tm_borders(col="#925552", lwd = 4) +
tm_text(text="TOWN",col="#925552",size = 3) +
tm_layout(fontfamily="蘋方-繁 中黑體",
legend.title.size = 2,
legend.text.size = 1.3)
```
## 地理資料前置處理
### 資料介紹
+ X：各村里所得資料
+ D：各村里人口資料
+ X1 & D1：上述資料合併為區資料
+ T：各區所得與人口資料(即X+D)
+ TY：各區所得與人口資料(有分年分)
+ tt：各區界圖 (shapefile檔讀進來)
+ mfriver：台中土石流潛勢溪圖 (shapefile檔讀進來)
+ landmark：台中地標 (shapefile檔讀進來)
+ tair：台灣網格資料 (raster檔讀進來)
### 讀取資料
```{r}
#shapefile檔案讀法 #sf packege
tt = st_read("路徑/檔案名稱.shp")
#raster檔案讀法 #raster packege
tair = raster("路徑/檔案名稱.tif")
```
### 併入地理資訊
+ 資料型態：sf data frame
+ 定義了地理資料邊界(多邊形)
```{r}
pacman::p_load(purrr,dplyr,ggplot2,tidyr,plotly,stringr,googleVis,sf,sp,tmap,grDevices,raster,RColorBrewer)
# 土石流潛勢溪資料
mfr = mfriver %>% st_transform(,crs=4326) # 將區界資料轉成wgs84的坐標系
st_crs(mfr) # 確認坐標系是否轉成功
class(mfr) # 確認是否轉成sf object
any(duplicated(tt1$TOWNCODE)) # 檢查是否有重複值
```
```
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
```
> [1] "sf" "data.frame"
> [2] FALSE
+ Coordinate Reference System(CRS)座標投射方式
+ crs=4326
+ 方法：wgs84
+ 適用於經緯度資料
+ crs=3826:
+ 方法：TWD97
+ 適用於台灣地圖(計算兩地平面投影距離較準確)
+ 例如：探討瓦斯管線和房價的關係
+ 如果資料有經緯度(Global)資料，應該先把它用wgs84投射，將兩種資料合併後，再一起轉換成各國(Local)的平面地圖(經緯度的投射其實是球型，非平面)
## 繪製地圖(tmap)
### 小熱身
+ plot
通常適用於快速看一下圖形狀況時的便捷寫法
能快速看各欄分布狀況
```{r}
plot(TC,max.plot = 17)
```
+ QTM
能快速疊圖
```{r}
qtm(TC)+qtm(tc_mfr)
```
### 主要功能
和ggplot類似，需要基底(tm_shape)，且可再疊上更多圖層元素(tm_fill or tm_dots)和地圖的邊界、美編與標示等。
#### 1. 地圖基底
```{r}
mapA=tm_shape(TC) + # 給予目標圖層
tm_fill() + # 填滿圖層工具
tm_borders() # 加邊界
# 可以用tm_polygon()取代tm_fill() + tm_borders(), 一次加兩個基底
```
<center>
</center>
#### 2. 地圖美學
+ tm_fill() 顏色填法(style)介紹
1) cat：類別值，顏色各自特別
2) fixed：自訂區間
3) sd：標準差分法
4) pretty：從0開始，四捨五入(預設)
5) quantile：四分位數
6) jenks：抓出數值相近，並分出區間差異
7) log10_pretty：取對數
8) cont：在連續色塊呈現大量變色，適合網格資料
9) equal：分開均勻分布的變量
+ 調色盤(palette)概念
Tmap()系統顏色會依據欄位而有所變化：
+ tm_layout() 可以設定字型
+ fontfamily="蘋方-繁 中黑體" (mac畫中文的好幫手~)
+ more detail: https://rdrr.io/cran/tmap/man/tm_layout.html
```{r}
map1 = tm_shape(TC) + # 給予目標圖層
tm_fill(col = "pop", # col 為顏色依據(資料或實際顏色)
alpha = 0.7, # 調整色彩濃度
style = "quantile", n=5, # style 可以調整圖例顏色區間規則(呈現的風格)，預設為"pretty"
palette = "Oranges", # 調色盤
title = "Population") # 為圖例命名
map2 = tm_shape(TC) +
tm_fill(col = "pop",
style = "fixed",breaks= c(1.5e5,2.5e5,3.5e5,4.5e5,5.5e5,6.5e5,7.5e5,8.5e5,9.5e5,1e6),
palette = "Oranges",
title = "Population")
tmap_arrange(map1, map2) # 組圖檢視
```
```{r}
mapx=tm_shape(TC)+
tm_borders(col ="grey40", # col 可設定顏色
lwd=2, # lwd 可設定邊界粗細
lty=2) # lty 可設定虛線
mapy=tm_shape(TC)+
tm_polygons("pop",
style = "quantile", # 自訂風格
n= 5, # 設定風格區間次數
palette = "BuGn", # 調色盤
title = "Population", # 新增標題
lwd=2,
lty=2) # tm_polygons 可包含tm_fill()與tm_borders的指令
```
```{r}
# 類別色塊
colA=tm_shape(TC) + tm_polygons("county")
# 連續色塊
colB=tm_shape(TC) + tm_polygons("pop", palette = "Blues")
# 顏色二分色塊
red_gray = brewer.pal(5,"RdGy") # 建立二分法的調色盤
colC=tm_shape(TC) + tm_polygons("pop", palette = red_gray)
```
+ 圖標/比例尺
```{r}
map2 +
tm_compass(type="8star",position = c("left","top"))+ # 增加指北針與設定位置
tm_scale_bar(breaks= c(0,10,25), text.size = 0.7) # 增加比例尺
```
+ 圖型周圍
```{r}
deco1= map2 + tm_layout(title = "Taichung Population") # 加上地圖標題
deco2= map2 + tm_layout(scale = 5) # 加粗所有框線
deco3= map2 + tm_layout(bg.color = "grey") # 背板上色
deco4= map2 + tm_layout(frame = F) # 移除外框
deco5= map2 + tm_layout(inner.margin = 0.2) # 圖型縮小
deco6= map2 + tm_layout(legend.position=c("left","bottom")) # 調整圖例位置
```
+ 主題風格
```{r}
sty1 = map2+ tm_style("bw" ) # 灰階
sty2 = map2+ tm_style("classic") # 復古
sty3 = map2+ tm_style("cobalt" ) # 深色底白字
sty4 = map2+ tm_style("col_blind" ) # 色盲顏色調整
```
#### 3. 疊圖
```{r}
TCG=tm_shape(TC)+tm_polygons() # 建立底圖
class(TCG)
```
> [1] "tmap"
```{r}
TCG +
tm_shape(landmark) + # 加入地標資料
tm_dots() + # 以點呈現 tm_dots()
```
```{r}
#將剛剛所學來畫幾圖吧
TC1=TC[,-c(1:2)]
#台中人口分布
tm_shape(TC1)+
tm_symbols(col = "red", border.col = "white", size = "pop")+
tm_borders(col = "grey60",lwd = 0.5)+ tmap_mode("plot")
```
```{r}
#台中地形與土石流潛勢區分布
TCH = TC %>% filter(county=="和平區")
tm_shape(TCH)+tm_polygons() +
tm_shape(tair)+ # 加入網格資料(tm_raster())
tm_raster(alpha=0.3) +
tm_shape(tc_mfr) + # 加入土石流潛勢溪圖
tm_lines(col="Potential",palette = "Blues")
```
```{r}
#台中市各區平均收入
tm_shape(TC1) +
tm_fill(col="Income_mean",
palette = "Greens" ,
style = "quantile",n=7 ,
title = "台中平均收入",
alpha = 0.5) +
tm_borders(col = "grey60",lwd = 0.5) +
tm_text("county",size = 0.6)
```
#### 4. 多圖呈現
+ 法1：tmap_arrange()
+ 法2：facet()
```{r}
#當需要看年分變化時，可使用年份分法的地圖繪製
TYC_100_103= TYC %>% filter(year %in% c(100,101,102,103))
tm_shape(TC)+tm_polygons()+
tm_shape(TYC_100_103)+
tm_symbols(col = "red", border.col = "white", size = "pop") + #size為以大小來顯示(也可用tm_bubble)
tm_facets(by = "year", # by 變量來分圖
nrow = 2, # 以幾列來呈現(以可設定ncol))
free.coords = F)
```
+ 法3：grid::viewport
```{r}
#將台中、台南、台北畫在同一張圖
TC_map = tm_shape(TC, projection = 4326) + tm_polygons() +
tm_layout(frame = FALSE,inner.margin = 0.2)
TN_map = tm_shape(TN) + tm_polygons() +
tm_layout(title = "Tainan", frame = FALSE, bg.color = NA,
title.position = c("LEFT", "BOTTOM"))
TP_map = tm_shape(TP) + tm_polygons() +
tm_layout(title = "Taipei", frame = FALSE, bg.color = NA)
```
```{r}
library(grid)
TC_map
print(TN_map, vp = grid::viewport(0.5, 0.25, width = 0.25, height = 0.25))
print(TP_map, vp = grid::viewport(0.25, 0.2, width = 0.3, height = 0.3))
```
<center>
</center>
### 互動式地圖
```{r}
subset(twVill, TOWN_ID %in% Town$TOWN_ID) %>%
tm_shape()+
tm_polygons("P_DEN", palette = "Greens", title="人口密度", alpha = 0.5)+
tm_shape(L) +
tm_lines(col='LRTID',lwd=3,palette=c("green","orange")) +
tm_shape(Houses) +
tm_dots(col='age',size=0.05, alpha=0.9 ,n=8, title="屋齡", palette = "RdPu") +
tm_shape(Town) + tm_borders(col="gray") +
tm_text(text="TOWN",col="black" ) +
tmap_mode("view") # 讓地圖如ggplotly般也能成為互動式
```
+ 同步縮放功能
```{r}
facets = c("pop", "male_ratio","female_ration")
tm_shape(TC) + tm_polygons(facets) +
tm_facets(nrow = 1, sync = TRUE)
```
Sign in with Wallet
Connect another wallet