百度地图

杜亚磊开发的baidumap包可以调用百度地图服务获取地理可视化的背景底图。百度是国内比较早的互联网公司，之前介绍的 Apache Echarts 也是百度可视化团队领头开发的，百度地理信息可视化平台使用的地理可视化组件来自开源的mapv。

在百度地图上，根据地理名称获得经纬度信息，获得经纬度后，可以将位置标记在地图上，下面以「中国矿业大学（北京）」为例。

remotes::install_github('badbye/baidumap')

BaiduMap_API_Key 是从百度 LBS 云服务申请到的地图 API 访问令牌，保存到本地 .Renviron 文件里 BaiduMap_API_Key=[key]，这样启动 R 软件时就会自动载入。

library(baidumap)
options(baidumap.key = Sys.getenv("BaiduMap_API_Key"))
getCoordinate('中国矿业大学（北京）学院路校区', formatted = T)
#   longtitude     latitude 
# 116.35376693  40.00382807

如图1，在地图上显示解析出来的位置，可见定位在矿大1号楼。

cumtb_map <- getBaiduMap("中国矿业大学（北京）学院路校区", zoom = 17)
cumtb_coordinate <- getCoordinate("中国矿业大学（北京）学院路校区", formatted = T)
cumtb_coordinate <- data.frame(t(cumtb_coordinate))
# 调 ggmap 绘制背景地图
ggmap::ggmap(cumtb_map) +
  ggplot2::geom_point(aes(x = longtitude, y = latitude),
    data = cumtb_coordinate, col = "red", size = 5
  )

图 1: 百度地图-中国矿业大学（北京）学院路校区

高德地图

高德和百度都是有甲级测绘资质的单位，意味着地图数据符合国家要求，来源权威可用。大家熟知的全球定位系统（GPS）采用 WGS 84，而 GCJ 02 由 WGS 84 加密后的坐标系，是由中国国家测绘局制定的地理坐标系统，又称火星坐标系。高德地图采用火星坐标系 GCJ 02（国测局 Guo Ce Ju，即国家测绘局），百度地图的坐标系 BD 09 在 GCJ 02 坐标系基础上再次加密，加密过程不可逆，是一种非线性加密方式，反向解密后的坐标在部分区域的定位差别很大。另外，值得注意的是对同一目标不同公司提供的定位可能是不同的，因为所选的地理参考系不同，常用的三种坐标系见文档坐标系。

高德提供很多调用地图服务的 API 接口，如将其它坐标转化为高德坐标的坐标转化服务，前提是先在高德开放平台上申请Web服务API类型KEY。继续以地标「中国矿业大学（北京）学院路校区」为例，下面在 R 语言环境中，使用 httr 包(Wickham 2022b)调用高德地图服务将百度坐标转化为高德坐标。

library(httr)
# 向高德地图 API 发 GET 请求
tmp <- GET(
  url = "https://restapi.amap.com/",
  path = "v3/assistant/coordinate/convert",
  query = list(
    # 原坐标，经纬度小数点后不得超过6位
    locations = "40.003828,116.353766", 
    # 原坐标参考系，还支持 GPS 坐标，取值 'gps'
    coordsys = "baidu", 
    # 返回数据格式，还支持 xml
    output = "json", 
    # 高德地图 WEB 服务 API 访问令牌
    key = Sys.getenv("AMAP_KEY") 
  )
)
# 抽取转化后的高德坐标
content(tmp)$locations
# [1] "39.997202126977,116.347817690225"

接着，采用 leaflet 包将高德坐标系下的经纬度绘制在地图上。由于 leaflet 包默认采用 OpenStreetMap 开放街道地图服务，需要切换一下地图服务，而郎大为开发的 leafletCN 包(Lang 2017) 正好用函数 amap() 封装了高德瓦片地图服务。如图2所示，定位在矿大学10楼，从百度地图转化过来，这有上百米的距离。一方面开放的高德地图 API 不支持转化高精度的经纬度坐标，另一方面转化是有损失的。

library(leaflet)
library(leafletCN)
# 绘图
leaflet(options = leafletOptions(
  minZoom = 4, maxZoom = 18, 
  zoomControl = FALSE
)) |> 
  amap() |> 
  setView(lng = 116.347817690225, lat = 39.997202126977, zoom = 18) |> 
  addCircles(
    data = data.frame(lon = 116.347817690225, lat = 39.997202126977, size = 10),
    lng = ~lon, lat = ~lat, radius = ~size, color = "red",
    fillOpacity = 0.55, stroke = T, weight = 1,
    label = htmltools::HTML("这里是 <b>学院路校区</b>, 中国矿业大学（北京）")
  )

图 2: 高德地图-中国矿业大学（北京）学院路校区

下面试试看，将原百度坐标不转化，直接绘制在高德地图上，如图3所示，可见定位到中国农业大学去了，真是失之毫厘，差之千里！

# 绘图
leaflet(options = leafletOptions(
  minZoom = 4, maxZoom = 18, 
  zoomControl = FALSE
)) |> 
  amap() |> 
  setView(lng = 116.35376693, lat = 40.00382807, zoom = 18) |> 
  addCircles(
    data = data.frame(lon = 116.35376693, lat = 40.00382807, size = 10),
    lng = ~lon, lat = ~lat, radius = ~size, color = "red",
    fillOpacity = 0.55, stroke = T, weight = 1,
    label = htmltools::HTML("这里是 <b>学院路校区</b>, 中国矿业大学（北京）")
  )

图 3: 高德地图-中国矿业大学（北京）

必应地图

与 ggmap 包相比，RgoogleMaps包 (Loecher and Ropkins 2015) 支持调用更多的地图服务，除了 Google 地图外，还支持必应地图，以及基于 lattice 包的可视化，借助 loa (Ropkins 2021)和 latticeExtra (Sarkar and Andrews 2019)还可以支持高级的自定义，如何使用见ggmap 使用案例、loa 文档和latticeExtra 文档。

library(RgoogleMaps) # 同时需要 RCurl
library(RCurl)

使用 Bing 地图服务也是需要 API_KEY 的，可先去 Bing 地图开发中心 用 Outlook 邮箱创建账户，然后申请免费的基础版，过程见说明。继续以地标「中国矿业大学（北京）学院路校区」为例，将百度坐标绘制在必应地图上，如图4所示，也定位在中国农业大学校园内了，而且显示的兴趣点 POI 少很多了！

# 保存背景地图
map_cumtb <- GetBingMap(
  center = c(lat = 40.00382807, lon = 116.35376693),
  zoom = 17, apiKey = Sys.getenv("BingMap_API_KEY"),
  verbose = 0, destfile = "cumtb-bing.png"
)
# 将点绘制到地图上
PlotOnStaticMap(map_cumtb,
  lat = 40.00382807, lon = 116.35376693,
  col = "red", FUN = points,
  pch = 19, cex = 2
)

图 4: 必应地图-中国矿业大学（北京）学院路校区

谷歌地图

David Kahle 开发的ggmap (Kahle and Wickham 2013) 包在 2011 年发布在 CRAN 上，是一个基于 ggplot2 的空间可视化包，背景地图来自谷歌地图。目前，由于 Google 地图服务策略调整，需要申请 Google 地图服务的访问令牌才可以使用地图服务。 可将申请到的 GGMAP_GOOGLE_API_KEY 保存到本地文件 ~/.Renviron，以供 ggmap 后续使用。

library(ggmap)
p <- ggmap::get_googlemap(
  center = c(lon = 116.35376693, lat = 40.00382807),
  zoom = 12, size = c(640, 640), scale = 2
)
ggmap::ggmap(p)

Google 提供一种瓦片地图，即由一张张小图片拼凑起来形成一张完整的地图，就好像盖房用的瓦片，精度由瓦片的分辨率决定，如图 5 所示。

图 5: ggmap 地形图

ggmap 也可以调用多种风格样式的地图，地形图、水彩图等，如图6 所示。

library(ggplot2)
library(ggmap)
library(patchwork)

us <- c(left = -125, bottom = 25.75, right = -67, top = 49)
# 从 tile server 服务器上下载 Stamen Maps 数据，组成一个图片
dat1 <- get_stamenmap(
  bbox = us, zoom = 5, maptype = "toner-lite"
  )
p1 <- ggmap(dat1)
# 地形图
dat2 <- get_stamenmap(
  bbox = us, zoom = 5, maptype = "terrain"
)
p2 <- ggmap(dat2)

dat3 <- get_stamenmap(
  bbox = us, zoom = 5, maptype = "watercolor"
)
p3 <- ggmap(dat3)

dat4 <- get_stamenmap(
  bbox = us, zoom = 5, maptype = "terrain-background"
)
p4 <- ggmap(dat4)

p1 / p2 | p3 / p4

图 6: ggmap 几种常用背景地图

网址定位

Os Keyes 开发的 R 包 rgeolocate (Keyes et al. 2021) 可以解析 IP 地址，根据 IP 地址数据库提取位置信息，其内置 GeoLite2-Country.mmdb 数据支持定位到国家，是一个 MaxMind 的 GeoIP2 数据库文件。

library(rgeolocate)
# IP 地址 111.203.130.69 在中国矿业大学（北京）校区内。
maxmind(
  ips = "111.203.130.69",
  file = system.file("extdata", "GeoLite2-Country.mmdb", package = "rgeolocate"),
  fields = c(
    "continent_name", "country_code", "country_name",
    "region_name", "city_name", "city_geoname_id",
    "timezone", "longitude", "latitude"
  )
)
#   continent_name country_code country_name region_name city_name city_geoname_id
# 1           Asia           CN        China        <NA>      <NA>              NA
#   timezone longitude latitude
# 1     <NA>        NA       NA

依次是洲、国家编码简称、国家名称、区域（省级）名称、城市名称、城市代码（不确定编码体系，类似国家行政区划编码）、时区、经度和纬度。内置的数据库定位精度有限，可以去MaxMind 官网下载免费的可以定位到城市的数据库。rgeolocate 包还提供连接其它 IP 定位服务的接口，比如前面提及的 https://ipinfo.io/。

ip_info('111.203.130.69')
#      city  region country              loc                                          org
# 1 Beijing Beijing      CN 39.9075,116.3972 AS4808 China Unicom Beijing Province Network
#        timezone hostname postal phone
# 1 Asia/Shanghai       NA     NA    NA

依次将城市、省份、国家、经纬度、网络服务提供商和时区都提供了。除了 rgeolocate 包，还可以使用 httr 包向 ipinfo.io 发 GET 请求，获取定位信息。

library(httr)
ip_geocode <- GET(url = "ipinfo.io")
content(ip_geocode)$loc
# [1] "39.9075,116.3972"

或者调用高德地图 IP 定位服务，定位到省、市，并返回一个两个坐标表示的矩形区域。

# 向高德地图 API 发 GET 请求
library(httr)
tmp <- GET(
    url = "https://restapi.amap.com/",
    path = "v3/ip",
    query = list(
      ip = "111.203.130.69", # IP 地址
      output = "json",       # 返回数据格式
      # 高德地图 WEB 服务 API 访问令牌
      key = Sys.getenv("AMAP_KEY")
    )
  )
# 查看全部返回信息
# content(tmp)
# 抽取位置区域
content(tmp)$rectangle
# [1] "116.0119343,39.66127144;116.7829835,40.2164962"

瓦片地图

使用高德静态地图遵守自定义地图服务协议， 应国家要求，国内的地图服务厂商必须是 GCJ-02 或 BD-09 坐标系统，用户必须拿着这两个坐标系下的空间数据，调用各大 Web 地图服务制图，比如地图瓦片，它们采用 Web 墨卡托投影，坐标参考系是 EPSG:3857 或某种再次转化加密过的 Web 墨卡托投影。

继续以地标「中国矿业大学（北京）学院路校区」为例，之前已经获得它的定位 116.347817,39.997202，调用高德静态地图服务，返回一个静态图片地址。

library(httr)
GET(
  url = "https://restapi.amap.com/",
  path = "v3/staticmap",
  query = list(
    location = "116.347817,39.997202", # 图的中心位置
    zoom = 16, # 缩放水平
    size = "600*400", # 图片长宽
    scale = 2, # 返回高清图
    labels = sprintf("%s,%s,%s,%s,%s,%s:%s",
      content = "矿大北京", # 标签内容
      font = 0, # 微软雅黑
      bold = 1, # 粗体
      fontSize = 24, # 字体大小
      fontColor = "0xFF0000",  # 字体颜色：红色
      background = "0xFFFFFF", # 字体背景：白色
      location = "116.347817,39.997202" # 标签位置
    ),
    # 高德地图 WEB 服务 API 访问令牌
    key = Sys.getenv("AMAP_KEY")
  )
)

将返回的链接输入浏览器，会得到一张如下的图片。

图 7: 中国矿业大学（北京）学院路校区

地理编码

地球上每一块区域都对应有一个唯一的 GeoHash 值，八位 GeoHash 码的定位精度可达400平米，更多详细的介绍可见维基百科Geohash，比起区、县、乡、镇行政单元，它是用来刻画规则的更小的网格单元，广泛用于本地生活服务场景，比如外卖、打车、配送等。geohashTools 包(Chirico 2020)可以非常高效地将 GeoHash 编码的地理区域转为经纬度坐标。

# 随意造几个位置
geohash_df <- data.frame(geohash = c(
  "wmd0k2ktr", "wmd0k2kv3", "wmd0k2kv7",
  "wmd0k2kvq", "wmd0k2mj0", "wmd0k2mj4"
))
# 转码为经纬度
dat <- geohashTools::gh_decode(geohashes = geohash_df$geohash)
# 将位置绘制在地图上
leaflet::leaflet() |> 
  leafletCN::amap() |> 
  leaflet::addCircles(lng = dat$longitude, lat = dat$latitude, radius = 3)

图 8: GeoHash 地理区域编码

而 lwgeom 包(E. Pebesma 2021)提供函数 st_geohash() 可将经纬度坐标转为 GeoHash 码。在数据规模很大的情况下，这编码解码的过程方便聚合 GeoHash 粒度的数据，并呈现在地图上去观察数据中的空间规律，比如一些局部热点效应。

library(sf)
# Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(lwgeom)
# Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.10.2, PROJ 8.2.1
st_geohash(st_sfc(st_point(c(116.35688, 40.00314)), st_point(c(0, 90))), 9)
# [1] "wx4exf27p" "upbpbpbpb"

坐标系统

地球既不是圆球的也不是标准的椭球，将三维的球面投影到二维的平面有不同的方法，每一种方法都对应一种坐标转化，选定坐标原点后即构成坐标参考系统（Coordinate Reference System，简称 CRS）。maps 包(Brownrigg 2021)内置了一份世界地图数据，它采用 WGS 84 坐标系，此坐标系的详细描述，读者可在 R 软件 控制台运行 sf::st_crs("EPSG:4326") 获取。下面以此为例介绍，本节空间数据的操作都用 sf包(E. Pebesma 2022)来完成。

library(sf)
library(maps)
# st_as_sf 将 map 类转化为 sf 类
world1 <- st_as_sf(map("world", plot = FALSE, fill = TRUE))
## 检索坐标参考系统
st_crs(world1)
# Coordinate Reference System:
#   User input: EPSG:4326 
#   wkt:
# GEOGCRS["WGS 84",
#     ENSEMBLE["World Geodetic System 1984 ensemble",
#         MEMBER["World Geodetic System 1984 (Transit)"],
#         MEMBER["World Geodetic System 1984 (G730)"],
#         MEMBER["World Geodetic System 1984 (G873)"],
#         MEMBER["World Geodetic System 1984 (G1150)"],
#         MEMBER["World Geodetic System 1984 (G1674)"],
#         MEMBER["World Geodetic System 1984 (G1762)"],
#         MEMBER["World Geodetic System 1984 (G2139)"],
#         ELLIPSOID["WGS 84",6378137,298.257223563,
#             LENGTHUNIT["metre",1]],
#         ENSEMBLEACCURACY[2.0]],
#     PRIMEM["Greenwich",0,
#         ANGLEUNIT["degree",0.0174532925199433]],
#     CS[ellipsoidal,2],
#         AXIS["geodetic latitude (Lat)",north,
#             ORDER[1],
#             ANGLEUNIT["degree",0.0174532925199433]],
#         AXIS["geodetic longitude (Lon)",east,
#             ORDER[2],
#             ANGLEUNIT["degree",0.0174532925199433]],
#     USAGE[
#         SCOPE["Horizontal component of 3D system."],
#         AREA["World."],
#         BBOX[-90,-180,90,180]],
#     ID["EPSG",4326]]

sf 包提供 st_as_sf() 函数将 map 类转化为 sf 类，st_crs() 函数检索地图数据中附带的坐标参考系信息。接下来，用 st_transform() 函数转化坐标，参数 crs 指定新的坐标系统，示例里 +proj=laea 表示 Lambert Azimuthal Equal Area 投影，而 +ellps=WGS84 即 WGS 1984，或称 EPSG:4326。+lon_0=155 和 +lat_0=-90 分别是经纬度，东经 155 度，南纬 90 度，也就是观测视角到南极去了。

# 坐标转化
world2 <- st_transform(
  x = world1,
  crs = st_crs("ESRI:102020")
)
## 检索坐标参考系统
st_crs(world2)
# Coordinate Reference System:
#   User input: ESRI:102020 
#   wkt:
# PROJCRS["South_Pole_Lambert_Azimuthal_Equal_Area",
#     BASEGEOGCRS["WGS 84",
#         DATUM["World Geodetic System 1984",
#             ELLIPSOID["WGS 84",6378137,298.257223563,
#                 LENGTHUNIT["metre",1]]],
#         PRIMEM["Greenwich",0,
#             ANGLEUNIT["Degree",0.0174532925199433]]],
#     CONVERSION["South_Pole_Lambert_Azimuthal_Equal_Area",
#         METHOD["Lambert Azimuthal Equal Area",
#             ID["EPSG",9820]],
#         PARAMETER["Latitude of natural origin",-90,
#             ANGLEUNIT["Degree",0.0174532925199433],
#             ID["EPSG",8801]],
#         PARAMETER["Longitude of natural origin",0,
#             ANGLEUNIT["Degree",0.0174532925199433],
#             ID["EPSG",8802]],
#         PARAMETER["False easting",0,
#             LENGTHUNIT["metre",1],
#             ID["EPSG",8806]],
#         PARAMETER["False northing",0,
#             LENGTHUNIT["metre",1],
#             ID["EPSG",8807]]],
#     CS[Cartesian,2],
#         AXIS["(E)",north,
#             MERIDIAN[90,
#                 ANGLEUNIT["degree",0.0174532925199433]],
#             ORDER[1],
#             LENGTHUNIT["metre",1]],
#         AXIS["(N)",north,
#             MERIDIAN[0,
#                 ANGLEUNIT["degree",0.0174532925199433]],
#             ORDER[2],
#             LENGTHUNIT["metre",1]],
#     USAGE[
#         SCOPE["Not known."],
#         AREA["Southern hemisphere."],
#         BBOX[-90,-180,0,180]],
#     ID["ESRI",102020]]

下面将两个不同坐标系统下的世界地图绘制出来，如图 9 所示。

library(ggplot2)
library(patchwork)
p1 <- ggplot() +
  geom_sf(data = world1)

p2 <- ggplot() +
  geom_sf(data = world2)

p1 + p2

图 9: WGS 84 和坐标系统

散点图

R 语言数据操作常常离不开 data.frame，有很多空间数据也是按照一张张二维表格存储的，为了能够使用 sf 包及其生态环境，将 data.frame 类型转化为空间数据类型 sf 是非常必要的。之前，介绍过如何用 sp 包将 data.frame 转化为 sp 类型，进一步，sf 包提供的泛型函数 st_as_sf() 支持将多种数据类型转化为 sf 类型。目前支持 12 种数据类型转化，基本覆盖日常所需。

methods("st_as_sf")
#  [1] st_as_sf.data.frame*   st_as_sf.lpp*          st_as_sf.map*         
#  [4] st_as_sf.owin*         st_as_sf.ppp*          st_as_sf.ppplist*     
#  [7] st_as_sf.psp*          st_as_sf.s2_geography* st_as_sf.sf*          
# [10] st_as_sf.sfc*          st_as_sf.Spatial*      st_as_sf.SpatVector*  
# [13] st_as_sf.wk_crc*       st_as_sf.wk_rct*       st_as_sf.wk_wkb*      
# [16] st_as_sf.wk_wkt*       st_as_sf.wk_xy*       
# see '?methods' for accessing help and source code

另外，值得一提的是sfheaders包(Cooley 2020)，它可以非常方便高效地将普通数据框 data.frame 类型转化为空间数据 sf（Simple feature）类型，据其官网介绍，转化性能可以和 data.table 媲美，上游依赖很少，甚至不需要依赖 sf 包。

library(sf)
# library(sfheaders)
# beijing_sf <- sf_point(beijing, x = "x", y = "y", keep = TRUE)
# 将 data.frame 转化为 sf 类型
beijing_sf <- st_as_sf(beijing, coords = c("x", "y"))
# 指定原数据的坐标参考系
st_crs(beijing) <- 3857
# 转到日常用的经纬度坐标参考系
beijing_sf <- st_transform(beijing_sf, crs = 4326)
# 快速绘图查看数据
plot(beijing_sf["den"], pch = 19, cex = 0.2)

图 10: 北京市热点分布图

ggplot2 早已支持 sf 数据类型，不强调精雕细琢，绘图过程也十分简单。

library(ggplot2)
ggplot() +
  geom_sf(data = beijing_sf, aes(color = den), cex = 0.2)

图 11: 北京市热点分布图

rasterly 渲染 35058 多个点毫无压力，使用方式和 ggplot2 非常相近，配合 plotly 包支持渲染交互式图形，如图12所示。rasterly 的基本想法来自 Python 模块datashader，将数据点映射到网格中，聚合每个网格中的点得到代表点，即像素点，接着照常绘图，简而言之，连续空间离散化以近似，空间点数据转栅格数据。

library(rasterly)
# 静态散点图
rasterly(data = beijing, aes(x = x, y = y)) |> 
  rasterly_points() 
# 静态散点图：点密度、图片长宽
rasterly(data = beijing, aes(x = x, y = y, color = den), 
         plot_width = 700, plot_height = 500) |> 
  rasterly_points()
# 也可以绘制交互式的热力图
LogTransform <- function(z) {
  log(z)
}
plotly::plot_ly(beijing, x = ~x, y = ~y) |>
  add_rasterly_heatmap(
    on = ~den, # 对 den 变量取对数
    reduction_func = "max", # 聚合方式是求极大
    scaling = LogTransform, # 对数变换
    plot_width = 600, plot_height = 400
  ) |> 
  plotly::layout(
    xaxis = list(
      title = "x"
    ),
    yaxis = list(
      title = "y"
    )
  )

图 12: 北京市热点分布图

最后，补充介绍一下 sp 包绘制图形的过程，短期内恐还有必要。先将 data.frame 数据类型转为 sp 空间数据类型，以便调用更加高效且丰富的数据处理和可视化函数。

# 调 sp 包转化为 sp 类型
library(sp)
beijing_sp <- beijing
# 指定坐标变量
coordinates(beijing_sp) <- ~ x + y
# 指定参考系 EPSG:3857
proj4string(beijing_sp) <- CRS("EPSG:3857")
# 查看数据
summary(beijing_sp)
# Object of class SpatialPointsDataFrame
# Coordinates:
#        min      max
# x 12855679 13079460
# y  4761527  4985306
# Is projected: TRUE 
# proj4string :
# [+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m
# +nadgrids=@null +wktext +no_defs]
# Number of points: 35058
# Data attributes:
#       den       
#  Min.   : 1.00  
#  1st Qu.: 1.00  
#  Median : 2.00  
#  Mean   : 2.48  
#  3rd Qu.: 3.00  
#  Max.   :81.00
# 转化坐标参考系 EPSG:4326
beijing_sp <- spTransform(beijing_sp, CRSobj = CRS("EPSG:4326"))
# 查看转化后的数据
summary(beijing_sp)
# Object of class SpatialPointsDataFrame
# Coordinates:
#      min    max
# x 115.48 117.49
# y  39.28  40.82
# Is projected: FALSE 
# proj4string : [+proj=longlat +datum=WGS84 +no_defs]
# Number of points: 35058
# Data attributes:
#       den       
#  Min.   : 1.00  
#  1st Qu.: 1.00  
#  Median : 2.00  
#  Mean   : 2.48  
#  3rd Qu.: 3.00  
#  Max.   :81.00

# 最后绘图
spplot(beijing_sp, "den",
  colorkey = TRUE,     # 连续型图例
  alpha = 0.4,         # 散点透明度
  key.space = "right", # 图例位置
  cex = 0.3,           # 点的大小
  aspect = 0.7,        # 图形长宽比例
  scales = list(
    draw = T,
    # 去掉图形上边、右边多余的刻度线
    x = list(alternating = 1, tck = c(1, 0)),
    y = list(alternating = 1, tck = c(1, 0))
  )
)

图 13: 北京市热点分布图

热力图

OpenStreetMap 是开放街道地图数据，由来自世界各地的贡献者维护，但不提供免费的地图 API。使用此地图服务，不会将自己的数据上传至 OpenStreetMap 的服务器上。而 Leaflet 是开源的交互式地理可视化 JS 库。Joe Cheng 开发维护的 leaflet 包(Cheng, Karambelkar, and Xie 2022)是 JavaScript 库 Leaflet 的 R 语言接口。一个提供基础地图服务，一个提供数据渲染可视化，一个提供封装好的 R 语言接口。

图 14: 开源交互式 JavaScript 库 Leaflet

library(leaflet)
# 提供 addHeatmap 函数绘制热力图 
library(leaflet.extras) 
leaflet() |>
  addTiles() |>
  addHeatmap(
    data = beijing,
    lng = ~x, lat = ~y,
    intensity = ~den, # 密度
    blur = 20, max = 0.05, radius = 15
  ) |>
  addScaleBar(position = "bottomleft")

图 15: 北京市热点分布图

mapdeck 包将 MapBox GL 的地图服务和deck.gl的大规模可视化能力封装到一起，功能还是很强的，函数 add_heatmap() 可用来绘制热力图，数万散点毫无压力，如图16。

mapdeck(style = mapdeck_style("dark"), pitch = 45) |> 
  add_heatmap(
    data = beijing, lat = "y", lon = "x",
    weight = "den", colour_range = hcl.colors(6)
  )

图 16: 北京市热点分布图

蜂窝图

除了热力图，mapdeck 调用 add_hexagon() 和 add_screengrid() 还可以绘制蜂窝图，如图17。

mapdeck(style = mapdeck_style("dark"), pitch = 45) |> 
  add_screengrid(
    data = beijing,
    lat = "y",
    lon = "x",
    weight = "den",
    layer_id = "screengrid_layer",
    cell_size = 5,
    opacity = 0.6,
    colour_range = hcl.colors(6)
  )

图 17: 北京市热点分布图

柱状图

咋一看，读者可能会疑惑，怎么柱状图也能归到地理可视化？它在蜂窝图的基础上，用柱子的高低表示密度值大小，整个地图是三维的，可拖拽。mapdeck 提供函数 add_grid() / add_hexagon() 绘制三维柱状图，如图18。

library(mapdeck)
mapdeck(style = mapdeck_style("dark"), pitch = 45) |> 
  add_hexagon(
    data = beijing,
    lat = "y", lon = "x", elevation = "den",
    layer_id = "hex_layer",
    elevation_scale = 30,
    colour_range = hcl.colors(6)
  )

图 18: 北京市热点分布图

飞线图

以 nycflights13 (Wickham 2021) 包提供航班数据做飞线图，数据由 Hadley Wickham 收集自OpenFlights Airports Database，从纽约三大机场 JFK, LGA 或 EWR 的准点起飞数据。

library(nycflights13)
data("flights") # 航班
data("airports") # 机场

提取2013年1月1日的航班数据，出发地 origin 和目的地 dest。数据集 airports 包含 FAA 机场代码，经纬度（经纬度精确到小数点后6位，坐标参考系 EPSG:4326）等信息。

# 提取部分航班数据
sub_flights <- subset(
  x = flights, subset = year == 2013 & month == 1 & day == 1,
  select = c("origin", "dest", "time_hour")
)
# 统计航班架次
sub_flights <- aggregate(formula = time_hour ~ origin + dest, data = sub_flights, FUN = length)
# 准备机场位置信息
sub_airports <- subset(
  x = airports, select = c("faa", "lon", "lat")
)

关联机场位置信息，获取起飞机场和落地机场的经纬度信息。

# 起飞机场
tmp <- merge(sub_flights, sub_airports, by.x = "origin", by.y = "faa")
colnames(tmp) <- c("origin", "dest", "cnt", "origin_lon", "origin_lat")
# 落地机场
tmp <- merge(tmp, sub_airports, by.x = "dest", by.y = "faa")
colnames(tmp) <- c(
  "origin", "dest", "cnt", "origin_lon", "origin_lat",
  "dest_lon", "dest_lat"
)

mapdeck 提供 add_greatcircle() 图层绘制飞线图。

mapdeck(style = "mapbox://styles/mapbox/dark-v9") |> 
  add_greatcircle(
    data = tmp,
    layer_id = "arc_layer",
    origin = c("origin_lon", "origin_lat"),
    destination = c("dest_lon", "dest_lat"),
    stroke_from = "cnt",
    stroke_width = 3, # 弧线宽度
    palette = "viridis", # 调色板 "plasma"
    legend = TRUE
  ) |>  
  mapdeck_view(
    # 视角位置
    location = c(-110, 48),
    # 设置缩放水平
    zoom = 2,
    # 俯仰角
    pitch = 45
  )

图 19: 2013年1月1日纽约机场出发的航班

leaflet

library(leaflet)
# 人口密度分段
bins <- c(0, 10, 20, 50, 100, 200, 500, 1000, Inf)
# 构造调色板
pal <- colorBin("YlOrRd", domain = states$density, bins = bins)
# 准备悬浮提示
labels <- sprintf(
  "<strong>%s</strong><br/>%g people / mi<sup>2</sup>",
  states$name, states$density
) |>  lapply(htmltools::HTML)

# 绘图
leaflet(states) |> 
  setView(lng = -96, lat = 37.8, zoom = 4) |> 
  # 添加默认的 OSM 瓦片服务
  addTiles() |> 
  addPolygons(
    color = "white", # 边界线颜色
    opacity = 1,     # 颜色透明度
    weight = 2,      # 边界线宽度
    fillOpacity = 0.6,   # 填充色透明度
    fillColor = ~ pal(density),  # 填充色
    label = labels
  ) |> 
  addLegend(
    pal = pal, 
    values = ~density,
    opacity = 0.7, 
    title = "人口密度",
    position = "bottomright"
  )

图 21: leaflet 包专题地图

LeafLet JS 在专题地图choropleth示例里调用了 Mapbox 瓦片服务和地图风格。下面在函数 addTiles() 里替换瓦片服务 URL 模版路径，leaflet 也可用上 Mapbox 瓦片地图。

leaflet(states) |> 
  setView(lng = -96, lat = 37.8, zoom = 4) |> 
  addTiles(urlTemplate = sprintf(
    paste(
      "https://api.mapbox.com/styles/v1/mapbox/light-v9",
      "tiles/{z}/{x}/{y}?access_token=%s",
      sep = "/"
    ),
    Sys.getenv("MAPBOX_TOKEN") # 准备个人访问令牌
  ), attribution = sprintf(
    paste(
      "Map data &copy; <a href='%s'>OpenStreetMap</a>",
      "contributors, Imagery © <a href='%s'>Mapbox</a>"
    ),
    "https://www.openstreetmap.org/copyright",
    "https://www.mapbox.com/"
  )) |> 
  addPolygons(
    color = "white", # 边界线颜色
    opacity = 1, # 颜色透明度
    weight = 2, # 边界线宽度
    fillOpacity = 0.6, # 填充色透明度
    fillColor = ~ pal(density), # 填充色
    label = labels
  ) |> 
  addLegend(
    pal = pal,
    values = ~density,
    opacity = 0.7,
    title = "人口密度",
    position = "bottomright"
  )

图 22: leaflet 包专题地图

借助leaflet-providers，Leaflet 已经支持很多瓦片服务，除了以上两种，比如夜景图23。

library(leaflet)
leaflet() |> 
  addProviderTiles("NASAGIBS.ViirsEarthAtNight2012") |> 
  addMarkers(
    lng = 116.35376693,
    lat = 40.00382807,
    popup = "中国矿业大学（北京）学院路校区"
  )

图 23: leaflet 包调 NASA 瓦片服务

mapdeck

与 leaflet 要求类似，mapdeck 也要求空间数据的坐标参考系是 EPSG:4326。

library(mapdeck)
states <- transform(states, density_log10 = log10(density))
mapdeck() |> 
  add_polygon(
    data = states,
    fill_colour = "density_log10",
    fill_opacity = 0.8,
    palette = "plasma",
    legend = TRUE
  )

图 24: mapdeck 包专题地图

sf

地理单元的编码具有唯一性，即一套编码体系下，同一块地方不能有两个编码值，可以将其视为主键，则每一块区域上的人口密度就是属性值。

plot(states["density"], logz = TRUE)

图 25: sf 包专题地图

ggplot2

library(ggplot2)
ggplot() +
  geom_sf(data = states, aes(geometry = geometry, fill = density), colour = NA) +
  scale_fill_gradientn(
    trans = "log10", colors = sf.colors(10),
    guide = guide_colourbar(
      barwidth = 25, barheight = 0.4,
      title.position = "left"
    )
  ) +
  theme(
    panel.grid = element_line(colour = "gray90"),
    legend.position = "top",
    panel.background = element_blank()
  )

图 26: ggplot2 包专题地图