用 R 來畫地形圖——以新竹縣市為例

# 用 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）出需要的範圍。 ![](https://saylordotorg.github.io/text_essentials-of-geographic-information-systems/section_11/a33268f6ff028c24152080d0aa3f2aad.jpg) 針對向量資料的地理處理方法。取自 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) 套件指定好字型的解析度，就可以避免這種情形，得到與繪圖指令相符的字體大小。 ![HCmap_231215.212206](https://hackmd.io/_uploads/ryQySRFLT.png) ### 加上陰影如果在地形圖上加入山體陰影（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`，就可以加入山體陰影，生成視覺上具有立體感的分層設色地形圖了。 ![HCmap_231215.214102](https://hackmd.io/_uploads/HypATCFL6.jpg) 其實從繪圖指令中，可以發現這些地圖的生成，就是把不同的地圖圖層調整（或不調整）透明度後依序疊在一起，最後加上圖例、指北針、比例尺等第圖要素，最後加上標題就完成了。地圖依據需要呈現或強調的資料不同，最後透過不同的製圖設計（[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>

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.