Leaflet & 高德｜试试全新的地图绘制思路
Leaflet & 高德｜试试全新的地图绘制思路
前言
统计学上常说：“一图胜千言”，而对于地图类可视化教程，大家介绍的并不是很多，本教程讲会教你如何绘制各类地图。首先介绍下地图的种类（个人经验仅作参考）：
-
填色地图 (Choropleth Maps)：饼图的“饼”被咬了几口，然后再拼在一起就是填色地图了。
-
散点地图 (Scatter on Maps)：散点图的 X 和 Y 轴改成经度和纬度，再使用图片（地图）作为背景。
-
气泡地图 (Bubble Maps)：和散点地图类似，只不过新增加了一个变量来定义点的大小。
-
路径地图 (Lines on Maps)： 线图的 X 轴和 Y 轴改成经度和纬度，再使用图片（地图）作为背景。
-
其他地图 (Other Maps)：如：房地产售楼的规划图，天气预报的云图，NASA 的城市灯光图等。
本文框架
1.材料准备
- 画统计图最重要的是啥？数据！没有数据怎么画？难不成去小卖部买个红蓝铅笔画个徐霞客同款？画地图通常需要两个数据。1. 经纬度数据（类似散点图中的 X 轴与 Y 轴）；2. 地图数据（类似散点图的背景图片）
下面我会一一讲解如何获取这两类数据。
- 得到数据之后，我们还需要相应软件来实现地图的可视化。老牌软件有 ArcGIS 和 PPT 插件等，但是正版费用较高。而 R 是开源软件，绘图也是它的强项，所以此教程我们将采用 R 语言的 leaflet 包进行地图的可视化。
2. 数据获取
2.1 经纬度数据
就拿今天讲的最简单的散点地图来说，需要通过经纬度来确定点的位置，然后再映射到地图上。
- 如果你手头有经纬度数据，恭喜你，你不用去找了，但是一定要确定你的经纬度数据使用的是哪种坐标系？（因为同一地点不同坐标系里的经纬度不一样，偏移大概有几百米，所以千万不要混用）
- WGS84坐标系：一般是谷歌等国外地图使用；
- GCJ02坐标系（加密的火星坐标系）：国内的高德地图和腾讯地图等使用；
- BD-09 坐标系（再次加密的火星坐标系）：国内的百度地图使用；
因为本教程为了适用性使用的是高德的底图（GCJ02坐标系），如果您是WGS84坐标系在后续代码中删除高德的底图就好（一定会面临主权问题）；如果您是 BD-09 坐标系，这个需要转换且比较复杂。
- 如果你手头没有数据，只有详细的地址，需要转换成经纬度。
简单的是自己去搜一下“坐标拾取”，然后借助百度地图等把地址转换成坐标，但是一次就只能拾取一个，并且频繁拾取还要验证码，如果地址比较多就太麻烦了。这时候可以使用高德提供的api进行批量查询地址对应的坐标。
2.1.1 批量查询经纬度（通过高德地图）
- 打开高德开放平台，并注册认证为个人开发者。
- 点击控制台，登录你的账户，打开左侧“应用管理——我的应用——创建新应用”，随便输入名称等信息创建应用，创建成功后点击右侧的“添加”，输入自定义名称，并设置服务平台为 Web 服务，IP 白名单有需求可以设置，没有就算了。
- 点击提交，此时就可以看到你项目的 key 了（下图红圈），复制下。
- 根据高德提供的地理编码开发文档，最终 url 应该是：
address <- '福建省厦门市厦门北站' ## 详细地址
city <- '厦门'
key <- 'Your key' ## 你刚刚复制的key
url <- paste0(
'https://restapi.amap.com/v3/geocode/geo?',
'&key=', key, ## 你刚刚复制的key
'&address=', address, ## 详细地址
'&output=', 'JSON',
'&city=', city ## 这个可以不要，填所在城市就好
)
- 使用 jsonlite 解析 api 地址
library(jsonlite)
temp_geo <- fromJSON(paste(readLines(url,warn = F, encoding = 'UTF-8'), collapse = ""))
temp_geo
## $status
## [1] "1"
##
## $info
## [1] "OK"
##
## $infocode
## [1] "10000"
##
## $count
## [1] "1"
##
## $geocodes
## formatted_address country province citycode city district township
## 1 福建省厦门市集美区厦门北站 中国 福建省 0592 厦门市 集美区 NULL
## neighborhood.name neighborhood.type building.name building.type adcode street
## 1 NULL NULL NULL NULL 350211 NULL
## number location level
## 1 NULL 118.074013,24.636101 兴趣点
- 分离经纬度
高德地图返回的经纬度信息放在了 geocodes 下面的 location，其他的返回参数可以参考高德地图提供的开发文档。
library(tidyr)
library(dplyr)
datafile <- data.frame( ## 提取经纬度
address = address,
city = temp_geo$geocodes$city,
geo = temp_geo$geocodes$location
) %>%
separate(col = geo, into = c('lng', 'lat'), sep = ',') ## 分割经度和纬度
datafile$lng <- as.numeric(datafile$lng)
datafile$lat <- as.numeric(datafile$lat)
head(datafile)
## address city lng lat
## 1 福建省厦门市厦门北站 厦门市 118.074 24.6361
2.2 地图数据
用过 ArcGIS 朋友会比较清楚，常用的格式有 shp，Geojson 等，这些大家可以去 Github 或者万能的淘宝找。如何处理和使用这些数据？我们会另作一期推文。
注意： 在收集的时候一定要注意主权完整，台湾省和南疆部分是中国领土，南海九段线是中国领海！此外，本教程不包括填色地图，所以不需要自定义地图数据。为了方便起见，本教程直接使用高德地图提供的底图。
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 来查询导航路径，详细可见：官方说明文档
- 通过地理编码获得起点和终点的经纬度，方法同上不再赘述了。
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)
)
小编有话说
受限于微信平台问题，所有的leaflet画的图都是以截图方式呈现，但是实际上leaflet生成的是交互式的地图，也就是你可以像导航软件里一样放大和缩小地图，不用受限于分辨率问题；上述代码中设置的label都是鼠标悬停显示。
另外，很多教程都没有提一件事：
根据《中华人民共和国测绘法》等有关法律、法规规定：在中华人民共和国境内公开出版地图、引进地图、展示、登载地图以及在生产加工的产品上附加的地图图形都需要经审核，审核通过之后编发审图号。
本教程使用的是高德的底图，所以可以直接使用高德提供的审图号。如果是来历不明的地图数据，无法提供审图号可能会引来一些不必要的麻烦。