用 R 建構網格系統

環境設定：安裝/ 讀入套件

建構網格系統

設定網格系統的範圍
產生網格系統
Select by location
Spatial join
輸出 shapefile

對有系統性的調查或監測而言，網格系統是選擇樣區時很常用到的工具。舉凡BBS、路殺社系統化調查等。除此之外，在時空資料的視覺化，網格系統也是很常用的輔助工具。在 QGIS 跟 ArcGIS 內都有相對應的工具來產生網格系統 (QGIS 內叫 Vector grid，ArcGIS 應該是 Fishnet)，不過當要產生的網格數量比較多的時候 (e.g. 台灣含東沙南沙的 1 x 1 公里網格)，用 R 來操作應該可以節省不少時間，至少把流程建在 code 中，當要重複做類似的事情時，不用一直在那邊點左鍵XD

以下範例為建置台灣範圍 (含離島) 坐標系統為 TWD97 121 的網格系統。為了讓最後呈現的圖示比較清楚，我們這邊建立的網格解析度設為 10 x 10 公里。
(註：本篇只探討方法，至於符不符合 TM2 這件事，就先不在這裡討論啦…:P)

環境設定：安裝/ 讀入套件

# 安裝套件
list.of.packages <- c("sf", "rgdal", "magrittr", "maptools")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

# 讀入套件
library(sf)       # for st_make_grid
library(rgdal)    # for readOGR, writeOGR
library(magrittr) # for pipe

建構網格系統

這篇示範的方法是以 sf 這個套件的函數為主，在這套建構網格系統方法中，有幾個必要的參數需要設定，包括網格系統的涵蓋範圍 (grid extent)、坐標系統 (coordinate system) 以及網格系統的解析度 (resolution)，解析度要換算成欄跟列的網格數量。

設定網格系統的範圍

網格系統的範圍 (grid extent) 是用 5 x 2 的矩陣來指定，順序依序為左上→右上→右下→左下→左上，我們先透過台灣縣市界縣的範圍來確定網格系統的範圍。

讀取台灣鄉鎮市區界圖層(shapefile)

TW.island <- readOGR("E:/for R markdown/mapdata", # 檔案路徑
                     "TOWN_MOI_1070516", # 檔案名稱
                     use_iconv = TRUE, # 允許設定encoding
                     encoding = "UTF-8") %>% # encoding設定為UTF-8
  # 將圖層坐標系統轉成與欲設立的網格座標系統相符 (TWD97_121, epsg代碼為3826)
  spTransform(CRS("+init=epsg:3826"))

確定shapefile圖層的範圍 (bbox)

TW.island@bbox

##          min       max
## x  -40418.36  606875.6
## y 2422004.77 2919551.3

設定網格系統範圍，範圍參照上面的數值往外延伸到以萬 ( 10 公里) 為單位

grid.extent <- 
  matrix(c(-50000,  2920000,   # 左上 (Xmin, Ymax)
           610000,  2920000,   # 右上 (Xmax, Ymax)
           610000, 2420000,    # 右下 (Xmax, Ymin)
           -50000, 2420000,    # 左下 (Xmin, Ymin)
           -50000,  2920000),  # 左上 (Xmin, Ymax) 
         byrow = TRUE, ncol = 2) %>%
  # 轉成list以符合st_polygon的需求
  list() %>% 
  # 產生網格系統範圍的polygon
  st_polygon() %>% 
  # 轉格式跟設置坐標系統
  st_sfc(., crs = 3826)

產生網格系統

這邊使用的是 sf package 的 st_make_grid()，製作出來的網格系統格式為 sfc，所以還要透過 st_sf() 將網格轉成帶有屬性表的 sf 格式。因為後續在網格選取以及 spatial join 是透過 rgdal 的函數，所以要再用 as() 將網格資料轉成 SpatialPolygonDataFrame。

Grid.sys <- 
  # 製作網格系統
  st_make_grid(grid.extent,  # 網格系統的範圍：grid.extent
               n = c(66, 50),   # 網格系統的解析度：66欄 x 50列 
               crs = 3826,  # 坐標系統：TWD97 121
               what = 'polygons') %>%  # 產出形式：polygon
  # 轉成sf格式，並且創造ID欄位
  st_sf('geometry' = ., data.frame('ID' = 1:length(.))) %>%
  # 轉成sp格式
  as(., 'Spatial') %>%
  # 再指定一次座標系統，確保與台灣範圍的圖層一致
  spTransform(CRS("+init=epsg:3826"))

Select by location

選取與台灣陸域範圍重疊的網格

Grid.TW <- 
  Grid.sys[subset(TW.island),] # 逗號不能省略

Spatial join

將台灣範圍圖層的鄉鎮資訊綁定到網格系統的屬性表中

TW_info <- 
  over(Grid.TW, TW.island)  # creat a data.frame of IDs in IBA for 1km grid

Grid.TW@data %<>% cbind(., TW_info)

新創建的網格系統範圍

網格系統屬性表的結構

	ID	TOWNID	TOWNCODE	COUNTYNAME	TOWNNAME	TOWNENG	COUNTYID	COUNTYCODE
28	28	T04	10013040	屏東縣	恆春鎮	Hengchun Township	T	10013
29	29	T04	10013040	屏東縣	恆春鎮	Hengchun Township	T	10013
37	37	V16	10014160	臺東縣	蘭嶼鄉	Lanyu Township	V	10014
93	93	T04	10013040	屏東縣	恆春鎮	Hengchun Township	T	10013
94	94	T04	10013040	屏東縣	恆春鎮	Hengchun Township	T	10013
95	95	T04	10013040	屏東縣	恆春鎮	Hengchun Township	T	10013
102	102	V16	10014160	臺東縣	蘭嶼鄉	Lanyu Township	V	10014
103	103	V16	10014160	臺東縣	蘭嶼鄉	Lanyu Township	V	10014
159	159	T32	10013320	屏東縣	獅子鄉	Shizi Township	T	10013
160	160	T32	10013320	屏東縣	獅子鄉	Shizi Township	T	10013

輸出 shapefile

writeOGR(Grid_island,
         "E:/for R markdown", # out put directoty
         "Grid_sys_tw", # name of shapefile
         driver = "ESRI Shapefile") #file format