近年来，基于 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 论文中对三种地理关联查询性能的对比：

拓扑查询的性能(毫秒)

可以看到使用 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 进行图层叠加与拓展。

地理单元

如果不是展示矢量元素，那常见的地理单元有三种选择：

分别是 H3、S2 和 geohash，在R中都提供了对应的实现：

其中，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/...

图层与元素管理

和坐标系所对应的就是底图，参考瓦片地图原理一文，你可以替换各种厂家的底图来适配地图的坐标系：

以 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
商业使用请联系作者。