前言

很多朋友说在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,width=-0.0001,capStyle="SQUARE")

leaflet

我们继续上面的例子，将空间对象绘制到高德地图上。

library(leaflet)
factpal <- colorFactor(c(rgb(1,1),rgb(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,maxZoom=17),group="高德地图"
  ) %>%
  setView(116.40,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",Sldf@bBox["x","max"] + 0.001,"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")

R空间数据处理与可视化

前言