[原] 解密 Uber 數據團隊的大規模地理數據可視化神器：Deck.gl 與 H3

背景

如何大規模可視化地理數據一直都是一個業界的難點，隨着2015年起 Uber 在這一領域的發力，構建了基於 Deck.gl + H3 (deckgl,h3r) 的大規模數據可視化方案。一方面，極大地知足了大規模地理數據可視化的需求。另外一方面，也極大地方便了數據科學家的可視化工做。在大規模空間軌跡分析、交通流量與供需預測等領域獲得普遍應用，突破了原來leaflet架構中數據量（一般不會超過10W個原始點）的瓶頸問題，實現百萬點繪製無壓力，而且能夠結合GPU實現加速渲染。

地理單元：H3

隨着互聯網出行公司的全球化擴張，愈來愈多的公司涌現出對地理單元劃分的需求。

一方面，傳統的地理單元好比 S2和geohash，在不一樣緯度的地區會出現地理單元單位面積差別較大的狀況，這致使業務指標和模型輸入的特徵存在必定的分佈傾斜和誤差，使用六邊形地理單元能夠減小指標和特徵normalization的成本。

另外一方面，在經常使用的地理範圍查詢中，基於矩形的查詢方法，存在8鄰域到中心網格的距離不相等的問題，也就是說六邊形網格與周圍網格的距離有且僅有一個，而四邊形存在兩類距離
，而六邊形的周圍鄰居到中心網格的距離倒是相等的，從形狀上來講更加接近於圓形。

因此，基於hexagon的地理單元已經成爲各大廠家的首選，好比 Uber 和 Didi 的峯時訂價服務。

在這樣的背景下 Uber 基於六邊形網格的地理單元開源解決方案 H3 應運而生，它使得部署 Hexagon 方案的成本很是低，經過UDF、R pacakge等方式能夠以很是低的成本大規模推廣。

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

轉化座標系後：

#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')

在 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)

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")

除了六邊形以外 Deck.gl 也支持其餘常見幾何圖形，好比 Grid、Arc、Contour、Polygon 等等。
更多信息能夠見官方文檔： https://crazycapivara.github....

地理儀表盤：結合 Shiny

Deck.gl 結合 Shiny 後，可將可視化結果輸出到儀表盤上：

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)

參考資料

• RStudio Spark/Leaflet 與 GIS 最佳實踐
• Uber H3 原理分析
• http://strimas.com/spatial/he...
• https://cran.r-project.org/we...
• https://en.wikipedia.org/wiki...
• http://www.pyxisinnovation.co...
• Large Scale Data Visualisation with Deck.gl and Shiny
• https://uber.github.io/h3/
• https://eng.uber.com/shan-he/
• https://eng.uber.com/keplergl/
• [譯] 解密 Uber 數據部門的數據可視化最佳實踐
• gpu-accelerated-aggregation-in-deck-gl
• awesome leaflet

做爲分享主義者(sharism)，本人全部互聯網發佈的圖文均聽從CC版權，轉載請保留做者信息並註明做者 Harry Zhu 的 FinanceR專欄: https://segmentfault.com/blog...，若是涉及源代碼請註明GitHub地址： https://github.com/harryprince。微信號: harryzhustudio 商業使用請聯繫做者。