前言
我们对 leaflet 包做了一期简单的入门教程：空间地理数据可视化之 leaflet 包及其拓展。之后，又将高德和该包相结合，介绍了前期需要准备的工作，见：Leaflet 与高德合并会擦出怎么样的火花？。这一期就到了绘制地图环节，下面将分享三类数据的绘制教程。
本文框架
这篇推文讲的是第三节内容
3. 绘制地图
3.1 散点地图绘制
- 将高德地图替换 leaflet 自带的底图
由于 leaflet 自带的底图不是很合规，所以我们使用高德地图进行替换。
library(leaflet) geo_map <- leaflet(width = '100%') %>% addTiles( urlTemplate = 'http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}', layerId = tileOptions(minZoom = 3, maxZoom = 17), attribution = NULL )
- 在地图上增加标记
使用addMarkers()在geo_map上增加散点图。
geo_map %>% addMarkers( data = datafile, ## 存放有坐标的文件 lat = ~lat, ## 变量纬度(latitude)所在列 lng = ~lng, ## 变量经度(longitude)所在列 label = ~address ## 变量标签所在列(这里设置了详细的地址作为标签) )
- 批量下载 api 地址
高德地图地理编码的 api 每天可以查询 30 万次，所以可以将上述查询过程写成一个 function 来进行批量查询。代码如下，代码中的注释已做详细解释。
## 加载packages library(dplyr) library(tidyr) library(leaflet) ## 设置文件 datafile <- data.frame( id = 1:9, address = c( "北京市东城区天安门广场人民英雄纪念碑", "上海市浦东新区东方明珠广播电视塔", "台湾省台北市南京东路台北体育馆", "海南省海口市海口站(出站口)", "香港特别行政区湾仔区岛湾仔博览道1号会议展览中心香港回归祖国纪念碑", "澳门特别行政区花王堂区澳门博物馆大炮台(北入口)", "西藏拉萨市城关区布达拉宫广场西藏和平解放纪念碑", "新疆克拉玛依市平南八路克拉玛依北站", "黑龙江省哈尔滨市道里区斯大林街20号哈尔滨市人民防洪胜利纪念塔" ) ) ## 自定义function gaodemap <- function(x){ address <- stringr::str_replace_all(datafile$address[x], "[^[:alnum:]]", "") ## 去除特殊字符，部分特殊字符会导致查询失败，所以在最开始直接去除 url <- paste0( 'https://restapi.amap.com/v3/geocode/geo?', '&key=', key, ## key '&address=', address, ## 详细地址 '&output=', 'JSON' ) temp_geo <- fromJSON(paste(readLines(url,warn = F, encoding = 'UTF-8'), collapse = "")) status <- temp_geo$status ## 判断查询是否成功 if (status == 1 & length(temp_geo$geocodes) != 0){ ## 有时候成功了但是没有收到数据，所以用了一个“且”进行判定 temp <- paste(datafile$address[x], temp_geo$geocodes$location, sep = ",") } else { temp <- paste(datafile$address[x], 'NULL,NULL', sep = ",") } return(temp) } ## 地理编码 df <- data.frame( geo = sapply(datafile$id, gaodemap, simplify = T) ) %>% separate(col = geo, into = c('address', 'lng', 'lat'), sep = ',') ## 分割function的结果 df$lng <- as.numeric(df$lng) df$lat <- as.numeric(df$lat) df ## address lng ## 1 北京市东城区天安门广场人民英雄纪念碑 116.39768 ## 2 上海市浦东新区东方明珠广播电视塔 121.50003 ## 3 台湾省台北市南京东路台北体育馆 121.54067 ## 4 海南省海口市海口站(出站口) 110.16203 ## 5 香港特别行政区湾仔区岛湾仔博览道1号会议展览中心香港回归祖国纪念碑 114.17738 ## 6 澳门特别行政区花王堂区澳门博物馆大炮台(北入口) 113.54620 ## 7 西藏拉萨市城关区布达拉宫广场西藏和平解放纪念碑 91.11885 ## 8 新疆克拉玛依市平南八路克拉玛依北站 85.08215 ## 9 黑龙江省哈尔滨市道里区斯大林街20号哈尔滨市人民防洪胜利纪念塔 126.61724 ## lat ## 1 39.90462 ## 2 31.23925 ## 3 25.05213 ## 4 20.02724 ## 5 22.28094 ## 6 22.19486 ## 7 29.65127 ## 8 45.56794 ## 9 45.78065
- 绘图
这里我们将上面获得的df数据进行绘制，可以得到以下图形
geo_map <- leaflet(width = '100%') %>% addTiles( urlTemplate = 'http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}', layerId = tileOptions(minZoom = 3, maxZoom = 17), attribution = NULL ) geo_map %>% addMarkers( data = df, lat = ~lat, lng = ~lng, label = ~address )
补充：如果数据量较大，可以考虑使用 parallel 和 foreach 等来实现并行访问和解析 api，因为 R 语言默认是单核运行的，所以会出现**“一核有难，多核围观”**的情形，使用并行运算可以使电脑发挥出多核的优势，提升数据处理速度。
您可能会发现高德限制每秒 api 访问量是 200 次，多核并行会超限，就我的经验而言每次访问和解析大概需要 0.1 秒，16 线程并行查询 api，一秒钟也就160次，更何况还存在网络波动，所以基本上不用担心超限问题。
3.2 路径地图绘制
- 按照画线图的经验，平面直角坐标系中的一条线的位置由两个点决定，而两个点位置由它们分别的坐标 (X, Y) 决定，同理路径地图上的线由起点和终点决定，起点和终点由它们对应的经纬度决定，这样子我们就可以知道绘制路径地图的数据至少需要 4 个值，分别表示起点的经纬度和终点的经纬度。
- 知道了绘图需要的基本数据，后面的就简单了，只要分别查询两个点的经纬度把他们合并到一个表就好了，这里就不赘述了。有数据的朋友也可以直接把数据整理下就行，下面的例子使用上面绘制散点地图的数据。
- 先把数据整理下，假如说希望画从北京到另外 8 个点的直线，我只要在上表中新增一个点作为线段起点就好
df_line <- df %>% mutate( lat_start = lat[1], ## 新增起点纬度 lng_start = lng[1] ## 新增起点经度 ) df_line ## address lng ## 1 北京市东城区天安门广场人民英雄纪念碑 116.39768 ## 2 上海市浦东新区东方明珠广播电视塔 121.50003 ## 3 台湾省台北市南京东路台北体育馆 121.54067 ## 4 海南省海口市海口站(出站口) 110.16203 ## 5 香港特别行政区湾仔区岛湾仔博览道1号会议展览中心香港回归祖国纪念碑 114.17738 ## 6 澳门特别行政区花王堂区澳门博物馆大炮台(北入口) 113.54620 ## 7 西藏拉萨市城关区布达拉宫广场西藏和平解放纪念碑 91.11885 ## 8 新疆克拉玛依市平南八路克拉玛依北站 85.08215 ## 9 黑龙江省哈尔滨市道里区斯大林街20号哈尔滨市人民防洪胜利纪念塔 126.61724 ## lat lat_start lng_start ## 1 39.90462 39.90462 116.3977 ## 2 31.23925 39.90462 116.3977 ## 3 25.05213 39.90462 116.3977 ## 4 20.02724 39.90462 116.3977 ## 5 22.28094 39.90462 116.3977 ## 6 22.19486 39.90462 116.3977 ## 7 29.65127 39.90462 116.3977 ## 8 45.56794 39.90462 116.3977 ## 9 45.78065 39.90462 116.3977
- 分组合并下数据
df_line <- data.frame( group = rownames(df_line), lat = c(df_line$lat, df_line$lat_start), lng = c(df_line$lng, df_line$lng_start) ) %>% group_by(group) %>% arrange(.by_group = TRUE) df_line ## # A tibble: 18 x 3 ## # Groups: group [9] ## group lat lng ## <chr> <dbl> <dbl> ## 1 1 39.9 116. ## 2 1 39.9 116. ## 3 2 31.2 122. ## 4 2 39.9 116. ## 5 3 25.1 122. ## 6 3 39.9 116. ## 7 4 20.0 110. ## 8 4 39.9 116. ## 9 5 22.3 114. ## 10 5 39.9 116. ## 11 6 22.2 114. ## 12 6 39.9 116. ## 13 7 29.7 91.1 ## 14 7 39.9 116. ## 15 8 45.6 85.1 ## 16 8 39.9 116. ## 17 9 45.8 127. ## 18 9 39.9 116.
- 绘图
geo_map <- leaflet(width = '100%') %>% addTiles( urlTemplate = 'http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}', layerId = tileOptions(minZoom = 3, maxZoom = 17), attribution = NULL ) geo_map %>% addPolylines( data = df_line, ## 存放有坐标文件 lng = ~lng, ## 变量经度(longitude)所在列 lat = ~lat, ## 变量纬度(latitude)所在列 group = ~group, ## 线条分组 weight = 1, ## 线条粗细 color = "red" ##线条颜色 )
3.3 导航路径图
这个是看高德地图的时候无意中发现的，高德还提供 api 来查询导航路径，详细可见：官方说明文档[1]
- 通过地理编码获得起点和终点的经纬度，方法同上不再赘述了。
address1 <- '福建省厦门市厦门高崎国际机场' ## 起点详细地址 address2 <- '福建省厦门市厦门大学思明校区' ## 终点详细地址 url1 <- paste0( ## 查询起点的url "https://restapi.amap.com/v3/geocode/geo?address=", address1, "&output=JSON&key=", key ) url2 <- paste0( ## 查询终点的url "https://restapi.amap.com/v3/geocode/geo?address=", address2, "&output=JSON&key=", key ) address1 <- fromJSON(paste(readLines(url1,warn = F, encoding = 'UTF-8'), collapse = "")) ## 起点坐标 address2 <- fromJSON(paste(readLines(url2,warn = F, encoding = 'UTF-8'), collapse = "")) ## 终点坐标 print(address1$geocodes) ## formatted_address country province citycode city district ## 1 福建省厦门市湖里区高崎国际机场 中国 福建省 0592 厦门市 湖里区 ## township neighborhood.name neighborhood.type building.name building.type ## 1 NULL NULL NULL NULL NULL ## adcode street number location level ## 1 350206 NULL NULL 118.137053,24.539473 兴趣点 print(address2$geocodes) ## formatted_address country province citycode city district ## 1 福建省厦门市思明区厦门大学思明校区 中国 福建省 0592 厦门市 思明区 ## township neighborhood.name neighborhood.type building.name building.type ## 1 NULL NULL NULL NULL NULL ## adcode street number location level ## 1 350203 NULL NULL 118.093470,24.437154 兴趣点
- 按照说明文档要求，生成url。
url <- paste0( "https://restapi.amap.com/v5/direction/bicycling?", ## 使用骑行路线的api，如果是步行或者驾车，需要按照说明文档修改即可 "&output=", "JSON", "&key=", key, '&origin=', address1$geocodes$location, ## 设置起点经纬度 '&destination=', address2$geocodes$location, ## 设置终点经纬度 '&show_fields=', 'polyline' ## 设置显示坐标点 ) ## 查询api address <- fromJSON(paste(readLines(url,warn = F, encoding = 'UTF-8'), collapse = ""))
- 导航路径的经纬度信息分成不同段的储存在
“route-paths-steps-polyline”中。
head(address$route$paths$steps[[1]]$polyline) ## [1] "118.13766;....
- 整理数据
## 经纬度分割开 datafile <- lapply(1:length(address[["route"]][["paths"]][["steps"]][[1]][["polyline"]]), function(x){ df <- strsplit(address[["route"]][["paths"]][["steps"]][[1]][["polyline"]][[x]], ";") return(df) }) ## 转化成表格 df <- data.frame( id = matrix(unlist(datafile)) ) %>% separate(id, c('lng', 'lat'), ",") %>% mutate_if(is.character, as.numeric)
- 绘制地图
leaflet(df, width = '100%') %>% addTiles( urlTemplate = 'http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}', layerId = tileOptions(minZoom = 3, maxZoom = 17), attribution = NULL ) %>% addPolylines( ## 绘制路线 lat = ~lat, lng = ~lng ) %>% addMarkers( ## 绘制起点所在点 lat = ~lat[1], lng = ~lng[1], label = paste0('起点:', address1$geocodes$formatted_address) )%>% addMarkers( ## 绘制终点所在点 lat = ~lat[nrow(df)], lng = ~lng[nrow(df)], label = paste0('终点:', address2$geocodes$formatted_address) )
