大规模地理数据可视化入门：Deck.gl 和 H3

如何大规模可视化地理数据一直都是一个业界的难点，随着2015年起 Uber 在这一领域的发力 (https://segmentfault.com/a/1190000005154321)，构建了基于 Deck.gl + H3 (deckgl (https://github.com/crazycapivara/deckgl), h3r (https://github.com/scottmmjackson/h3r/)) 的大规模数据可视化方案。一方面，极大地满足了日常前端开发者的需求。另一方面，也极大地方便了数据科学家的可视化工作。在大规模空间轨迹分析、交通流量与供需预测等领域这一方案正得到广泛应用，突破了传统方法中数据量（通常不会超过10W个原始点）的瓶颈问题，实现百万点绘制无压力，并且可以结合 GPU 实现加速渲染。

举例

具体而言，比如，在移动互联网中常见的一个场景，预测城市中每个区域的供给和需求。在这个过程中，通常需要将预测的结果进行可视化以追踪模型的表现和问题。

Uber Movement 展示供需情况

在传统的数据量下（比如以行政区或街道为单位）通常是1000个多边形或者点的数量级进行城市数据的汇总或者采样，在这种数据量条件下直接查询和渲染数据并不会出现瓶颈问题。

但是，随着移动互联网的发展，涌现出了精度更高的需求（10米或100米），以求更精细地追踪用户在末端网络上的需求。通常这种场景都是至少百万个多边形或者点的数据量级，在这种条件下粗暴地直接查询和渲染是完全不可行的。

为了实现大规模地理数据可视化中，需要从 空间划分 和 渲染引擎 两方面入手。

空间划分

什么是空间划分

通过空间划分，建立数据库索引可以实现高效的实时查询服务，比如周边车辆位置查询(https://segmentfault.com/a/1190000008657566#articleHeader1)，周边餐厅位置查询 (https://tech.meituan.com/2014/09/05/lucene-distance.html) 等。

传统的空间划分方法主要分为两类：

1. 以 R-tree 为代表的动态单元

   (https://en.wikipedia.org/wiki/R-tree)

2. 以 Geohash 为代表的静态单元 

   (https://en.wikipedia.org/wiki/Geohash)

R-tree 简单来说是将空间划分为若干个不规则的边界框矩形的b+树索引，适用于面、线数据，查询时间复杂度为 O(n)。

R-tree 示意图

Geohash 简单来说是将二维的经纬度转换成字符串的四叉树线性索引，适用于点数据，查询时间复杂度为 O(log(N))

Geohash 示意图

下面是动态单元与静态单元的对比：

动态单元通常应用于对精度要求苛刻的场景，比如共享单车的违停区域判罚。

在本例中，供需预测对于几何边缘并没有高精度要求，所以牺牲一部分尽可能小的索引精度来换取计算性能是被允许的，此时，静态单元是优于动态单元的选择。

静态单元对比

常见的静态单元除 geohash 以外还有 S2。其中 S2 和 geohash 非常类似，也是基于四叉树的一种方法，只是在填充空间时使用了希尔伯特曲线而不是 geohash 中的Z阶曲线使得索引更加稳定，二者的详细原理分析见高效的多维空间点索引算法—— Geohash 和 Google S2 (https://halfrost.com/go_spatial_search/)。

静态地理单元特点对比

但是，一方面，传统的地理单元比如 S2 和 geohash，在国际化业务中却存在一个致命缺陷：在不同纬度的地区会出现地理单元单位面积差异较大的情况，比如北京和新加坡的 geohash 对应面积有将近30%的差异。这导致业务指标和模型输入的特征存在一定的分布倾斜和偏差，使用等面积、等形状的六边形地理单元可以减少指标和特征 normalization 的成本。

另一方面，在常用的地理范围查询中，基于矩形的查询方法，存在8邻域到中心网格的距离不相等的问题，四边形存在两类长度不等的距离，而六边形的周围邻居到中心网格的距离却是有且仅有一个，从形状上来说更加接近于圆形。

所以，基于hexagon的地理单元已经成为各大厂家的首选，比如 Uber 和 Didi 的峰时定价服务。

图示为 Uber 峰时定价服务示意图

在这样的背景下 Uber 基于六边形网格的地理单元开源解决方案 H3 应运而生，它使得部署 Hexagon 方案的成本非常低，通过 UDF、R pacakge 等方式可以以非常低的成本大规模推广。

什么是 H3

DDGS

H3 的前身其实是 DDGS(Discrete global grid systems) 中的 ISEA3H，其原理是把无限的不规则但体积相等的六棱柱从二十面体中心延伸，这样任何半径的球体都会穿过棱镜形成相等的面积 cell，基于该标准使得每一个地理单元的面积大小就可以保证几乎相同。

然而原生的 ISEA3H 方案在任意级别中都存在12个五边形，H3 的主要改进是通过坐标系的调整将其中的五边形都转移到水域上，这样就不影响大多数业务的开展。

下面是 ISEA3H 五边形问题的示例：

 1# 加载相关包
 2library(dggridR)
 3library(dplyr)
 4# 构建公里网格
 5dggs <- dgconstruct(spacing = 1000, metric = FALSE, resround = "down")
 6# 加载测试数据集
 7data(dgquakes)
 8# 获取每个震源中心对应的网格
 9dgquakes$cell <- dgGEO_to_SEQNUM(dggs, dgquakes$lon, dgquakes$lat)$seqnum
10# 将 SEQNUM 转为网格中心
11cellcenters <- dgSEQNUM_to_GEO(dggs, dgquakes$cell)
12# 获取每个单元的地震次数
13quakecounts <- dgquakes %>% group_by(cell) %>% summarise(count = n())
14# 获取地震网格单元边界
15grid <- dgcellstogrid(dggs, quakecounts$cell, frame = TRUE, wrapcells = TRUE)
16# 更新网格单元的地震次数
17grid <- merge(grid, quakecounts, by.x = "cell", by.y = "cell")
18# Normarlize 指标便于展示
19grid$count <- log(grid$count)
20cutoff <- quantile(grid$count, 0.9)
21grid <- grid %>% mutate(count = ifelse(count > cutoff, cutoff, count))
22# 获取每个国家的多边形
23countries <- map_data("world")
24# 绘制地图
25p <- ggplot() +
26  geom_polygon(data = countries, aes(x = long, y = lat, group = group), fill = NA, color = "black") +
27  geom_polygon(data = grid, aes(x = long, y = lat, group = group, fill = count), alpha = 0.4) +
28  geom_path(data = grid, aes(x = long, y = lat, group = group), alpha = 0.4, color = "white") +
29  geom_point(aes(x = cellcenters$lon_deg, y = cellcenters$lat_deg)) +
30  scale_fill_gradient(low = "blue", high = "red")
31p

转化坐标系后：

1# 重新在球坐标上绘制
2p + coord_map("ortho", orientation = c(-38.49831, -179.9223, 0)) +
3  theme(
4    axis.ticks.x = element_blank(),
5    axis.ticks.y = element_blank(),
6    axis.text.x = element_blank(),
7    axis.text.y = element_blank()
8  ) +
9  labs(x = "", y = "", title = "Your data could look like this")

可以看到此时在若干个六边形中存在五边形的情形。

在 H3 开源后，也可以使用 h3r 实现六边形的编码与解码：

 1# 以亮马桥地铁站为例
 2devtools::install_github("scottmmjackson/h3r") #或者 devtools::install_github("harryprince/h3r", ref="bug-fix/Makefile")
 3library(h3r)
 4df <- h3r::getBoundingHexFromCoords(39.949958, 116.46343, 11) %>% # 单边长为24米
 5  purrr::transpose() %>%
 6  purrr::simplify_all() %>%
 7  data.frame()
 8df %>%
 9  bind_rows(
10    df %>% head(1)
11  ) %>%
12  leaflet::leaflet() %>%
13  leafletCN::amap() %>%
14  leaflet::addPolylines(lng = ~lon, lat = ~lat)

H3 中还提供了类似 S2 的六边形压缩技术，使得数据的存储空间可以极大压缩，在处理大规模稀疏数据时将体现出优势：

H3 数据压缩技术

渲染引擎

什么是 Leaflet

在使用 Deck.gl 之前，业界主流的解决方案通常是另一个开源的轻量级地理数据可视化框架 Leaflet。Leaflet 经过十余年的积累已经拥有足够成熟的生态，支持各式各样的插件扩展。

leaflet: 

https://github.com/rstudio/leaflet

leaflet.opacity: 

https://github.com/be-marc/leaflet.opacity

leaflet.extras: 

https://github.com/bhaskarvk/leaflet.extras

leaflet.esri: 

https://github.com/bhaskarvk/leaflet.esri

mapview: 

https://github.com/r-spatial/mapview

mapedit: 

https://github.com/r-spatial/mapedit

leafletCN: 

https://github.com/Lchiffon/leafletCN

虽然 Leaflet 功能强大，不过工业界的发展也暴露出一些新的问题。如何更好地支持诸如 轨迹、风向、三维空间、六边形网格的交互式可视化此前没有好的解决方案。好在近年来 Mapbox 和 Deck.gl 正在着手改变这一现状。

什么是 Deck.gl

Deck.gl (http://deck.gl) 基于 WebGL 的大规模数据可视化框架，通过响应式编程和GPU并行加速的方式进行高效地 WebGL 渲染，与 Mapbox GL 深度结合能够呈现非凡的 3D 视觉效果。

下面是一个具体的例子，如何以 Hexagon 可视化百万个样本点：

1# 初始化
2library(mapdeck)
3# 生成 百万数据样本点
4df = data.frame(lat = rnorm(1000000,40,1),lng =rnorm(1000000,160,1)) # 以二维正态生成随机数据
5# 渲染
6mapdeck::mapdeck(style = "mapbox://styles/mapbox/dark-v9",token = "pk.eyJ1IjoidWJlcmRhdGEiLCJhIjoiY2poczJzeGt2MGl1bTNkcm1lcXVqMXRpMyJ9.9o2DrYg8C8UWmprj-tcVpQ") %>%
7  mapdeck::add_hexagon(lon = "lng",lat="lat",data = df,elevation_scale = 1000)

除了六边形之外 Deck.gl 也支持其他常见几何图形，比如 Grid、Arc、Contour、Polygon 等等。

更多信息可以见官方文档：https://crazycapivara.github.io/deckgl/

地理仪表盘：结合 Shiny

Deck.gl 结合 Shiny 后，可将可视化结果输出到仪表盘上，举个例子：

 1library(mapdeck)
 2library(shiny)
 3library(shinydashboard)
 4library(jsonlite)
 5ui <- dashboardPage(
 6  dashboardHeader(),
 7  dashboardSidebar(), dashboardBody(
 8    mapdeckOutput(
 9      outputId = "myMap"
10    ),
11    sliderInput(
12      inputId =
13        "longitudes", label =
14        "Longitudes", min =
15        -180, max =
16        180, value = c(-90, 90)
17    ), verbatimTextOutput(
18      outputId = "observed_click"
19    )
20  )
21)
22server <- function(input, output) {
23  set_token("pk.eyJ1IjoidWJlcmRhdGEiLCJhIjoiY2poczJzeGt2MGl1bTNkcm1lcXVqMXRpMyJ9.9o2DrYg8C8UWmprj-tcVpQ") ## 如果token 过期了，需要去Mapbox官网免费申请一个
24  origin <- capitals[capitals$country == "Australia", ]
25  destination <- capitals[capitals$country != "Australia", ]
26  origin$key <- 1L
27  destination$key <- 1L
28  df <- merge(origin, destination, by = "key", all = T)
29  output$myMap <- renderMapdeck({
30    mapdeck(style = mapdeck_style("dark"))
31  })
32  ## plot points & lines according to the selected longitudes
33  df_reactive <- reactive({
34    if (is.null(input$longitudes)) return(NULL)
35    lons <- input$longitudes
36    return(
37      df[df$lon.y >= lons[1] & df$lon.y <= lons[2], ]
38    )
39  })
40  observeEvent({
41    input$longitudes
42  }, {
43    if (is.null(input$longitudes)) return()
44    mapdeck_update(map_id = "myMap") %>%
45      add_scatterplot(
46        data =
47          df_reactive(), lon =
48          "lon.y", lat =
49          "lat.y", fill_colour =
50          "country.y", radius =
51          100000, layer_id = "myScatterLayer"
52      ) %>%
53      add_arc(
54        data =
55          df_reactive(), origin =
56          c("lon.x", "lat.x"), destination =
57          c("lon.y", "lat.y"), layer_id =
58          "myArcLayer", stroke_width = 4
59      )
60  })
61  ## observe clicking on a line and return the text
62  observeEvent(input$myMap_arc_click, {
63    event <- input$myMap_arc_click
64    output$observed_click <- renderText({
65      jsonlite::prettify(event)
66    })
67  })
68}
69shinyApp(ui, server)

总结

目前，在空间划分上，H3 正在超越 S2/Geohash 成为新标准，相关生态也趋于成熟。在渲染引擎上，Deck.gl 在特定领域已经全面领先 Leaflet，相关产品不断涌现，比如对标 carto (https://carto.com) 的地理数据分析工具 kepler (http://kepler.gl) 和毫秒级 OLAP 交互式分析工具 OmniSci (https://www.omnisci.com/)。

OmniSci NYC Taxi Rides Demo

参考资料

• RStudio Spark/Leaflet 与 GIS 最佳实践 (https://segmentfault.com/a/1190000015048613)

• Uber H3 原理分析 (https://qiita.com/gshirato/items/d8cc928c4131f3292b14)

• http://strimas.com/spatial/hexagonal-grids/

• https://cran.r-project.org/web/packages/dggridR/vignettes/dggridR.html

• https://en.wikipedia.org/wiki/Discrete_Global_Grid

• http://www.pyxisinnovation.com/pyxwiki/index.php?title=How_PYXIS_Works

• Large Scale Data Visualisation with Deck.gl and Shiny (https://www.youtube.com/watch?v=s0k_Hwn5Slg)

• https://uber.github.io/h3/

• https://eng.uber.com/shan-he/

• https://eng.uber.com/keplergl/

• 译 解密 Uber 数据团队的车辆定位查询算法 (https://segmentfault.com/a/1190000008657566)

• 译 解密 Uber 数据部门的数据可视化最佳实践 (https://segmentfault.com/a/1190000005154321)

• gpu-accelerated-aggregation-in-deck-gl (https://medium.com/vis-gl/gpu-accelerated-aggregation-in-deck-gl-7e2c7d701fb0)

作者：朱俊辉

日期：2019-01-05

统计之都：专业、人本、正直的中国统计学社区。

关注方式：扫描下图二维码。或查找公众号，搜索 统计之都 或 CapStat 即可。

往期推送：进入统计之都会话窗口，点击右上角小人图标，查看历史消息即可。

继续阅读