数据中心运维管理 2019-06-27
近年来,基于 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/Hive 其实是目前业内处理大数据的最佳选择和事实标准,如何将数据分析框架迁移到 Spark/Hive 架构成为不二选择。
从 PostGIS 到 Spark 可以分为两步:
空间分析中通常包含几个核心要素:
这些核心要素都包含在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/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()
地理数据在计算之外最重要的一个应用就是可视化技术,而R一直以地图数据可视化而闻名。国内有 leafletCN
/geoChina
/REmap
等包,国外有 leaflet
/leaflet.extra
/mapview
/ggmap
等,支持各种矢量/栅格数据。
可视化技术核心有下面几个要素:
mapview
、leaflet
、R notebook
是地图交互式分析最好用的三个工具。
在一个 R notebook chunk 中将 sf 对象通过 mapview 直接可视化,再通过 leaflet 进行图层叠加与拓展。
如果不是展示矢量元素,那常见的地理单元有三种选择:
分别是 H3、S2 和 geohash,在R中都提供了对应的实现:
方法 | 特点 | 厂家 |
---|---|---|
H3 | 美观、等面积、等形状、Normlized | Uber |
S2 | 层次丰富、精度高 | |
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
作为分享主义者(sharism),本人所有互联网发布的图文均遵从CC版权,转载请保留作者信息并注明作者 Harry Zhu 的 FinanceR专栏:https://segmentfault.com/blog...,如果涉及源代码请注明GitHub地址:https://github.com/harryprince。微信号: harryzhustudio
商业使用请联系作者。