R空間數據處理與可視化

前言

不少朋友說在R裏無法使用高德地圖，這裏給出一個基於leaflet包的解決方法。

library(leaflet)

# 添加高德地圖
m <- leaflet() %>%
      addTiles(
        'http://webrd0{s}.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',
        options=tileOptions(tileSize=256, minZoom=9, maxZoom=17, subdomains="1234"),
        attribution = '&copy; <a href="http://ditu.amap.com/">高德地圖</a>',
        group="Road Map"
      ) %>% 
      setView(116.40,39.90, zoom = 10)
m

固然，除了告訴你們怎麼在R裏調用高德地圖外，本文還想作的深刻一些，嘗試對空間可視化的基礎知識作一個簡單的介紹。

空間數據基礎知識

shp文件

空間數據最經常使用的格式是shp，主要由三個文件組成：shp文件用於存儲位置幾何信息，dbf文件用於存儲attribute，shx用於存儲位置幾何信息與attribute的對照表。位置幾何信息主要有如下幾類：points，multipoints，lines，polygons等。

WKT與WKB

WKT(Well-known text)是開放地理空間聯盟OGC（Open GIS Consortium ）制定的一種文本標記語言，用於表示矢量幾何對象、空間參照系統及空間參照系統之間的轉換。舉例以下：

• 點（Point）："POINT(1 1)"

• 線（Line）："LINESTRING(0 0,1 1,2 2)"

• 多邊形（Polygon）："POLYGON((0 0,3 0,3 3,0 3,0 0),(1 1,2 1,2 2,1 2,1 1))"

WKB(well-known binary) 是WKT的二進制表示形式，解決了WKT表達方式冗餘的問題，便於傳輸和在數據庫中存儲相同的信息.

R的空間數據處理與可視化

空間數據處理與可視化，須要解決三個問題，一是怎麼在R中表示空間數據，二是怎麼對空間對象進行計算；三是怎麼在R中繪製空間數據/地圖。sp用於解決第一個問題，rgeos用於解決第二個問題，leaflet用於解決第三個問題。

sp

sp包的功能是在R中提供對象表示shp文件。SpatialPoints，SpatialMultiPoints，SpatialLines，SpatialPolygons等用於表示位置幾何信息。attribute通常以表格形式存在，因此sp包用dataframe對齊進行表示。爲前面提到的SpatialXXX添加dataframe後獲得諸如SpatialPointsDataFrame，SpatialMutilPointsDataFrame，SpatialLinesDataFrame，SpatialPolygonsDataFrame等類。在這些類中，位置幾何信息與attribute的對照關係經過Spatial類的ID與dataframe的rownames進行匹配獲得。

SpatialXXDataFrame的結構示意圖以下（出處：http://neondataskills.org/R/）：

下面舉一個例子，怎麼從dataframe數據變爲sp對象。

library(splitstackshape)
library(sp)
library(dplyr)
library(tidyr)

# 準備測試數據
link_id <- c("road_one", "road_two")  # 兩條道路，道路1和道路2
coors <- c("116.44469451904297,39.890071868896484:116.44451141357422,39.891361236572266", "116.44499969482422,39.887630462646484:116.44469451904297,39.890071868896484")  # 道路1的經緯度座標序列和道路2的經緯度座標序列
status <- c("congest", "uncongest")  # 道路1處於擁堵狀態，道路2處於暢通狀態

link_coors <- data.frame(link_id, coors, status)
lon_lat_df <- cSplit(link_coors %>% select(link_id, coors), 
       c("coors"), sep=":", direction="long") %>% 
  separate(coors, c("lng", "lat"),  sep=",", convert=TRUE)

# 轉化函數
df2sp <- function(route_df) {
  # 將df的一行轉化爲一個Lines
  xy2sp <- function(route_df) {
    coors <- route_df$coors
    link_id <- route_df$link_id
    line <- coors %>%
      stringr::str_split(pattern=":", simplify=T) %>%
      t() %>%
      stringr::str_split(pattern=",", simplify=T) %>%
      apply(2, as.numeric) %>%
      Line() %>%
      list() %>%
      Lines(ID=link_id)
    return(line)  
  }
  
  # 幾何信息join屬性信息
  splines2splinesdf <- function(splines,  data, id_field)  {
    ids <- data.frame(names(splines), stringsAsFactors =F)
    colnames(ids) <- id_field
    join_name <- dplyr::inner_join(ids, data)
    row.names(join_name ) <- join_name[, id_field]
    splinesdf <- SpatialLinesDataFrame(splines, data=join_name)
    proj4string(splinesdf ) <- CRS("+init=epsg:4326") # 設置投影座標系，leaflet能夠不用設置
    return(splinesdf)
  }
  
  route_list <- plyr::alply(route_df, 1, xy2sp)
  attributes(route_list) <- NULL  # 必須設置，不然leaflet不可識別
  spline <- SpatialLines(route_list)
  return(splines2splinesdf(spline, route_df, "link_id"))
}

# data.frame轉化爲sp
Sldf <- df2sp(link_coors)

str(Sldf)

rgeos

空間處理，主要用來作一些空間運算，好比計算兩個空間對象的位置關係：相交，重疊，包含等等。再好比，根據空間對象建立新的空間對象。此外，rgeos還可以完成WKT與sp對象的相互轉換。

library(rgeos)

# 建立外擴與內縮buffer，演示WKT的讀寫
dilated_buffer <- gBuffer(Sldf, byid=TRUE, width=0.0002, capStyle="FLAT")
dilated_buffer_wkt <- readWKT(writeWKT(dilated_buffer, byid = FALSE))
eroded_buffer <- gBuffer(dilated_buffer, byid=TRUE, width=-0.0001, capStyle="SQUARE")

leaflet

咱們繼續上面的例子，將空間對象繪製到高德地圖上。

library(leaflet)
factpal <- colorFactor(c(rgb(1,0,0,1),rgb(0, 1, 0, 1)), domain=c("congest", "uncongest"))

m <- leaflet() %>%
  addTiles(
    'http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',
    options = tileOptions(tileSize=256, minZoom=9, maxZoom=17),
    attribution = '&copy; <a href="http://ditu.amap.com/">高德地圖</a>',
    group="高德地圖"
  ) %>%
  setView(116.40,39.90, zoom = 10)  %>%
  addPolylines(color=~factpal(status), weight=3,opacity=1,  data=Sldf, group="實時路況") %>%
  addPolygons(data=dilated_buffer_wkt, group="空間計算") %>% 
  addPolygons(data=eroded_buffer, color="black", group="空間計算") %>%
  addLayersControl(
    overlayGroups = c("高德地圖", "實時路況", "空間計算"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  addLegend("bottomleft", pal = factpal, values = Sldf@data$status,
            title = "實時交通", 
            opacity = 1
  ) %>%
  fitBounds(Sldf@bbox["x", "min"] - 0.001, Sldf@bbox["y", "min"] - 0.001, 
            Sldf@bbox["x", "max"] + 0.001, Sldf@bbox["y", "max"] + 0.001
  )

m

後記

R空間處理的第三方庫主要由sp提供了R的存儲結構，rgdal提供了讀寫操做，rgeos提供了運算操做。然而，sp的實現只實現了空間數據處理標準Simple Feature的一個子集，且處理效率較低。爲此，一個全新的sf包正在開發中，目標是替換掉sp，並把rgdal和rgeos的功能整合進來。目前已經具有了基本的使用功能，咱們可使用sf來完成從數據框建立空間對象的操做，能夠看到代碼簡單了不少：

link_wkt <- link_coors %>%
  mutate(wkt_prefix="LINESTRING(",
         wkt_content=str_replace_all(coors, ",", " ") %>% str_replace_all(":", ", "),
         wkt_posix=")",
         geom=str_c(wkt_prefix, wkt_content, wkt_posix)
         ) %>%
  select(link_id, status, geom)

link_sf <- st_as_sf(link_wkt, wkt="geom")

# sf提供的buffer運算功能還不夠完善，咱們仍須要使用rgeos來完成相應的操做；rgeos只認sp，因此要作一次從sf到sp的轉化。
Sldf <- as(link_sf, "Spatial")