前言
很多朋友說在R里沒法使用高德地圖，這里給出一個基於leaflet包的解決方法。
library(leaflet)
# 添加高德地圖
m <- leaflet() %>%
addTiles(
'http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',
tileOptions(tileSize=256, minZoom=9, maxZoom=17),
attribution = '© <a href="http://ditu.amap.com/">高德地圖</a>',
) %>%
setView(116.40,39.90, zoom = 10)
m
空間數據基礎知識
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)
# 自定義函數
points2spline <- function(df, id_field, lng_field, lat_field){
df <- as.data.frame(df)
data <- as.matrix(df[,c(lng_field, lat_field)])
id = df[1, id_field]
Lines(list(Line(data)), ID=id)
}
splines2splinesdf <- function(splines, data, id_field) {
ids <- data.frame(names(splines))
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)
}
# 准備數據
link_id <- c("road_one", "road_two")
coors <- c("116.44469451904297,39.890071868896484:116.44451141357422,39.891361236572266", "116.44499969482422,39.887630462646484:116.44469451904297,39.890071868896484")
status <- c("congest", "uncongest")
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)
# data.frame轉化為sp
link_list <- split(lon_lat_df, lon_lat_df$link_id)
names(link_list) <- NULL
Sl <- SpatialLines(plyr::llply(link_list, points2spline, "link_id", "lng", "lat"))
Sldf <- splines2splinesdf(Sl, link_coors, "link_id")
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}',
tileOptions(tileSize=256, minZoom=9, maxZoom=17),
attribution = '© <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
)
本篇博客系轉載，原文請見： https://segmentfault.com/a/1190000006023703