tags: 計量經濟 地理資訊繪圖 Econometrics Visualization tmap sf Map Geometry

使用Tmap繪製地理資料地圖

地理資料的構成

1. 點：每一點(如房價、地標)
2. 線：線段(如道路、河流)
3. 面：區域(如村里、工業區等)
4. 網格：多使用於地形、自然資源，由每一網格所構成的一張圖(類似解析度)

而地理資料若加上當地區域或與四大資料有關的外接資料時，此時所能呈現在地圖上的效果將非常強大，而今日主角tmap便是讓功能放大的好套件。

TGOS(地址轉經緯度服務平台)

• 申請帳號後，可以申請批次地址轉經緯度服務
• 批准後會得到一個API Key
• 再到批次門牌地址比對服務貼上API Key
• 依照範本將地址資料整理成csv檔並上傳
• 即可進行批次比對

Tmap套件介紹

範例：高雄輕軌捷運沿線房價研究

# 台灣地理經緯度資料
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檔讀進來)

讀取資料

#shapefile檔案讀法  #sf packege
tt = st_read("路徑/檔案名稱.shp")

#raster檔案讀法  #raster packege
tair = raster("路徑/檔案名稱.tif")

併入地理資訊

• 資料型態：sf data frame
  • 定義了地理資料邊界(多邊形)

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

通常適用於快速看一下圖形狀況時的便捷寫法
能快速看各欄分布狀況

plot(TC,max.plot = 17)

• QTM

能快速疊圖

qtm(TC)+qtm(tc_mfr)

主要功能

和ggplot類似，需要基底(tm_shape)，且可再疊上更多圖層元素(tm_fill or tm_dots)和地圖的邊界、美編與標示等。

1. 地圖基底

mapA=tm_shape(TC) +  # 給予目標圖層
	tm_fill() +      # 填滿圖層工具
	tm_borders()     # 加邊界
	
# 可以用tm_polygon()取代tm_fill() + tm_borders(), 一次加兩個基底

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

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) # 組圖檢視

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的指令

# 類別色塊
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)

• 圖標/比例尺

map2 +
  tm_compass(type="8star",position = c("left","top"))+ # 增加指北針與設定位置
  tm_scale_bar(breaks= c(0,10,25), text.size = 0.7)    # 增加比例尺

• 圖型周圍

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"))      # 調整圖例位置

• 主題風格

sty1 = map2+ tm_style("bw" )             # 灰階
sty2 = map2+ tm_style("classic")         # 復古
sty3 = map2+ tm_style("cobalt" )         # 深色底白字
sty4 = map2+ tm_style("col_blind" )      # 色盲顏色調整

3. 疊圖

TCG=tm_shape(TC)+tm_polygons()             # 建立底圖
class(TCG)

[1] "tmap"

TCG +
  tm_shape(landmark) +  # 加入地標資料     
  tm_dots() +           # 以點呈現 tm_dots()

#將剛剛所學來畫幾圖吧
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")

#台中地形與土石流潛勢區分布
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")

#台中市各區平均收入
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()

#當需要看年分變化時，可使用年份分法的地圖繪製
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

#將台中、台南、台北畫在同一張圖
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)

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))

互動式地圖

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般也能成為互動式

• 同步縮放功能

facets = c("pop", "male_ratio","female_ration")
tm_shape(TC) + tm_polygons(facets) + 
  tm_facets(nrow = 1, sync = TRUE)