[原] RStudio Spark/Leaflet 與 GIS 最佳實踐

近年來，基於 Spark 的大數據並行計算方案日漸成熟，在GIS領域有了不少最佳實踐。過去，大多數數據分析師可能都是基於Excel/Hive進行分析工做，可是隨着數據分析架構的成熟，基於 RStudio 和 Spark/Leaflet 的數據分析環境正在變得更加易用和富有生產力。本文將分享 R語言社區最前沿的 Spark/Leaflet 和 GIS 數據處理方法。

在GIS分析領域，咱們最經常使用的數據分析工具就是 PostGIS/QGIS。一方面，PostGIS 自己是一個開源的、單機部署的。這意味着，在物聯網和大數據時代，隨着傳感器和數據量的迅速增加，單機內存、磁盤、網絡等將成爲GIS數據分析的瓶頸。另外一方面，QGIS則是一個基於客戶端非交互式的可視化工具，它沒法快速將本身的想法轉化爲儀表盤生產力。對於互聯網下半場的今天，美團、滴滴、摩拜、高德、Airbnb等擁有豐富GIS數據的公司而言，簡化GIS數據分析門檻變得十分有意義。

通過長期實踐發現，目前基於RStuido 和 Spark/Leaflet 的方案足夠成熟，具備性能穩定、快速移植、低學習成本的特色。

從 PostGIS 到 Spark

在 PostGIS 中，因爲數據庫能夠自建索引，因此在作數據分析時，性能會很是好。可是在數據量很是龐大時，每每也無力招架，而 Spark/Hive 實際上是目前業內處理大數據的最佳選擇和事實標準，如何將數據分析框架遷移到 Spark/Hive 架構成爲不二選擇。

從 PostGIS 到 Spark 能夠分爲兩步：

1. 從 PostGIS 數據庫到 內存計算
2. 從 單機內存 到 Spark 分佈式內存

從 PostGIS 數據庫到 內存計算

空間分析中一般包含幾個核心要素：

1. 空間幾何構建
2. 空間關係計算
3. 空間幾何計算
4. 空間IO
5. 座標系系統

空間幾何轉化

這些核心要素都包含在PostGIS/sf/ESRI 的空間計算函數中，若是你對 PostGIS 很是熟悉，那麼你在學習 R 中的 sf 包就很是簡單不過了，若是你還在入門GIS，那能夠參考這個PostGIS 手冊。 sf包提供了相似 PostGIS的空間計算函數，也一樣包含來上述空間分析的核心要素，很幸運，它們的函數名都是同樣的，可是它能夠直接在R的內存中進行，沒必要要將數據存在數據庫中。簡單、一致、學習成本低，這也是R語言受到GISer垂青的緣由之一。如下面一個例子就能明白它們之間的異同：

PostGIS：

# sc is a postgis connection

sc %>%
tbl(sql("SELECT st_contains(st_polygon(1,1, 1,4, 4,4, 4,1), st_point(2, 3) as res"))

sf：

list(matrix(c(1,1, 1,4, 4,4, 4,1),ncol=2, byrow=TRUE)) %>%
st_polygon() %>% 
st_contains(st_point(c(2,3)))

從 PostGIS 到 sf，分佈式計算的路已經走了一半，起碼已經從 數據庫 到了內存。

從 單機內存 到 Spark 分佈式內存

在 Spark/Hive 生態中，有一種工具叫作 UDF（User Define Function），經過它，能夠自定義一些空間函數來進行GIS分析，而這些工具已經有社區提供，沒必要再造車輪。

引用自 Query the planet: Geospatial big data analytics at Uber 中地理圍欄示意圖

利用 ESRI UDF/GeoSpark 提供的一系列的 ST_* 函數，它能夠將你的 sf 代碼與 sparklyr 很好的結合在一塊兒。下面是一個具體圍欄關聯查詢的例子：

# devtools::install_github("harryprince/geospark")
library(geospark)
library(sparklyr)
config = sparklyr::spark_config()
                                   
sc <- sparklyr::spark_connect(master = "yarn-client", version = "2.2.0",config = config)

# 添加 geospark UDF
register_gis(sc)

# (optional) 添加 ESRI UDF
# DBI::dbGetQuery(sc,"add jar esri-geometry-api-2.1.0.jar") 
# DBI::dbGetQuery(sc,"add jar spatial-sdk-hive-2.1.0-SNAPSHOT.jar")
# DBI::dbGetQuery(sc,"add jar spatial-sdk-json-2.1.0-SNAPSHOT.jar")

# 地理圍欄計算

dplyr::tbl(sc,"ai.financer_tbl") %>% 
    filter(dt=="20180420") %>% 
    mutate(point_a = st_points(c(lat,lng)) %>%
    mutate(polygon_b = st_polygon(c(lat1,lng1,lat2,lng2))) %>%
    transmute(geohash= st_contains(polygon_b, point_a))

它們和 PostGIS/ESRI UDF/GeoSpark/ R 中的 sf 包都繼承了同一套GIS數據處理語言。若是你掌握了 PostGIS 那麼 sf 包和 ESRI UDF/GeoSpark 天然不在話下。
而且 GeoSpark 在ESRI UDF 的基礎上還實現了SQL中的空間Join查詢優化，經過對地理圍欄進行四叉樹的預見索引，總體提高空間Join查詢的性能。

性能對比

下面是 GeoSpark 論文中對三種地理關聯查詢性能的對比：

序號	測試用例	記錄條數
1	SELECT IDCODE FROM zhenlongxiang WHERE ST_Disjoint(geom,ST_GeomFromText(‘POLYGON((517000 1520000,619000 1520000,619000 2530000,517000 2530000,517000 1520000))’));	85,236 rows
2	SELECT fid FROM cyclonepoint WHERE ST_Disjoint(geom,ST_GeomFromText(‘POLYGON((90 3,170 3,170 55,90 55,90 3))’,4326))	60,591 rows

拓撲查詢的性能(毫秒)

序號	PostGIS/PostgreSQL	GeoSpark SQL	ESRI Spatial Framework for Hadoop
1	9631	480	40,784
2	110872	394	64,217

能夠看到使用 GeoSpark SQL 的方式，能夠在數據規模和計算性能上同時超過 PG 和 ESRI UDF，而經過sparklyr的語法糖，能夠完美移植 local的sf代碼到 spark 上。

使用技巧

使用 Spark 的內存並行計算和傳統單機計算最大的差別就是須要利用 惰性計算與異步計算 來優化性能。

首先介紹一下什麼是 惰性計算和異步計算：

• 惰性計算

由於大數據計算中，數據量很是大，若是每次寫代碼，數據就直接開始計算會對計算引擎帶來沒必要要的壓力。

爲了應對大數據計算，spark 和 dplyr 都繼承了 Lazy Computation的思想，也就是把定義計算步驟的圖結構和得到數據結果分開。經過 tbl_df %>% collect() 或者 tbl_df %>% compute() 來實現。

• 異步計算

以一個簡單的例子來比喻就是 咱們早晨起牀的時候，會一邊燒開水，一邊刷牙。燒水並不須要一直守在水壺邊上，燒水的同時作點別的事情，等水燒好來再把水倒出來，這就是簡單的一個異步計算的例子。

在 R 中，咱們經過 future和promises 這兩個包來輔助實現。

具體例子：

pipeline = df %>% # 定義計算圖
head(100)

pipeline 
collect() %>% # 惰性求值，自動優化計算邏輯
View()

library(future)
library(promises)

future(head(df,100)) %...>% # 使用 %...>% 代替 %>% 自動切換爲異步計算模式
    collect() %...>%
    View()

從 Leaflet 到 Dashboard

地理數據在計算以外最重要的一個應用就是可視化技術，而R一直以地圖數據可視化而聞名。國內有 leafletCN/geoChina/REmap 等包，國外有 leaflet/leaflet.extra/mapview/ggmap等，支持各類矢量/柵格數據。

IDE 交互式分析

可視化技術核心有下面幾個要素：

1. 地理單元選擇
2. 座標系選擇
3. BBox/Zoom level
4. 圖層管理
5. 元素管理

mapview、leaflet、R notebook 是地圖交互式分析最好用的三個工具。

在一個 R notebook chunk 中將 sf 對象經過 mapview 直接可視化，再經過 leaflet 進行圖層疊加與拓展。

地理單元

若是不是展現矢量元素，那常見的地理單元有三種選擇：

分別是 H三、S2 和 geohash，在R中都提供了對應的實現：

方法	特色	廠家
H3	美觀、等面積、等形狀、Normlized	Uber
S2	層次豐富、精度高	Google
geohash	成熟實現、各端SDK齊全	-

其中，H3 最適用於機器學習場景，避免了地理單元由於面積、形狀不一致帶來的統計指標的誤差。geohash在北京和深圳的面積偏差就很是嚴重。

座標系

常見的座標系有三種：

分別是 國測局座標系（GCJ-02）、地球座標 (WGS84) 和 百度座標 (BD-09)，在R中的 geoChina包都提供了對應的轉換實現，以亮馬橋地鐵站位置爲例：

# 經緯度轉化
geogcj = data.frame(lat=39.949958,lng=116.46343)
geowgs = geoChina::gcj2wgs(39.949958, 116.46343)
geobd = geoChina::gcj2bd(39.949958,116.46343)

leaflet.mapbox::leafletMapbox(access_token=MAPBOX_TOKEN) %>%
 addCircleMarkers(lat = geogcj$lat,lng = geogcj$lng,color = "red") %>% 
 addCircleMarkers(lat = geowgs$lat,lng = geowgs$lng,color = "green") %>% 
 addCircleMarkers(lat = geobd$lat,lng = geobd$lng,color = "blue") %>% 
 # leafletCN::amap()
 leaflet.mapbox::addMapboxTileLayer(mapbox.classicStyleIds$dark)
REmap::remapB(markPointData = data.frame(a=c("gcj","wgs","bd"),
color=c("red","green","blue")),color = "Blue",geoData = geogcj %>% 
rbind(geowgs) %>%
rbind(geobd) %>% 
select(lon = lng,lat= lat) %>% 
cbind(city=c("gcj","wgs","bd")))

圖中，紅點表示 高德座標系；藍點表示百度座標系；綠點表示 Google座標系。
如圖所示，在國內 高德地圖（左一白）底圖中，紅點能夠準確表示 亮馬橋地鐵站；
在國內百度地圖（中藍）底圖中，藍點能夠準確表示亮馬橋地鐵站位置。
在國際Google地圖（右黑）底圖中，綠點能夠準確表示亮馬橋地鐵站位置。

對於初學者最容易犯的一個錯誤就是將座標系混淆使用，這一般會帶來上百米的精度偏差。雖然，各個廠家標準不一，但好在各類語言但開源轉化方案也不少，好比：https://github.com/googollee/...

圖層與元素管理

和座標系所對應的就是底圖，參考瓦片地圖原理一文，你能夠替換各類廠家的底圖來適配地圖的座標系：

地圖商	瓦片編碼	圖層	底圖URI
高德地圖	谷歌XYZ	道路	http://webrd02.is.autonavi.co...
高德地圖	谷歌XYZ	衛星	http://webst04.is.autonavi.co...
谷歌地圖	谷歌XYZ	道路	http://mt2.google.cn/vt/lyrs=...
谷歌地圖	谷歌XYZ	衛星	http://mt2.google.cn/vt/lyrs=...
谷歌地圖	谷歌XYZ	地形	http://mt0.google.cn/vt/lyrs=...
OpenStreetMap	谷歌XYZ	道路	http://a.tile.openstreetmap.o...
騰訊地圖	TMS	道路	http://rt1.map.gtimg.com/real...
Bing地圖	QuadTree	道路	http://r1.tiles.ditu.live.com...
百度地圖	百度XYZ	道路	http://online4.map.bdimg.com/...;styles=pl&scaler=1&udt=20170406
百度地圖	百度XYZ	交通	http://its.map.baidu.com:8002/traffic/TrafficTileService?level=19&x=99052&y=20189&time=1373790856265&label=web2D&;v=017

以 leaflet 官方文檔爲例：

outline <- quakes[chull(quakes$long, quakes$lat),]

leaflet(quakes) %>%
  # Base groups
  addTiles(group = "OSM (default)") %>%
  addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
  addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
  # Overlay groups
  addCircles(~long, ~lat, ~10^mag/5, stroke = F, group = "Quakes") %>%
  addPolygons(data = outline, lng = ~long, lat = ~lat,
    fill = F, weight = 2, color = "#FFFFCC", group = "Outline") %>%
  # Layers control
  addLayersControl(
    baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
    overlayGroups = c("Quakes", "Outline"),
    options = layersControlOptions(collapsed = FALSE)
  )

底圖、點、多邊形等元素在圖層上都獲得統一的管理，更多 leaflet 及其插件的操做能夠參考 時空維度挖掘之 leaflet。

儀表盤構建

在咱們寫完一段臨時的分析代碼，只要稍做排版上的改進，咱們就能夠將它發佈爲一個 dashboard，這極大釋放地理分析的生產力，這主要得益於 Rmarkdown、flexdashboard、shiny@async 這三個工具。

---
title: "FinanceR"
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    vertical_layout: fill
runtime: shiny
---

/```{r setup, include=FALSE}
library(flexdashboard)
/```

Column {data-width=650}
-----------------------------------------------------------------------

### Chart A

/```{r}
library(dplyr)
library(leaflet)
outline <- quakes[chull(quakes$long, quakes$lat),]

leaflet(quakes) %>%
  # Base groups
  addTiles(group = "OSM (default)") %>%
  addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
  addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
  # Overlay groups
  addCircles(~long, ~lat, ~10^mag/5, stroke = F, group = "Quakes") %>%
  addPolygons(data = outline, lng = ~long, lat = ~lat,
    fill = F, weight = 2, color = "#FFFFCC", group = "Outline") %>%
  # Layers control
  addLayersControl(
    baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
    overlayGroups = c("Quakes", "Outline"),
    options = layersControlOptions(collapsed = FALSE)
  )
/```

在原來的分析圖表基礎上，稍微作一些 markdown 樣式的配置就能夠生成一個dashboard而且發佈。更多樣式配置的方法能夠參考 打造數據產品的快速原型：如何使用 flexdashboard 製做dashboard

總結

1. 數據操做徹底複用 dplyr 接口，無縫切換 PostGIS/MySQL到Spark。
2. R本地代碼可直接經過 sparklyr::spark_apply 分佈式執行，不須要手動管理R包依賴，目前算法分發最簡單模式。
3. ESRI UDF 直接遷移本地 sf 空間計算函數到 Spark上，因爲 dplyr 直接兼容 hive UDF，能夠直接複用，經過 dbplyr::sql_render 生成對應SQL 在PostGIS/Hive 可繼續使用。
4. Rmarkdown + sparklyr 下降本地調試和調度部署切換的成本，經過 params 區分 yarn-client 和 yarn-cluster模式。
5. 經過 future + promises 實現 spark 無縫遷移到異步計算模式
6. mapview::mapview + leaflet + shiny@async + flexdashboard 可視化方案，默認配置足夠簡單，調用快速，拓展性強。
7. 經過 Rcpp/sparklyr/reticulate/V8 能夠和 C++/JVM/Python/JS等其餘語言通信，拓展性強。

參考資料

• https://rstudio.github.io/lea...
• http://spark.rstudio.com/
• 瓦片地圖原理
• 打造數據產品的快速原型：如何使用 flexdashboard 製做dashboard
• ESRI UDF
• Query the planet: Geospatial big data analytics at Uber
• 互聯網在線地圖平臺對比分析

推薦產品

• ESRI
• ElasticSearch

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