使用Tmap繪製地理資料地圖

--- 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套件介紹 ## 範例：高雄輕軌捷運沿線房價研究 ![高雄輕軌捷運沿線房價研究](https://i.imgur.com/0MEHExc.jpg) ```{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) ``` ![](https://i.imgur.com/oyhvkRX.png) + QTM 能快速疊圖 ```{r} qtm(TC)+qtm(tc_mfr) ``` ![](https://i.imgur.com/7z4ErSb.png) ### 主要功能和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> ![](https://i.imgur.com/wvDKNQL.png =60%x) </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()系統顏色會依據欄位而有所變化： ![](https://i.imgur.com/hgJiBTt.png) + 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) # 組圖檢視 ``` ![](https://i.imgur.com/VdqtAB9.png) ```{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的指令 ``` ![](https://i.imgur.com/khgWY3K.png) ```{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) ``` ![](https://i.imgur.com/dfdSs3L.png) + 圖標/比例尺 ```{r} map2 + tm_compass(type="8star",position = c("left","top"))+ # 增加指北針與設定位置 tm_scale_bar(breaks= c(0,10,25), text.size = 0.7) # 增加比例尺 ``` ![](https://i.imgur.com/JtwJuBa.png) + 圖型周圍 ```{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")) # 調整圖例位置 ``` ![](https://i.imgur.com/hCKQ9MZ.png) + 主題風格 ```{r} sty1 = map2+ tm_style("bw" ) # 灰階 sty2 = map2+ tm_style("classic") # 復古 sty3 = map2+ tm_style("cobalt" ) # 深色底白字 sty4 = map2+ tm_style("col_blind" ) # 色盲顏色調整 ``` ![](https://i.imgur.com/Dcy5awm.png) #### 3. 疊圖 ```{r} TCG=tm_shape(TC)+tm_polygons() # 建立底圖 class(TCG) ``` > [1] "tmap" ```{r} TCG + tm_shape(landmark) + # 加入地標資料 tm_dots() + # 以點呈現 tm_dots() ``` ![](https://i.imgur.com/3Q80xyz.png) ```{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") ``` ![](https://i.imgur.com/EB9D1hs.png) ```{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") ``` ![](https://i.imgur.com/0CpuI6v.png) ```{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) ``` ![](https://i.imgur.com/xtLqJ6A.png) #### 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) ``` ![](https://i.imgur.com/GUgk22m.png) + 法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> ![](https://i.imgur.com/T9dJsgJ.png =60%x) </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般也能成為互動式 ``` ![](https://i.imgur.com/ADqEUiM.png) + 同步縮放功能 ```{r} facets = c("pop", "male_ratio","female_ration") tm_shape(TC) + tm_polygons(facets) + tm_facets(nrow = 1, sync = TRUE) ``` ![](https://i.imgur.com/Abuhg70.png)

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.