[原] 解密 Uber 数据团队的大规模地理数据可视化神器:Deck.gl 与 H3

clipboard.png

背景

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

地理单元:H3

clipboard.png

随着互联网出行公司的全球化扩张,越来越多的公司涌现出对地理单元划分的需求。

一方面,传统的地理单元比如 S2和geohash,在不同纬度的地区会出现地理单元单位面积差异较大的情况,这导致业务指标和模型输入的特征存在一定的分布倾斜和偏差,使用六边形地理单元可以减少指标和特征normalization的成本。

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

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

clipboard.png

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

clipboard.png

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

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

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

#Include libraries
library(dggridR)
library(dplyr)

#Construct a global grid with cells approximately 1000 miles across
dggs <- dgconstruct(spacing=1000, metric=FALSE, resround='down')

#Load included test data set
data(dgquakes)

#Get the corresponding grid cells for each earthquake epicenter (lat-long pair)
dgquakes$cell <- dgGEO_to_SEQNUM(dggs,dgquakes$lon,dgquakes$lat)$seqnum

#Converting SEQNUM to GEO gives the center coordinates of the cells
cellcenters <- dgSEQNUM_to_GEO(dggs,dgquakes$cell)

#Get the number of earthquakes in each cell
quakecounts <- dgquakes %>% group_by(cell) %>% summarise(count=n())

#Get the grid cell boundaries for cells which had quakes
grid <- dgcellstogrid(dggs,quakecounts$cell,frame=TRUE,wrapcells=TRUE)

#Update the grid cells' properties to include the number of earthquakes
#in each cell
grid <- merge(grid,quakecounts,by.x="cell",by.y="cell")

#Make adjustments so the output is more visually interesting
grid$count <- log(grid$count)
cutoff <- quantile(grid$count,0.9)
grid <- grid %>% mutate(count=ifelse(count>cutoff,cutoff,count))

#Get polygons for each country of the world
countries <- map_data("world")

#Plot everything on a flat map
p<- ggplot() + 
 geom_polygon(data=countries, aes(x=long, y=lat, group=group), fill=NA, color="black") +
 geom_polygon(data=grid, aes(x=long, y=lat, group=group, fill=count), alpha=0.4) +
 geom_path (data=grid, aes(x=long, y=lat, group=group), alpha=0.4, color="white") +
 geom_point (aes(x=cellcenters$lon_deg, y=cellcenters$lat_deg)) +
 scale_fill_gradient(low="blue", high="red")
p

clipboard.png

转化坐标系后:

#Replot on a spherical projection
p+coord_map("ortho", orientation = c(-38.49831, -179.9223, 0))+
  xlab('')+ylab('')+
  theme(axis.ticks.x=element_blank())+
  theme(axis.ticks.y=element_blank())+
  theme(axis.text.x=element_blank())+
  theme(axis.text.y=element_blank())+
  ggtitle('Your data could look like this')

clipboard.png

在 H3 开源后,你也可以使用 h3r 实现:

# 以亮马桥地铁站为例
devtools::install_github("scottmmjackson/h3r")
library(h3r)

df <- h3r::getBoundingHexFromCoords(39.949958,116.46343,11) %>% # 单边长为24米
 purrr::transpose() %>% 
 purrr::simplify_all() %>%
 data.frame()

df %>% bind_rows(
 df %>% head(1)
) %>% 
 leaflet::leaflet() %>% 
 leafletCN::amap() %>% 
 leaflet::addPolylines(lng = ~lon,lat=~lat)

clipboard.png

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

地理数据可视化:Deck.gl

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

不过随着 Leaflet 也暴露出一些新的问题,比如如何大规模渲染地理数据,支持诸如 轨迹、风向、六边形网格的可视化。好在近年来 Mapbox 和 Deck.gl 正在着手改变这一现状。

下面是一个具体的例子,如何可视化Hexagon:

# 初始化
devtools::install_github("crazycapivara/deckgl")

library(deckgl)

# 设置 Mapbox token,过期需要免费在 Mapbox 官网申请
Sys.setenv(MAPBOX_API_TOKEN = "pk.eyJ1IjoidWJlcmRhdGEiLCJhIjoiY2poczJzeGt2MGl1bTNkcm1lcXVqMXRpMyJ9.9o2DrYg8C8UWmprj-tcVpQ")


# 数据集合
sample_data <- paste0(
  "https://raw.githubusercontent.com/",
  "uber-common/deck.gl-data/",
  "master/website/sf-bike-parking.json"
)

properties <- list(
  pickable = TRUE,
  extruded = TRUE,
  cellSize = 200,
  elevationScale = 4,
  getPosition = JS("data => data.COORDINATES"),
  getTooltip = JS("object => object.count")
)

# 可视化
deckgl(zoom = 11, pitch = 45) %>%
  add_hexagon_layer(data = sample_data, properties = properties) %>%
  add_mapbox_basemap(style = "mapbox://styles/mapbox/light-v9") 

clipboard.png

除了六边形之外 Deck.gl 也支持其他常见几何图形,比如 Grid、Arc、Contour、Polygon 等等。
更多信息可以见官方文档: https://crazycapivara.github....

clipboard.png

地理仪表盘:结合 Shiny

Deck.gl 结合 Shiny 后,可将可视化结果输出到仪表盘上:

clipboard.png

library(mapdeck)
library(shiny)
library(shinydashboard)
library(jsonlite)
ui <- dashboardPage(
    dashboardHeader()
    , dashboardSidebar()
    , dashboardBody(
        mapdeckOutput(
            outputId = 'myMap'
            ),
        sliderInput(
            inputId = "longitudes"
            , label = "Longitudes"
            , min = -180
            , max = 180
            , value = c(-90, 90)
        )
        , verbatimTextOutput(
            outputId = "observed_click"
        )
    )
)
server <- function(input, output) {
    
    set_token('pk.eyJ1IjoidWJlcmRhdGEiLCJhIjoiY2poczJzeGt2MGl1bTNkcm1lcXVqMXRpMyJ9.9o2DrYg8C8UWmprj-tcVpQ') ## 如果token 过期了,需要去Mapbox官网免费申请一个
    
    origin <- capitals[capitals$country == "Australia", ]
    destination <- capitals[capitals$country != "Australia", ]
    origin$key <- 1L
    destination$key <- 1L
    
    df <- merge(origin, destination, by = 'key', all = T)
    
    output$myMap <- renderMapdeck({
        mapdeck(style = mapdeck_style('dark')) 
    })
    
    ## plot points & lines according to the selected longitudes
    df_reactive <- reactive({
        if(is.null(input$longitudes)) return(NULL)
        lons <- input$longitudes
        return(
            df[df$lon.y >= lons[1] & df$lon.y <= lons[2], ]
        )
    })
    
    observeEvent({input$longitudes}, {
        if(is.null(input$longitudes)) return()
        
        mapdeck_update(map_id = 'myMap') %>%
            add_scatterplot(
                data = df_reactive()
                , lon = "lon.y"
                , lat = "lat.y"
                , fill_colour = "country.y"
                , radius = 100000
                , layer_id = "myScatterLayer"
            ) %>%
            add_arc(
                data = df_reactive()
                , origin = c("lon.x", "lat.x")
                , destination = c("lon.y", "lat.y")
                , layer_id = "myArcLayer"
                , stroke_width = 4
            )
    })
    
    ## observe clicking on a line and return the text
    observeEvent(input$myMap_arc_click, {
        
        event <- input$myMap_arc_click
        output$observed_click <- renderText({
            jsonlite::prettify( event )
        })
    })
}
shinyApp(ui, server)

参考资料

作为分享主义者(sharism),本人所有互联网发布的图文均遵从CC版权,转载请保留作者信息并注明作者 Harry Zhu 的 FinanceR专栏:https://segmentfault.com/blog...,如果涉及源代码请注明GitHub地址:https://github.com/harryprince。微信号: harryzhustudio
商业使用请联系作者。

FinanceR
循环写作,持续更新,形成闭环,贵在坚持
2.2k 声望
2.2k 粉丝
0 条评论
推荐阅读
[译] 层次时间序列预测法
大多数关于时间序列预测的文章都侧重于特定的聚合程度。但是,当我们能够深入分析聚合的数据,以便在更细粒度的层次上观察同一个序列时,挑战就出现了。在这种情况下,我们往往会发现,对较低水平的预测与总体预...

HarryZhu阅读 2.7k

封面图
在Canvas中绘制Geojson数据
需求分析在做地图开发的时候遇到一个需求,是在 canvas 中绘制 Geojson 数据数据格式为 EPSG:4326 的 Polygon:三维数组每一项都是由经纬度组成的第一个点和最后一个点相同,表示 Polygon 是闭合的 {代码...} 需...

uccs阅读 466

封面图
ArcGIS QGIS学习二:图层如何只显示需要的部分几何面数据(附最新坐标边界下载全国省市区县乡镇)
当我们用GIS软件打开一个SHP文件的时候,会显示出里面全部的几何图形,假如我只想要其中的一部分数据显示出来,其他的均不要显示,有那么几种操作方法。

高坚果阅读 408

RPA界面元素智能自适应定位与操控技术
什么是RPA?RPA(Robotic Process Automation,机器人流程自动化)是通过特定的、可模拟人类在计算机界面上进行操作的技术,它可以按规则自动执行相应的流程任务,代替或辅助人类完成相关的计算机操作,从而节约...

达观数据阅读 335

封面图
ALLUVIAL DIAGRAM(冲积图)详解和R语言实现实例
冲积图是最初开发用来表示随时间变化的网络结构的一种流程图。为了兼顾它们的视觉外观和对流动的重视,冲积图是以流水堆积的土壤自然形成的冲积扇命名的。变量分配给平行的垂直轴。值由每个轴上的块表示。块的高...

拓端tecdat阅读 314

leaflet底图样式修改及热力图加载、地图打点
{代码...}

云胡不喜阅读 343

Y 分钟速成 R
源代码下载: learnr-zh.rR 是一门统计语言。它有很多数据分析和挖掘程序包。可以用来统计、分析和制图。 你也可以在 LaTeX 文档中运行 R 命令。 {代码...} 获得 R从 [链接] 获得安装包和图形化界面RStudio 是另...

小X学技术阅读 273

2.2k 声望
2.2k 粉丝
宣传栏