技術標籤：r語言 merge

之前做過兩個空間資料的疊加：

從那之後就想做一個系列，主要想研究sf的工具能不能替代我在ArcGIS裡頭常用的工具。本次先嚐試做一些簡單工作，看看可行性。

這裡先宣告一下，本文中出現的warnings並沒有被解決，一些輸出資料可能不精確。

1. 讀入shp檔案並加入FID

程式碼：

grid<-read_sf('BJ_1km_132x138.shp', fid_column_name = 'FID')
grid_my<-read_sf('GRID_WGS84_ID.shp')
 

說明：

fid_column_name為字串則會在資料框裡頭新增一列以該字串命名的FID，shp檔案的FID會從0開始編號。如果要把處理資料跟原始檔案匹配的話，建議加上FID。

2. 實現ArcGIS中的join功能

程式碼：

  raw<-fread(filepath, 
             colClasses = c('character',rep('numeric',2)),
             col.names = c('grid','hour','value'))
  
  grid_my2<-merge(grid_my, raw, by.x = c('ID'), by.y = c('grid'), all.y = T)
  
 

說明：

可以使用merge函式來進行分類匹配，與一般資料的merge操作相同。

3. 計算面要素中心點

程式碼：

grid_my2<-grid_my2 %>% st_centroid()

說明：

允許管道操作符。st_centroid()計算之後的sf物件裡，geometry會從面要素變成點要素，顯示的是中心點座標。

這裡報了個warning是沒有計算精確的經緯度，但我看結果應該是精確的，不清楚怎麼回事，姑且把warning貼出來：

>   grid_my2<-grid_my2 %>% st_centroid()
Warning messages:
1: In st_centroid.sf(.) :
  st_centroid assumes attributes are constant over geometries of x
2: In st_centroid.sfc(st_geometry(x), of_largest_polygon = of_largest_polygon) :
  st_centroid does not give correct centroids for longitude/latitude data
 

4. Dissolve

程式碼：

st_union(grid_my2, by_feature = F)

說明：

應該是跟ArcGIS中Dissolve的功能一樣，但似乎我原始資料有些問題，結果不太好，因此說明僅供參考。若by_feature為TRUE，則會按除geometry之外的列分組dissolve。

5. 轉換座標系

程式碼：

grid_wgs84<-st_transform(grid, crs = 4326)
st_crs(grid_cmaq_wgs84)

說明：

傳座標系迴轉換後的sf物件，crs是EPSG code，4326為WGS84座標系。st_crs()顯示sf物件的座標系資訊。

6. intersect

上一篇文章已經示範過了，見文章頂部連結。這裡要補充的是intersect的兩個物件的座標系要相同。這裡同樣也報了個warning：

> test<-st_intersection(grid,grid_my2)
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning message:
attribute variables are assumed to be spatially constant throughout all geometries

7. 轉換成其他物件

程式碼：

test<-data.table(test)[,geometry:=NULL][,ID_12_13:=NULL]
test<-test[,lapply(.SD, sum), by = c('FID')]
test<-merge(grid, test, by = 'FID', all.x = T)
test<-select(test,FID,variable,value,geometry)

說明：

sf可以轉換為其他物件做資料處理。FID的功能在這裡體現，資料處理後再與sf物件merge。

8. 輸出shp檔案

程式碼：

st_write(test, 'grid_new.shp')

說明：跟一般輸出資料的方法相同。

以上是幾個ArcGIS軟體功能在sf包裡頭的實現方法，其他關於sf包的內容還是推薦看官方文件：

這個包對我的意義在於一些重複性的空間資料操作能通過R程式迴圈實現，厭煩用滑鼠來回做重複性操作的同學們可以試試。