# 用 R 來畫地形圖——以新竹縣市為例
###### tags: `地理資訊系統` `用R做`
R 可以完成許多資料的視覺化，地理資料的分析與視覺化也可以由 R 來做。如果熟悉 R 的操作，就可以直覺地調整繪圖參數與算繪（rendering），不用在視覺化軟體上拉來拉去，平白消耗電腦效能。
在 R 當中有許多套件可以做簡單的地理資訊處理（geoprocessing），例 [raster](https://cran.r-project.org/web/packages/raster/index.html)、[sf](https://cran.r-project.org/web/packages/sf/index.html)（sf 的[使用簡介](https://r-spatial.github.io/sf/)）；視覺化套件主要有 [tmap](https://cran.r-project.org/web/packages/tmap/index.html)（tmap 的[使用簡介](https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html)）以及搭配 ggplot2 使用的 [ggspatial](https://cran.r-project.org/web/packages/ggspatial/index.html)（ggspatial 的[使用簡介](https://paleolimbot.github.io/ggspatial/articles/ggspatial.html)）。
在這邊我們以 raster 與 sf 分別讀入網格與向量資料，並以 tmap 為視覺化的主要套件來製圖。
這篇用的程式碼參考與修改自[此](https://github.com/GorkyFlorez/Relieve_Brazil/blob/main/Relieve%20de%20Brazil.R?fbclid=IwAR1vTLkvTxw-xWRRAnUOAs7EZ6w7pE53Rf2PW7k8wGYoo3SuuwmDLUHskig)。
:::info
:::spoiler {state="open"}
[TOC]
:::
---
## 取得資料
由政府開放平台取得[直轄市、縣市界線](https://data.gov.tw/dataset/7442)與最新的[網格數值地形模型（DTM）資料](https://data.gov.tw/dataset/160361)（這份資料有挖掉[樂山雷達站](https://zh.wikipedia.org/zh-tw/%E6%A8%82%E5%B1%B1%E9%9B%B7%E9%81%94%E7%AB%99)，可以在[這份資料集](https://data.gov.tw/dataset/35430?fbclid=IwAR2h-T8bP0E4iUWf0OAzR5byuSLgw9482sSFeEnppZqkRYfKj_fADwzKtBY)取得完整的 DEM）。
要注意兩份資料集的座標系統不同。
座標系統的內涵包括大地基準（datum）與座標格式（coordinate format）；大地基準是用作大地測量時，坐標計算的參考依據，而座標格式則是座落位置的不同表示方法。
建立地圖的流程，首先要選擇一個適合的球體／橢球體／基準來近似地球形狀，接著投影到平面，最後展開為地圖，因此人們會依據需求選擇適合的基準與投影方法，搭配合適的座標格式來記錄位置。
縣市界資料集使用 WGS84 基準與經緯度格式（即 [**EPSG:4326**](https://epsg.io/4326)），而
DEM 資料集則使用 TWD97 基準與 121 二度分帶經緯度格式（TWD97 / TM2 zone 121，即[**EPSG:3826**](https://epsg.io/3826)）。
進行疊圖與相關的計算之前，需要將所有圖層轉換成一致的座標系統後，方可操作。
## 讀取資料與地理處理
向量檔 shapefile 使用 `sf` 套件來讀取為 [simple feature](https://en.wikipedia.org/wiki/Simple_Features) 格式（可以想像為點／線／面的地理資料 + 屬性表 attribute table）以利後續縣市資料的抽取，而網格的 DEM 資料，則使用 `raster` 套件讀取。所有圖層的
因為想做大新竹地區的地形圖，所以在讀入縣市邊界資料後，以縣市代號篩出新竹縣市，方法就跟平常處理資料框 dataframe 相同，接著就使用各種地理處理（geoprocessing）方法，首先取新竹縣市邊界的聯集（union），再以這個範圍當作邊界，於 DEM 的網格圖層上切割（clip）出需要的範圍。
針對向量資料的地理處理方法。取自 Saylor Academy 出版的 [Essentials of Geographic Information Systems](https://saylordotorg.github.io/text_essentials-of-geographic-information-systems/index.html)。
```{r}=
border <- st_read("./mapdata202301070205/COUNTY_MOI_1090820.shp") %>%
st_transform(crs = "+init=EPSG:3826")
HC <- st_union(border[border$COUNTYID == c("J", "O"),]) #保留縣市界用st_combine
HC <- as_Spatial(HC)
plot(HC)
DEM <- raster("./不分幅_全台及澎湖/dem_20m.tif", crs = "+init=EPSG:3826")
HC = spTransform(x = HC, CRSobj = proj4string(DEM))
coop <- mask(DEM, HC) #用crop的話會變成剪下HC範圍xlim,ylim的方形
coop <- trim(coop, values = NA)
plot(coop)
plot(HC, add = T)
#writeRaster(coop, "./HCtopomap/DEM_HsinChu_20m_crs.tif", overwrite = TRUE)
```
如果要加入苗栗縣的邊界，若照著上面的方法做，會發現大新竹與苗栗的縣市界在 union 之後，仍然有殘存的邊界。透過一個環域（buffer）的偷吃步，將環域的距離設為 0，就可以完成融合大新竹地區與苗栗縣的輪廓。
操作有時候會出錯，可以加上 `sf_use_s2(FALSE)` 指令，強制運算過程中，將物件從空間幾何圖形視為平面幾何圖形。
```{R}=
HM <- st_union(HC, border[border$COUNTYID == "K",])
sf_use_s2(FALSE)
HM <- st_buffer(HM, 0)
HM <- as_Spatial(HM)
```
## 地圖視覺化
### 分層設色圖
把資料準備好之後，就可以進行疊圖與地圖的視覺化了。視覺化使用 `tmap` 套件進行，也可以用 `ggplot2` 系列套件來繪製。兩個套件都是在 `grid` 圖形系統的基礎上發展而來，所以繪圖指令很類似。
```{R}=
colss <-c("#41641F","#6B9F3A","#EDAE5F","#955235","#9F3B03","#843907","#471600","#4E1C05","#50362F", "#391306", "purple")
HC_map = tm_shape(HC, projection = "+init=epsg:3826") +
tm_borders(lwd = 2) + tm_graticules(n.x = 4, alpha = 0.1, labels.size = 1) +
tm_shape(coop) +
tm_raster(palette = gray(0:10 / 10), n = 100, legend.show = FALSE, alpha = 0.7) +
tm_shape(coop) +
tm_raster(alpha = 0.6, palette = colss ,n = 8, style = "cont",
legend.show = T, title = "Elevation(m)") +
tm_scale_bar(breaks = c(0, 10, 20), text.size = 0.8, text.color = "black", color.dark = "black",
position = c(.01, 0.005), lwd = 1, color.light= "white")+
tm_compass(type = "8star", position=c(.01, .1), size = 7.5, text.color = "black")+
tm_layout( bg.color = "white",
legend.title.size = 1,
legend.position = c("right", "bottom") ,
legend.text.size = 0.8,
fontface = "bold",
legend.format = c(text.align = "right", text.separator = "-"),
inner.margins = c(.1,.1,.1,.22)) +
tm_credits("Elevation Map \n of Hsinchu", position = c(.58, .82), col = "black",
fontface = "bold", size = 2.4, align = "right", fontfamily = "serif")
#儲存圖檔
showtext::showtext_opts(dpi = 1200)
tmap_save(HC_map, paste0("HCmap_", format(Sys.time(), "%y%m%d_%H%M%S"), ".png"),
dpi = 1200, width = 9, height = 7.5)
```
最後可以輸出 1200dpi、9 $\times$ 7 吋的地圖。
有些電腦如果只在 `ggsave` 或者 `tmap_save` 等指令中指定解析度，最後生成的圖片，字會很小；在輸出圖片前，先以 [showtext](https://cran.rstudio.com/web/packages/showtext/vignettes/introduction.html) 套件指定好字型的解析度，就可以避免這種情形，得到與繪圖指令相符的字體大小。
### 加上陰影
如果在地形圖上加入山體陰影（hillshade），可以讓原本的分層設色地形圖，在視覺上更有立體感。
要模擬山體陰影，首先要計算坡度（slope）與坡向（aspect），之後指定太陽的仰角與入射方位角，計算出特定條件下的山體陰影。
```{R}=
slope = terrain(coop , opt = "slope")
aspect = terrain(coop , opt = "aspect")
hill = hillShade(slope, aspect, angle = 45, direction = 135)
```
有了山體陰影的圖層，就可以加入原本的繪圖指令中：
```{R}=
tm_shape(hill) +
tm_raster(palette = gray(0:10 / 10), n = 100, legend.show = FALSE, alpha = 0.7)
```
把原本的 `coop` 改成 `hill`，就可以加入山體陰影，生成視覺上具有立體感的分層設色地形圖了。
其實從繪圖指令中，可以發現這些地圖的生成，就是把不同的地圖圖層調整（或不調整）透明度後依序疊在一起，最後加上圖例、指北針、比例尺等第圖要素，最後加上標題就完成了。地圖依據需要呈現或強調的資料不同，最後透過不同的製圖設計（[cartographic design](https://en.wikipedia.org/wiki/Cartographic_design)），呈現出的地圖美學，之後會嘗試製作不同種類的地圖來介紹。
#### 參考資料與延伸閱讀
1. 主要參考[GorkyFlorez 的程式碼](https://github.com/GorkyFlorez/Relieve_Brazil/blob/main/Relieve%20de%20Brazil.R?fbclid=IwAR01VpkRzctcKUVT2adu2sXd8Eiap9ORTbeDxPKsFAA3A7Gw5xfufpDmUWI)
2. 查詢座標系統的細節：[epsg.io](https://epsg.io/)
3. [世界大地測量系統（WGS）](https://zh.wikipedia.org/zh-tw/%E4%B8%96%E7%95%8C%E5%A4%A7%E5%9C%B0%E6%B5%8B%E9%87%8F%E7%B3%BB%E7%BB%9F)
4. 台師大地理系周學政老師對[座標系統與台灣的座標系統介紹](https://sites.google.com/site/ntnugeo0126/Home/projection/twprojection)
5. [台灣常用的坐標系統及EPSG代碼](https://gis.rchss.sinica.edu.tw/qgis/archives/2823/)
6. 開源地理空間基金會對[台灣大地基準及座標系統的介紹](https://wiki.osgeo.org/wiki/Taiwan_datums)
7. 上河文化的[大地座標系統漫談](https://www.sunriver.com.tw/grid_tm2.htm)
8. [Overview of Coordinate Reference Systems (CRS) in R](https://www.nceas.ucsb.edu/sites/default/files/2020-04/OverviewCoordinateReferenceSystems.pdf)
9. Dennis 的 [[R語言專題]用ggplot畫地圖 – 基礎篇](https://r-lover.com/tutorial/r-seminar/ggplot-map-basic/)
10. 台大昆蟲系邱名鍾老師的[繪圖應用：地圖繪製與物種分布](/TN9b5fzLTWiLqIOQW6qTGg?view)
11. 張安瑜的 [Vector 的操作](https://kemushi54.github.io/R-basic/sp_and_rgdal.html)、
[用 R 建構網格系統](https://kemushi54.github.io/R-basic/fishnet.html)
12. 台大地理系溫在弘老師[《空間分析：方法與應用》書中的範例程式碼](https://wenlab501.github.io/SpatialAnalysis/)、
[空間分析課程的資料](https://wenlab501.github.io/GEOG2017/)
13. Martijn Tennekes & Jakub Nowosad (2021). [Elegant and informative maps with tmap](https://r-tmap.github.io/tmap-book/).
14. Rita Tang 的[使用Tmap繪製地理資料地圖](/1Xpfqu4ITqOoPTY-EOmWqA)
15. 舒乐乐，孟宪红，陈昊，李照国，赵林（2023）。[R在地球科学中的应用](https://www.shud.xyz/bookr/)。
16. 國立臺北大學經濟學系林茂廷老師的[經濟資料視覺化處理](https://bookdown.org/tpemartin/108-1-ntpu-datavisualization/)
<span style="font-size:30px!">🐕🦺</span><font color="dcdcdc">2023.12.24</font>
Sign in with Wallet
Connect another wallet