Vector 的操作 (一)
Kemushi54
2018年4月7日
這一篇主要目的是整理用 R 來操作平常在 ArcGIS 或是 QGIS 上操作的基本功能，在 GIS 的概念部分不會著墨太多，所以如果對 GIS 沒有概念的話，建議還是先上一下相關的課程或是爬個估狗了，另外這一篇是針對 vector 的操作，跟 raster 有關的會放在另一篇。
在 R 裡面可以讀取及操作 vector 的 package 有很多，我自己習慣用的是 sp 以及 rgdal (gdal念作“goodle”) 這個系列，其中套件 sp 提供了空間資料的資料結構，classes 與 methods，在 R 裡面常見的資料格式包含SpatialPoint/ SpatialLine/ SpatialPolygon/ SpatialPointDataFrame/ SpatialLineDataFrame/ SpatialPolygonDataFrame。而套件 rgdal 則是將 GDAL (Geospatial Data Abstraction Library) 與 R 綁在一起。另外還有 sf 系列，sf 在空間資料處理上又是另一種結構，以後有機會再來介紹。
本篇整理的功能如下：
- Vector 的讀寫
- 座標系統轉換
- Dissolve
- Buffer
- Select by location
- Spatial join
以下的 vector 的檔案格式都會用相對普遍的 shapefile (.shp) 為代表。
安裝與讀入套件
# 安裝套件
list.of.packages <- c("sp", "rgdal", "rgeos", "raster")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)# 讀入套件
library(rgdal)讀入 rgdal 時 sp 也會一併讀入，所以不用特別去讀。
Vector 的讀寫
在 rgdal 中讀取 shapefile 的函數為 readOGR(dsn, layer, ...)，其中 dsn (data source name) 為 存放 shapefile 的資料夾路徑，不包含檔案名稱；layer 為 shapefile 的檔案名稱，不需要副檔名，這兩個為必要的參數。另外可以用 p4s = 來指定讀入檔案的座標系統，預設是 NULL，會自動讀取 shapefile 的 .prj 檔；用 encoding = 來設定讀入檔案的書寫編碼，記得同時要設定 use_iconv = TRUE 來開啟設定編碼的功能。
這邊使用的範例檔為政府開放資料平臺的臺灣縣市界線向量檔，可以按這裡下載。
# 讀入 shapefile
TW <- readOGR("F:/R markdown/COUNTY", "TW_COUNTY")## OGR data source with driver: ESRI Shapefile
## Source: "F:/R markdown/COUNTY", layer: "TW_COUNTY"
## with 22 features
## It has 4 fields# 簡單畫個圖看一下shapefile長什麼樣
plot(TW)
# 看一下在R裡面是什麼樣子
summary(TW)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 118.1437 124.56115
## y 21.8956 26.38528
## Is projected: FALSE
## proj4string : [+proj=longlat +ellps=GRS80 +no_defs]
## Data attributes:
## COUNTYID COUNTYCODE COUNTYNAME COUNTYENG
## A : 1 09007 : 1 宜蘭縣 : 1 Changhua County: 1
## B : 1 09020 : 1 花蓮縣 : 1 Chiayi City : 1
## C : 1 10002 : 1 金門縣 : 1 Chiayi County : 1
## D : 1 10004 : 1 南投縣 : 1 Hsinchu City : 1
## E : 1 10005 : 1 屏東縣 : 1 Hsinchu County : 1
## F : 1 10007 : 1 苗栗縣 : 1 Hualien County : 1
## (Other):16 (Other):16 (Other):16 (Other) :16用readOGR()讀進來的縣市界線檔資料格式為 SpatialPolygonDataFrame，包含了空間資 polygons (每個 polygon 的空間資訊), bbox (整個 vector 的範圍), proj4string (座標系統) 以及 data，shapefile 的 attribute table。
若要各別讀取上述的資料可以在變數後面加上 @ 來讀取。例如下面我們用@data來擷取 attribute table 的部分。
# 擷取attribute table並且存到TW.attr
TW.attr <- TW@data
# 看一下擷取的資料
head(TW.attr)## COUNTYID COUNTYCODE COUNTYNAME COUNTYENG
## 0 Z 09007 連江縣 Lienchiang County
## 1 W 09020 金門縣 Kinmen County
## 2 G 10002 宜蘭縣 Yilan County
## 3 N 10007 彰化縣 Changhua County
## 4 M 10008 南投縣 Nantou County
## 5 P 10009 雲林縣 Yunlin County
如果要修改 attribute table，可以用 @data 擷取出資料後修改，再指定回 @data 就可以了。像下面這樣~
# 在剛才擷取的TW.attr增加一欄ID.test，填入每列的序號
TW.attr$ID.test <- 1:nrow(TW)
# 將TW.attr指定為TW的attribute table (TW@data)
TW@data <- TW.attr
# 看一下新的attribute table
head(TW@data)## COUNTYID COUNTYCODE COUNTYNAME COUNTYENG ID.test
## 0 Z 09007 連江縣 Lienchiang County 1
## 1 W 09020 金門縣 Kinmen County 2
## 2 G 10002 宜蘭縣 Yilan County 3
## 3 N 10007 彰化縣 Changhua County 4
## 4 M 10008 南投縣 Nantou County 5
## 5 P 10009 雲林縣 Yunlin County 6
另外一種在 GIS 軟體讀 vector 的方式是將有坐標的 csv 檔或 txt 檔，將坐標的部分轉成空間資訊，在 R 裡面也可以做到。概念上是先將資料以 data.frame 或 data.table 的方式讀進 R，再把坐標的部分轉換成空間資訊。這邊拿之前用過的 BBS 資料來作示範。
# 讀入需要的套件
library(data.table)
library(magrittr)
# 讀取BBS資料
BBS <-
fread("F:/R markdown/BBS2009-2015/occurrence.txt",
sep = "\t", encoding = "UTF-8") %>%
# ↓刪除全都是NA值的欄位
Filter(function(x)!all(is.na(x)), .) %>%
# ↓簡單清理資料，刪除坐標為NA的紀錄
.[!is.na(decimalLongitude) & !is.na(decimalLatitude)]
將 data.frame/ data.table 轉成 SpatialPointDataFram，在 BBS 的資料中，x 與 y 的資訊分別是 decimalLongitude 跟 decimalLatitude 這兩個欄位。
coordinates(BBS) <- ~ decimalLongitude + decimalLatitude
# 前面的~不能省略，xy之間用+來連接
proj4string(BBS) <- CRS("+init=epsg:4326")
# 用上面方法轉成的SpatialPointDataFrame沒有坐標系統，需要利用proj4string()來指定，BBS的坐標系統是WGS84，EPSG代碼為4326，在這邊用"+init=epsg:4326"來指定。
# 看一下轉格式後的BBS資料格式
class(BBS)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
如果要將轉換成空間資料形式的 BBS 資料輸出成 shapefile，在 rgdal 內可以用 writeOGR()。 writeOGR(obj, dsn, layer, driver)，obj 是要儲存的資料；dsn 是要存放的資料夾路徑；layer 是檔案名稱，不用副檔名；driver 指定要儲存的檔案格式，要存成 shapefile 的話為 driver = "ESRI Shapefile"。
writeOGR(BBS, "F:/R markdown", "BBS", driver = "ESRI Shapefile")坐標系統轉換
當兩個 vector 要合併或是 spatial join 等操作時，通常會需要一致的坐標系統，上面講到的 proj4string() 是用在當空間資料 (sp 格式) 沒有指定坐標系統的時候的指定方式，但是如果資料本身是有指定坐標系統，就需要做轉換了。這邊用的函數是 spTransform(data, CRS)，data 放 sp 格式的空間資料，CRS 為座標系統，可以用上面 epsg 的方式指定，也可以寫 proj4string 的形式。
把前面的台灣縣市界圖層轉成跟 BBS 一致的 WGS84。
TW %<>% spTransform(., CRS("+init=epsg:4326"))在 R 裡面，即使實際上是一樣的坐標系統，會因為在 proj4string 上的寫法不一樣，在 R 裡面被當作不同的坐標系統，導致進一步的操作 error 不斷。要確保兩個圖層的坐標系統一致的方法，其中一個做法是像上面一樣，都用一樣的方式 (“init=epsg:4326”) 來指定坐標系統，另外一種做法，是以其中一個資料的坐標系統為基準，去指定另一個資料的坐標格式。
以 BBS 的坐標系統為基準，把台灣縣市界的坐標系統轉為跟 BBS 資料一致
TW %<>% spTransform(., CRS(proj4string(BBS)))Dissolve
將在 attribute table 特定欄位具有相同值的多個 polygon 合併為一個。這邊示範的方法是 raster 套件的 aggregate() 函數。aggregate(sp, by = , FUN)，sp 為要 dissolve 的空間資料，用 by = 來指定合併依據的 attribute table 欄位，FUN 則是 dissolve 時特定欄位的計算，例如計算面積加總、保留 dissolve 的 polygon中某個欄位的最大值等。
將所有 TW 的 polygon dissolve 在一起
## Loading required namespace: rgeos
依據 COUNTYID 來 dissolve
TW.county <- aggregate(TW, by = "COUNTYID")
# 不過這邊原本一個COUNTYID下也只有一個polygon，結果沒差就不顯示了...Buffer
以 point/ line/ polygon 為中心，像外延伸一定距離產生新的 polygon。這邊用的是 rgeos 套件的 gBuffer() 函數 (注意 B 要大寫)。 gBuffer(spgeom, byid=FALSE, id=NULL, width=1.0)，spgeom 是我們原本的空間資料，byid = TRUE/FALSE 指定 buffer 是要分各別的 polygon 還是 buffer 後全部 dissolve 起來；id = 指定 buffer 後的 attribute table，預設 NULL 會填入原本 vector 的 attribute table，但如果前面的 byid 是設 FALSE，attribute table 則不會填入任何東西；width 指定要 buffer 的距離，坐標系統的依據是原本 vector 的坐標系統，所以當 vector 坐標系統為 WGS84，width = 1指的是 buffer 寬度為 1 度。如果想 buffer 的距離是公尺或公里的話，建議先將 vector 的座標系統轉成 TWD97 (epsg:3826) 或TWD67 (epsg:3828)。
畫圖看一下結果
# 原本的縣市界線 (黑色)
plot(TW)
# 剛才的buffer (藍色)
plot(TW.buffer, border = "blue", add = TRUE) 
Select by location
選取落在 vector A 裡面的 vector B 的資料，在 R 裡面可以用 subset， vector_A [ subset ( vector_B, condition ) , ]。中括號前面放要篩選的 vector (vector_A)。subset() 裡面放篩選條件的 vector (vector_B)，第一個位置擺 vector，後面放 vector 內 attribute table 的欄位與條件 (不設條件就是選取整個 vector_B 涵蓋到的 vector_A)
選取在台灣縣市界範圍內的BBS資料
BBS.sel <- BBS[subset(TW), ]選取在南投縣內的BBS資料
BBS.Nantou <- BBS[subset(TW, COUNTYNAME == "南投縣"), ]上面第二個例子中是指選 COUNTYNAME 為南投縣的 polygon 內的資料。小括號後面還有一個逗號是不可以省略的！而且很容易忘記加，要特別注意~！
畫圖看一下剛剛的結果
# 台灣縣市界線 (黑線)
plot(TW)
# 剛才選出來的 BBS 點位 (藍圈圈)
points(BBS.Nantou, col = "blue")
Spatial join
這篇最後一個要介紹的功能是 spatial join，在 QGIS 裡面的名字是 Join attributes by location，也就是把同位置 vector_A 的 attribute table 加到 vector_B 上面，例如把所有 BBS 的資料加上縣市界線裡面的縣市名稱資訊。這邊用的是 sp 套件的 over 函是 data.frame。
擷取 TW 的所有 attribute table
TW.attr.BBS <- over(BBS, TW)也可以指定只擷取特定的欄位，例如只擷取 COUNTYNAME
county.BBS <- over(BBS, TW[, "COUNTYNAME"])
# 看一下擷取出的資料
head(county.BBS)## COUNTYNAME
## 1 彰化縣
## 2 彰化縣
## 3 彰化縣
## 4 彰化縣
## 5 彰化縣
## 6 彰化縣# 再把上面用 over 擷取出的資料指定回BBS的attribute table內就可以了
BBS$county <- county.BBS$COUNTYNAME
# 左右兩邊的資料格式要對等，$COUNTYNAME不可以省略!!以上是這一次的介紹，空間資訊處理上的方法太多了，所以考慮分篇來寫。另外跟前面幾篇一樣，再 R 裡面要做到同一件事情，有很多函數可以使用也有很多種寫法，這邊介紹的都只是我在拜了眾多估狗大神後選擇用的寫法，不過因為當初多少還是有比較過不同方法，所以上面提供的這些方法效率應該也不會太差，但如果有更有效率的寫法拜託告訴我！ 礙於篇幅太長大家可能會沒耐心看下去，裡面幾乎每個函數都還有很多參數的設定是這裡沒有提到的，使用前還是建議善用 R 內建的 help，先看過使用說明書再來使用。